Rev 855 | Rev 901 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
671 | werner | 2 | /******************************************************************************************** |
3 | ** iLand - an individual based forest landscape and disturbance model |
||
4 | ** http://iland.boku.ac.at |
||
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl |
||
6 | ** |
||
7 | ** This program is free software: you can redistribute it and/or modify |
||
8 | ** it under the terms of the GNU General Public License as published by |
||
9 | ** the Free Software Foundation, either version 3 of the License, or |
||
10 | ** (at your option) any later version. |
||
11 | ** |
||
12 | ** This program is distributed in the hope that it will be useful, |
||
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
||
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
||
15 | ** GNU General Public License for more details. |
||
16 | ** |
||
17 | ** You should have received a copy of the GNU General Public License |
||
18 | ** along with this program. If not, see <http://www.gnu.org/licenses/>. |
||
19 | ********************************************************************************************/ |
||
20 | |||
595 | werner | 21 | #include "sapling.h" |
22 | #include "model.h" |
||
23 | #include "species.h" |
||
24 | #include "resourceunit.h" |
||
25 | #include "resourceunitspecies.h" |
||
26 | #include "tree.h" |
||
27 | |||
28 | /** @class Sapling |
||
697 | werner | 29 | @ingroup core |
595 | werner | 30 | Sapling stores saplings per species and resource unit and computes sapling growth (before recruitment). |
31 | http://iland.boku.ac.at/sapling+growth+and+competition |
||
32 | Saplings are established in a separate step (@sa Regeneration). If sapling reach a height of 4m, they are recruited and become "real" iLand-trees. |
||
33 | Within the regeneration layer, a cohort-approach is applied. |
||
34 | |||
35 | */ |
||
36 | |||
37 | double Sapling::mRecruitmentVariation = 0.1; // +/- 10% |
||
38 | |||
39 | Sapling::Sapling() |
||
40 | { |
||
41 | mRUS = 0; |
||
42 | clearStatistics(); |
||
663 | werner | 43 | mAdded = 0; |
595 | werner | 44 | } |
45 | |||
663 | werner | 46 | // reset statistics, called at newYear |
47 | void Sapling::clearStatistics() |
||
48 | { |
||
49 | // mAdded: removed |
||
50 | mRecruited=mDied=mLiving=0; |
||
51 | mSumDbhDied=0.; |
||
52 | mAvgHeight=0.; |
||
53 | mAvgAge=0.; |
||
54 | mAvgDeltaHPot=mAvgHRealized=0.; |
||
55 | } |
||
56 | |||
595 | werner | 57 | /// get the *represented* (Reineke's Law) number of trees (N/ha) |
58 | double Sapling::livingStemNumber(double &rAvgDbh, double &rAvgHeight, double &rAvgAge) const |
||
59 | { |
||
60 | double total = 0.; |
||
61 | double dbh_sum = 0.; |
||
62 | double h_sum = 0.; |
||
63 | double age_sum = 0.; |
||
64 | const SaplingGrowthParameters &p = mRUS->species()->saplingGrowthParameters(); |
||
65 | for (QVector<SaplingTree>::const_iterator it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
66 | float dbh = it->height / p.hdSapling * 100.f; |
||
67 | if (dbh<1.) // minimum size: 1cm |
||
68 | continue; |
||
69 | double n = p.representedStemNumber(dbh); // one cohort on the pixel represents that number of trees |
||
70 | dbh_sum += n*dbh; |
||
71 | h_sum += n*it->height; |
||
72 | age_sum += n*it->age.age; |
||
73 | total += n; |
||
74 | } |
||
75 | if (total>0.) { |
||
76 | dbh_sum /= total; |
||
77 | h_sum /= total; |
||
78 | age_sum /= total; |
||
79 | } |
||
80 | rAvgDbh = dbh_sum; |
||
81 | rAvgHeight = h_sum; |
||
82 | rAvgAge = age_sum; |
||
83 | return total; |
||
84 | } |
||
85 | |||
86 | /// maintenance function to clear dead/recruited saplings from storage |
||
87 | void Sapling::cleanupStorage() |
||
88 | { |
||
89 | QVector<SaplingTree>::iterator forw=mSaplingTrees.begin(); |
||
90 | QVector<SaplingTree>::iterator back; |
||
91 | |||
92 | // seek last valid |
||
93 | for (back=mSaplingTrees.end()-1; back>=mSaplingTrees.begin(); --back) |
||
94 | if ((*back).isValid()) |
||
95 | break; |
||
96 | |||
97 | if (back<mSaplingTrees.begin()) { |
||
98 | mSaplingTrees.clear(); // no valid trees available |
||
99 | return; |
||
100 | } |
||
101 | |||
102 | while (forw < back) { |
||
103 | if (!(*forw).isValid()) { |
||
104 | *forw = *back; // copy (fill gap) |
||
105 | while (back>forw) // seek next valid |
||
106 | if ((*--back).isValid()) |
||
107 | break; |
||
108 | } |
||
109 | ++forw; |
||
110 | } |
||
111 | if (back != mSaplingTrees.end()-1) { |
||
112 | // free resources... |
||
113 | mSaplingTrees.erase(back+1, mSaplingTrees.end()); |
||
114 | } |
||
115 | } |
||
116 | |||
117 | // not a very good way of checking if sapling is present |
||
118 | // maybe better: use also a (local) maximum sapling height grid |
||
119 | // maybe better: use a bitset: |
||
120 | // position: index of pixel on LIF (absolute index) |
||
121 | bool Sapling::hasSapling(const QPoint &position) const |
||
122 | { |
||
123 | const QPoint &offset = mRUS->ru()->cornerPointOffset(); |
||
124 | int index = (position.x()- offset.x())*cPxPerRU + (position.y() - offset.y()); |
||
802 | werner | 125 | if (index<0) |
126 | qDebug() << "Sapling error"; |
||
595 | werner | 127 | return mSapBitset[index]; |
128 | /* |
||
129 | float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y()); |
||
130 | QVector<SaplingTree>::const_iterator it; |
||
131 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
132 | if (it->pixel==target) |
||
133 | return true; |
||
134 | } |
||
135 | return false; |
||
136 | */ |
||
137 | } |
||
138 | |||
600 | werner | 139 | /// retrieve the height of the sapling at the location 'position' (given in LIF-coordinates) |
140 | /// this is quite expensive and only done for initialization |
||
141 | double Sapling::heightAt(const QPoint &position) const |
||
142 | { |
||
143 | if (!hasSapling(position)) |
||
144 | return 0.; |
||
145 | // ok, we'll have to search through all saplings |
||
146 | QVector<SaplingTree>::const_iterator it; |
||
147 | float *lif_ptr = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y()); |
||
148 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
149 | if (it->isValid() && it->pixel == lif_ptr) |
||
150 | return it->height; |
||
151 | } |
||
152 | return 0.; |
||
595 | werner | 153 | |
600 | werner | 154 | } |
155 | |||
156 | |||
695 | werner | 157 | void Sapling::setBit(const QPoint &pos_index) |
158 | { |
||
159 | int index = (pos_index.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(pos_index.y() - mRUS->ru()->cornerPointOffset().y()); |
||
160 | mSapBitset.set(index,true); // set bit: now there is a sapling there |
||
161 | } |
||
162 | |||
163 | /// a a sapling at given position (index on the LIF grid, i.e. 2x2m) |
||
595 | werner | 164 | void Sapling::addSapling(const QPoint &pos_lif) |
165 | { |
||
166 | // adds a sapling... |
||
167 | mSaplingTrees.push_back(SaplingTree()); |
||
168 | SaplingTree &t = mSaplingTrees.back(); |
||
780 | werner | 169 | t.height = 0.05f; // start with 5cm height |
595 | werner | 170 | Grid<float> &lif_map = *GlobalSettings::instance()->model()->grid(); |
171 | t.pixel = lif_map.ptr(pos_lif.x(), pos_lif.y()); |
||
695 | werner | 172 | setBit(pos_lif); |
595 | werner | 173 | mAdded++; |
174 | } |
||
175 | |||
176 | /// clear saplings on a given position (after recruitment) |
||
177 | void Sapling::clearSaplings(const QPoint &position) |
||
178 | { |
||
179 | float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y()); |
||
180 | QVector<SaplingTree>::const_iterator it; |
||
181 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
182 | if (it->pixel==target) { |
||
183 | // trick: use a const iterator to avoid a deep copy of the vector; then do an ugly const_cast to actually write the data |
||
662 | werner | 184 | //const SaplingTree &t = *it; |
185 | //const_cast<SaplingTree&>(t).pixel=0; |
||
186 | clearSapling(const_cast<SaplingTree&>(*it), false); // kill sapling and move carbon to soil |
||
595 | werner | 187 | } |
188 | } |
||
189 | int index = (position.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(position.y() - mRUS->ru()->cornerPointOffset().y()); |
||
190 | mSapBitset.set(index,false); // clear bit: now there is no sapling on this position |
||
191 | |||
192 | } |
||
193 | |||
662 | werner | 194 | /// clear saplings within a given rectangle |
195 | void Sapling::clearSaplings(const QRectF &rectangle, const bool remove_biomass) |
||
196 | { |
||
197 | QVector<SaplingTree>::const_iterator it; |
||
198 | FloatGrid *grid = GlobalSettings::instance()->model()->grid(); |
||
199 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
200 | if (rectangle.contains(grid->cellCenterPoint(grid->indexOf(it->pixel)))) { |
||
201 | clearSapling(const_cast<SaplingTree&>(*it), remove_biomass); |
||
202 | } |
||
203 | } |
||
204 | } |
||
205 | |||
206 | void Sapling::clearSapling(SaplingTree &tree, const bool remove) |
||
207 | { |
||
208 | tree.pixel=0; |
||
209 | if (!remove) { |
||
210 | // killing of saplings: |
||
211 | // if remove=false, then remember dbh/number of trees (used later in calculateGrowth() to estimate carbon flow) |
||
212 | mDied++; |
||
213 | mSumDbhDied+=tree.height / mRUS->species()->saplingGrowthParameters().hdSapling * 100.; |
||
214 | } |
||
215 | |||
216 | } |
||
217 | |||
595 | werner | 218 | /// growth function for an indivudal sapling. |
219 | /// returns true, if sapling survives, false if sapling dies or is recruited to iLand. |
||
220 | /// see also http://iland.boku.ac.at/recruitment |
||
221 | bool Sapling::growSapling(SaplingTree &tree, const double f_env_yr, Species* species) |
||
222 | { |
||
223 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(tree.pixel); |
||
224 | |||
225 | // (1) calculate height growth potential for the tree (uses linerization of expressions...) |
||
802 | werner | 226 | double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height); // TODO check if this can be source of crashes (race condition) |
595 | werner | 227 | double delta_h_pot = h_pot - tree.height; |
228 | |||
229 | // (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel. |
||
230 | double lif_value = *tree.pixel; |
||
231 | double h_height_grid = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(p.x()/cPxPerHeight, p.y()/cPxPerHeight).height; |
||
232 | if (h_height_grid==0) |
||
233 | throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(p.x()).arg(p.y())); |
||
234 | |||
816 | werner | 235 | double rel_height = tree.height / h_height_grid; |
236 | |||
595 | werner | 237 | double lif_corrected = mRUS->species()->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height |
683 | werner | 238 | |
595 | werner | 239 | double lr = mRUS->species()->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index) |
240 | |||
241 | double delta_h_factor = f_env_yr * lr; // relative growth |
||
242 | |||
243 | if (h_pot<0. || delta_h_pot<0. || lif_corrected<0. || lif_corrected>1. || delta_h_factor<0. || delta_h_factor>1. ) |
||
244 | qDebug() << "invalid values in Sapling::growSapling"; |
||
245 | |||
246 | // check mortality of saplings |
||
247 | if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) { |
||
248 | tree.age.stress_years++; |
||
249 | if (tree.age.stress_years > species->saplingGrowthParameters().maxStressYears) { |
||
250 | // sapling dies... |
||
662 | werner | 251 | clearSapling(tree, false); // false: put carbon to the soil |
595 | werner | 252 | return false; |
253 | } |
||
254 | } else { |
||
255 | tree.age.stress_years=0; // reset stress counter |
||
256 | } |
||
257 | DBG_IF(delta_h_pot*delta_h_factor < 0.f || delta_h_pot*delta_h_factor > 2., "Sapling::growSapling", "inplausible height growth."); |
||
258 | |||
259 | // grow |
||
260 | tree.height += delta_h_pot * delta_h_factor; |
||
261 | tree.age.age++; // increase age of sapling by 1 |
||
262 | |||
263 | // recruitment? |
||
264 | if (tree.height > 4.f) { |
||
265 | mRecruited++; |
||
266 | |||
267 | ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru()); |
||
268 | float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f; |
||
269 | // the number of trees to create (result is in trees per pixel) |
||
270 | double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh); |
||
271 | int to_establish = (int) n_trees; |
||
863 | werner | 272 | |
595 | werner | 273 | // if n_trees is not an integer, choose randomly if we should add a tree. |
274 | // e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%. |
||
707 | werner | 275 | if (drandom() < (n_trees-to_establish) || to_establish==0) |
595 | werner | 276 | to_establish++; |
277 | |||
278 | // add a new tree |
||
279 | for (int i=0;i<to_establish;i++) { |
||
280 | Tree &bigtree = ru->newTree(); |
||
281 | bigtree.setPosition(p); |
||
282 | // add variation: add +/-10% to dbh and *independently* to height. |
||
707 | werner | 283 | bigtree.setDbh(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
284 | bigtree.setHeight(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
||
595 | werner | 285 | bigtree.setSpecies( species ); |
286 | bigtree.setAge(tree.age.age,tree.height); |
||
287 | bigtree.setRU(ru); |
||
288 | bigtree.setup(); |
||
855 | werner | 289 | const Tree *t = &bigtree; |
290 | mRUS->statistics().add(t, 0); // count the newly created trees already in the stats |
||
595 | werner | 291 | } |
292 | // clear all regeneration from this pixel (including this tree) |
||
662 | werner | 293 | clearSapling(tree, true); // remove this tree (but do not move biomass to soil) |
294 | ru->clearSaplings(p); // remove all other saplings on the same pixel |
||
855 | werner | 295 | |
595 | werner | 296 | return false; |
297 | } |
||
298 | // book keeping (only for survivors) |
||
299 | mLiving++; |
||
300 | mAvgHeight+=tree.height; |
||
301 | mAvgAge+=tree.age.age; |
||
302 | mAvgDeltaHPot+=delta_h_pot; |
||
303 | mAvgHRealized += delta_h_pot * delta_h_factor; |
||
304 | return true; |
||
305 | |||
306 | } |
||
307 | |||
608 | werner | 308 | |
662 | werner | 309 | /** main growth function for saplings. |
310 | Statistics are cleared at the beginning of the year. |
||
311 | */ |
||
595 | werner | 312 | void Sapling::calculateGrowth() |
313 | { |
||
314 | Q_ASSERT(mRUS); |
||
821 | werner | 315 | if (mSaplingTrees.count()==0) |
595 | werner | 316 | return; |
317 | |||
318 | ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru() ); |
||
319 | Species *species = const_cast<Species*>(mRUS->species()); |
||
320 | |||
321 | // calculate necessary growth modifier (this is done only once per year) |
||
322 | mRUS->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration |
||
323 | double f_env_yr = mRUS->prod3PG().fEnvYear(); |
||
324 | |||
325 | mLiving=0; |
||
326 | QVector<SaplingTree>::const_iterator it; |
||
327 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
328 | const SaplingTree &tree = *it; |
||
329 | if (tree.height<0) |
||
330 | qDebug() << "Sapling::calculateGrowth(): h<0"; |
||
821 | werner | 331 | // if sapling is still living check execute growth routine |
595 | werner | 332 | if (tree.isValid()) { |
821 | werner | 333 | // growing (increases mLiving if tree did not die, mDied otherwise) |
595 | werner | 334 | if (growSapling(const_cast<SaplingTree&>(tree), f_env_yr, species)) { |
335 | // set the sapling height to the maximum value on the current pixel |
||
336 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(tree.pixel); |
||
337 | ru->setMaxSaplingHeightAt(p,tree.height); |
||
338 | } |
||
339 | } |
||
340 | } |
||
341 | if (mLiving) { |
||
342 | mAvgHeight /= double(mLiving); |
||
343 | mAvgAge /= double(mLiving); |
||
344 | mAvgDeltaHPot /= double(mLiving); |
||
345 | mAvgHRealized /= double(mLiving); |
||
346 | } |
||
347 | // calculate carbon balance |
||
608 | werner | 348 | CNPair old_state = mCarbonLiving; |
595 | werner | 349 | mCarbonLiving.clear(); |
608 | werner | 350 | |
595 | werner | 351 | CNPair dead_wood, dead_fine; // pools for mortality |
352 | // average dbh |
||
353 | if (mLiving) { |
||
354 | // calculate the avg dbh number of stems |
||
355 | double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.; |
||
356 | double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh ); |
||
357 | // woody parts: stem, branchse and coarse roots |
||
358 | double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh); |
||
359 | double foliage = species->biomassFoliage(avg_dbh); |
||
360 | double fineroot = foliage*species->finerootFoliageRatio(); |
||
361 | |||
362 | mCarbonLiving.addBiomass( woody_bm*n, species->cnWood() ); |
||
363 | mCarbonLiving.addBiomass( foliage*n, species->cnFoliage() ); |
||
364 | mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot() ); |
||
608 | werner | 365 | |
595 | werner | 366 | // turnover |
802 | werner | 367 | if (mRUS->ru()->snag()) |
368 | mRUS->ru()->snag()->addTurnoverLitter(species, foliage*species->turnoverLeaf(), fineroot*species->turnoverRoot()); |
||
595 | werner | 369 | |
370 | // calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers |
||
371 | // from Reinekes formula. |
||
372 | // |
||
373 | if (avg_dbh>1.) { |
||
374 | double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.; |
||
375 | double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) ); |
||
376 | if (n<n_before) { |
||
377 | dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() ); |
||
378 | dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage() ); |
||
379 | dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot() ); |
||
380 | } |
||
381 | } |
||
382 | |||
383 | } |
||
384 | if (mDied) { |
||
385 | double avg_dbh_dead = mSumDbhDied / double(mDied); |
||
386 | double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead ); |
||
387 | // woody parts: stem, branchse and coarse roots |
||
388 | |||
389 | dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood() ); |
||
390 | double foliage = species->biomassFoliage(avg_dbh_dead)*n; |
||
391 | |||
392 | dead_fine.addBiomass( foliage, species->cnFoliage() ); |
||
393 | dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot() ); |
||
394 | } |
||
395 | if (!dead_wood.isEmpty() || !dead_fine.isEmpty()) |
||
802 | werner | 396 | if (mRUS->ru()->snag()) |
397 | mRUS->ru()->snag()->addToSoil(species, dead_wood, dead_fine); |
||
595 | werner | 398 | |
608 | werner | 399 | // calculate net growth: |
400 | // delta of stocks |
||
625 | werner | 401 | mCarbonGain = mCarbonLiving + dead_fine + dead_wood - old_state; |
402 | if (mCarbonGain.C < 0) |
||
403 | mCarbonGain.clear(); |
||
608 | werner | 404 | |
595 | werner | 405 | if (mSaplingTrees.count() > mLiving*1.3) |
406 | cleanupStorage(); |
||
407 | |||
408 | mRUS->statistics().add(this); |
||
615 | werner | 409 | GlobalSettings::instance()->systemStatistics()->saplingCount+=mLiving; |
410 | GlobalSettings::instance()->systemStatistics()->newSaplings+=mAdded; |
||
663 | werner | 411 | mAdded = 0; // reset |
412 | |||
595 | werner | 413 | //qDebug() << ru->index() << species->id()<< ": (living/avg.height):" << mLiving << mAvgHeight; |
414 | } |
||
415 | |||
416 | /// fill a grid with the maximum height of saplings per pixel (2x2m). |
||
417 | /// this function is used for visualization only |
||
418 | void Sapling::fillMaxHeightGrid(Grid<float> &grid) const |
||
419 | { |
||
420 | QVector<SaplingTree>::const_iterator it; |
||
421 | for (it = mSaplingTrees.begin(); it!=mSaplingTrees.end(); ++it) { |
||
422 | if (it->isValid()) { |
||
423 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(it->pixel); |
||
424 | if (grid.valueAtIndex(p)<it->height) |
||
425 | grid.valueAtIndex(p) = it->height; |
||
426 | } |
||
427 | } |
||
428 | |||
429 | } |
||
600 | werner | 430 | |
608 | werner | 431 | |
432 | |||
695 | werner | 433 |