Rev 593 | Rev 600 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
595 | werner | 2 | #include "sapling.h" |
3 | #include "model.h" |
||
4 | #include "species.h" |
||
5 | #include "resourceunit.h" |
||
6 | #include "resourceunitspecies.h" |
||
7 | #include "tree.h" |
||
8 | |||
9 | /** @class Sapling |
||
10 | Sapling stores saplings per species and resource unit and computes sapling growth (before recruitment). |
||
11 | http://iland.boku.ac.at/sapling+growth+and+competition |
||
12 | Saplings are established in a separate step (@sa Regeneration). If sapling reach a height of 4m, they are recruited and become "real" iLand-trees. |
||
13 | Within the regeneration layer, a cohort-approach is applied. |
||
14 | |||
15 | */ |
||
16 | |||
17 | double Sapling::mRecruitmentVariation = 0.1; // +/- 10% |
||
18 | |||
19 | Sapling::Sapling() |
||
20 | { |
||
21 | mRUS = 0; |
||
22 | clearStatistics(); |
||
23 | } |
||
24 | |||
25 | /// get the *represented* (Reineke's Law) number of trees (N/ha) |
||
26 | double Sapling::livingStemNumber(double &rAvgDbh, double &rAvgHeight, double &rAvgAge) const |
||
27 | { |
||
28 | double total = 0.; |
||
29 | double dbh_sum = 0.; |
||
30 | double h_sum = 0.; |
||
31 | double age_sum = 0.; |
||
32 | const SaplingGrowthParameters &p = mRUS->species()->saplingGrowthParameters(); |
||
33 | for (QVector<SaplingTree>::const_iterator it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
34 | float dbh = it->height / p.hdSapling * 100.f; |
||
35 | if (dbh<1.) // minimum size: 1cm |
||
36 | continue; |
||
37 | double n = p.representedStemNumber(dbh); // one cohort on the pixel represents that number of trees |
||
38 | dbh_sum += n*dbh; |
||
39 | h_sum += n*it->height; |
||
40 | age_sum += n*it->age.age; |
||
41 | total += n; |
||
42 | } |
||
43 | if (total>0.) { |
||
44 | dbh_sum /= total; |
||
45 | h_sum /= total; |
||
46 | age_sum /= total; |
||
47 | } |
||
48 | rAvgDbh = dbh_sum; |
||
49 | rAvgHeight = h_sum; |
||
50 | rAvgAge = age_sum; |
||
51 | return total; |
||
52 | } |
||
53 | |||
54 | /// maintenance function to clear dead/recruited saplings from storage |
||
55 | void Sapling::cleanupStorage() |
||
56 | { |
||
57 | QVector<SaplingTree>::iterator forw=mSaplingTrees.begin(); |
||
58 | QVector<SaplingTree>::iterator back; |
||
59 | |||
60 | // seek last valid |
||
61 | for (back=mSaplingTrees.end()-1; back>=mSaplingTrees.begin(); --back) |
||
62 | if ((*back).isValid()) |
||
63 | break; |
||
64 | |||
65 | if (back<mSaplingTrees.begin()) { |
||
66 | mSaplingTrees.clear(); // no valid trees available |
||
67 | return; |
||
68 | } |
||
69 | |||
70 | while (forw < back) { |
||
71 | if (!(*forw).isValid()) { |
||
72 | *forw = *back; // copy (fill gap) |
||
73 | while (back>forw) // seek next valid |
||
74 | if ((*--back).isValid()) |
||
75 | break; |
||
76 | } |
||
77 | ++forw; |
||
78 | } |
||
79 | if (back != mSaplingTrees.end()-1) { |
||
80 | // free resources... |
||
81 | mSaplingTrees.erase(back+1, mSaplingTrees.end()); |
||
82 | } |
||
83 | } |
||
84 | |||
85 | // not a very good way of checking if sapling is present |
||
86 | // maybe better: use also a (local) maximum sapling height grid |
||
87 | // maybe better: use a bitset: |
||
88 | // position: index of pixel on LIF (absolute index) |
||
89 | bool Sapling::hasSapling(const QPoint &position) const |
||
90 | { |
||
91 | const QPoint &offset = mRUS->ru()->cornerPointOffset(); |
||
92 | int index = (position.x()- offset.x())*cPxPerRU + (position.y() - offset.y()); |
||
93 | return mSapBitset[index]; |
||
94 | /* |
||
95 | float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y()); |
||
96 | QVector<SaplingTree>::const_iterator it; |
||
97 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
98 | if (it->pixel==target) |
||
99 | return true; |
||
100 | } |
||
101 | return false; |
||
102 | */ |
||
103 | } |
||
104 | |||
105 | |||
106 | void Sapling::addSapling(const QPoint &pos_lif) |
||
107 | { |
||
108 | // adds a sapling... |
||
109 | mSaplingTrees.push_back(SaplingTree()); |
||
110 | SaplingTree &t = mSaplingTrees.back(); |
||
111 | t.height = 0.05; // start with 5cm height |
||
112 | Grid<float> &lif_map = *GlobalSettings::instance()->model()->grid(); |
||
113 | t.pixel = lif_map.ptr(pos_lif.x(), pos_lif.y()); |
||
114 | int index = (pos_lif.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(pos_lif.y() - mRUS->ru()->cornerPointOffset().y()); |
||
115 | mSapBitset.set(index,true); // set bit: now there is a sapling there |
||
116 | mAdded++; |
||
117 | } |
||
118 | |||
119 | /// clear saplings on a given position (after recruitment) |
||
120 | void Sapling::clearSaplings(const QPoint &position) |
||
121 | { |
||
122 | float *target = GlobalSettings::instance()->model()->grid()->ptr(position.x(), position.y()); |
||
123 | QVector<SaplingTree>::const_iterator it; |
||
124 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
125 | if (it->pixel==target) { |
||
126 | // trick: use a const iterator to avoid a deep copy of the vector; then do an ugly const_cast to actually write the data |
||
127 | const SaplingTree &t = *it; |
||
128 | const_cast<SaplingTree&>(t).pixel=0; |
||
129 | } |
||
130 | } |
||
131 | int index = (position.x() - mRUS->ru()->cornerPointOffset().x()) * cPxPerRU +(position.y() - mRUS->ru()->cornerPointOffset().y()); |
||
132 | mSapBitset.set(index,false); // clear bit: now there is no sapling on this position |
||
133 | |||
134 | } |
||
135 | |||
136 | /// growth function for an indivudal sapling. |
||
137 | /// returns true, if sapling survives, false if sapling dies or is recruited to iLand. |
||
138 | /// see also http://iland.boku.ac.at/recruitment |
||
139 | bool Sapling::growSapling(SaplingTree &tree, const double f_env_yr, Species* species) |
||
140 | { |
||
141 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(tree.pixel); |
||
142 | |||
143 | // (1) calculate height growth potential for the tree (uses linerization of expressions...) |
||
144 | double h_pot = species->saplingGrowthParameters().heightGrowthPotential.calculate(tree.height); |
||
145 | double delta_h_pot = h_pot - tree.height; |
||
146 | |||
147 | // (2) reduce height growth potential with species growth response f_env_yr and with light state (i.e. LIF-value) of home-pixel. |
||
148 | double lif_value = *tree.pixel; |
||
149 | double h_height_grid = GlobalSettings::instance()->model()->heightGrid()->valueAtIndex(p.x()/cPxPerHeight, p.y()/cPxPerHeight).height; |
||
150 | if (h_height_grid==0) |
||
151 | throw IException(QString("growSapling: height grid at %1/%2 has value 0").arg(p.x()).arg(p.y())); |
||
152 | double rel_height = 4. / h_height_grid; |
||
153 | |||
154 | double lif_corrected = mRUS->species()->speciesSet()->LRIcorrection(lif_value, rel_height); // correction based on height |
||
155 | // Note: difference to trees: no "LRIcorrection" |
||
156 | double lr = mRUS->species()->lightResponse(lif_corrected); // species specific light response (LUI, light utilization index) |
||
157 | |||
158 | double delta_h_factor = f_env_yr * lr; // relative growth |
||
159 | |||
160 | if (h_pot<0. || delta_h_pot<0. || lif_corrected<0. || lif_corrected>1. || delta_h_factor<0. || delta_h_factor>1. ) |
||
161 | qDebug() << "invalid values in Sapling::growSapling"; |
||
162 | |||
163 | // check mortality of saplings |
||
164 | if (delta_h_factor < species->saplingGrowthParameters().stressThreshold) { |
||
165 | tree.age.stress_years++; |
||
166 | if (tree.age.stress_years > species->saplingGrowthParameters().maxStressYears) { |
||
167 | // sapling dies... |
||
168 | tree.pixel=0; |
||
169 | mDied++; |
||
170 | mSumDbhDied+=tree.height / mRUS->species()->saplingGrowthParameters().hdSapling * 100.; |
||
171 | return false; |
||
172 | } |
||
173 | } else { |
||
174 | tree.age.stress_years=0; // reset stress counter |
||
175 | } |
||
176 | DBG_IF(delta_h_pot*delta_h_factor < 0.f || delta_h_pot*delta_h_factor > 2., "Sapling::growSapling", "inplausible height growth."); |
||
177 | |||
178 | // grow |
||
179 | tree.height += delta_h_pot * delta_h_factor; |
||
180 | tree.age.age++; // increase age of sapling by 1 |
||
181 | |||
182 | // recruitment? |
||
183 | if (tree.height > 4.f) { |
||
184 | mRecruited++; |
||
185 | |||
186 | ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru()); |
||
187 | float dbh = tree.height / species->saplingGrowthParameters().hdSapling * 100.f; |
||
188 | // the number of trees to create (result is in trees per pixel) |
||
189 | double n_trees = species->saplingGrowthParameters().representedStemNumber(dbh); |
||
190 | int to_establish = (int) n_trees; |
||
191 | // if n_trees is not an integer, choose randomly if we should add a tree. |
||
192 | // e.g.: n_trees = 2.3 -> add 2 trees with 70% probability, and add 3 trees with p=30%. |
||
193 | if (drandom() < (n_trees-to_establish) || to_establish==0) |
||
194 | to_establish++; |
||
195 | |||
196 | // add a new tree |
||
197 | for (int i=0;i<to_establish;i++) { |
||
198 | Tree &bigtree = ru->newTree(); |
||
199 | bigtree.setPosition(p); |
||
200 | // add variation: add +/-10% to dbh and *independently* to height. |
||
201 | bigtree.setDbh(dbh * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
||
202 | bigtree.setHeight(tree.height * nrandom(1. - mRecruitmentVariation, 1. + mRecruitmentVariation)); |
||
203 | bigtree.setSpecies( species ); |
||
204 | bigtree.setAge(tree.age.age,tree.height); |
||
205 | bigtree.setRU(ru); |
||
206 | bigtree.setup(); |
||
207 | } |
||
208 | // clear all regeneration from this pixel (including this tree) |
||
209 | ru->clearSaplings(p); |
||
210 | return false; |
||
211 | } |
||
212 | // book keeping (only for survivors) |
||
213 | mLiving++; |
||
214 | mAvgHeight+=tree.height; |
||
215 | mAvgAge+=tree.age.age; |
||
216 | mAvgDeltaHPot+=delta_h_pot; |
||
217 | mAvgHRealized += delta_h_pot * delta_h_factor; |
||
218 | return true; |
||
219 | |||
220 | } |
||
221 | |||
222 | void Sapling::calculateGrowth() |
||
223 | { |
||
224 | Q_ASSERT(mRUS); |
||
225 | if (mLiving==0 && mAdded==0) |
||
226 | return; |
||
227 | |||
228 | clearStatistics(); |
||
229 | ResourceUnit *ru = const_cast<ResourceUnit*> (mRUS->ru() ); |
||
230 | Species *species = const_cast<Species*>(mRUS->species()); |
||
231 | |||
232 | // calculate necessary growth modifier (this is done only once per year) |
||
233 | mRUS->calculate(true); // calculate the 3pg module (this is done only if that did not happen up to now); true: call comes from regeneration |
||
234 | double f_env_yr = mRUS->prod3PG().fEnvYear(); |
||
235 | |||
236 | mLiving=0; |
||
237 | QVector<SaplingTree>::const_iterator it; |
||
238 | for (it = mSaplingTrees.constBegin(); it!=mSaplingTrees.constEnd(); ++it) { |
||
239 | const SaplingTree &tree = *it; |
||
240 | if (tree.height<0) |
||
241 | qDebug() << "Sapling::calculateGrowth(): h<0"; |
||
242 | if (tree.isValid()) { |
||
243 | // growing |
||
244 | if (growSapling(const_cast<SaplingTree&>(tree), f_env_yr, species)) { |
||
245 | // set the sapling height to the maximum value on the current pixel |
||
246 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(tree.pixel); |
||
247 | ru->setMaxSaplingHeightAt(p,tree.height); |
||
248 | } |
||
249 | } |
||
250 | } |
||
251 | if (mLiving) { |
||
252 | mAvgHeight /= double(mLiving); |
||
253 | mAvgAge /= double(mLiving); |
||
254 | mAvgDeltaHPot /= double(mLiving); |
||
255 | mAvgHRealized /= double(mLiving); |
||
256 | } |
||
257 | // calculate carbon balance |
||
258 | mCarbonLiving.clear(); |
||
259 | CNPair dead_wood, dead_fine; // pools for mortality |
||
260 | // average dbh |
||
261 | if (mLiving) { |
||
262 | // calculate the avg dbh number of stems |
||
263 | double avg_dbh = mAvgHeight / species->saplingGrowthParameters().hdSapling * 100.; |
||
264 | double n = mLiving * species->saplingGrowthParameters().representedStemNumber( avg_dbh ); |
||
265 | // woody parts: stem, branchse and coarse roots |
||
266 | double woody_bm = species->biomassWoody(avg_dbh) + species->biomassBranch(avg_dbh) + species->biomassRoot(avg_dbh); |
||
267 | double foliage = species->biomassFoliage(avg_dbh); |
||
268 | double fineroot = foliage*species->finerootFoliageRatio(); |
||
269 | |||
270 | mCarbonLiving.addBiomass( woody_bm*n, species->cnWood() ); |
||
271 | mCarbonLiving.addBiomass( foliage*n, species->cnFoliage() ); |
||
272 | mCarbonLiving.addBiomass( fineroot*n, species->cnFineroot() ); |
||
273 | // turnover |
||
274 | mRUS->ru()->snag()->addTurnoverLitter(species, foliage*species->turnoverLeaf(), fineroot*species->turnoverRoot()); |
||
275 | |||
276 | // calculate the "mortality from competition", i.e. carbon that stems from reduction of stem numbers |
||
277 | // from Reinekes formula. |
||
278 | // |
||
279 | if (avg_dbh>1.) { |
||
280 | double avg_dbh_before = (mAvgHeight - mAvgHRealized) / species->saplingGrowthParameters().hdSapling * 100.; |
||
281 | double n_before = mLiving * species->saplingGrowthParameters().representedStemNumber( qMax(1.,avg_dbh_before) ); |
||
282 | if (n<n_before) { |
||
283 | dead_wood.addBiomass( woody_bm * (n_before-n), species->cnWood() ); |
||
284 | dead_fine.addBiomass( foliage * (n_before-n), species->cnFoliage() ); |
||
285 | dead_fine.addBiomass( fineroot * (n_before-n), species->cnFineroot() ); |
||
286 | } |
||
287 | } |
||
288 | |||
289 | } |
||
290 | if (mDied) { |
||
291 | double avg_dbh_dead = mSumDbhDied / double(mDied); |
||
292 | double n = mDied * species->saplingGrowthParameters().representedStemNumber( avg_dbh_dead ); |
||
293 | // woody parts: stem, branchse and coarse roots |
||
294 | |||
295 | dead_wood.addBiomass( ( species->biomassWoody(avg_dbh_dead) + species->biomassBranch(avg_dbh_dead) + species->biomassRoot(avg_dbh_dead)) * n, species->cnWood() ); |
||
296 | double foliage = species->biomassFoliage(avg_dbh_dead)*n; |
||
297 | |||
298 | dead_fine.addBiomass( foliage, species->cnFoliage() ); |
||
299 | dead_fine.addBiomass( foliage*species->finerootFoliageRatio(), species->cnFineroot() ); |
||
300 | } |
||
301 | if (!dead_wood.isEmpty() || !dead_fine.isEmpty()) |
||
302 | mRUS->ru()->snag()->addToSoil(species, dead_wood, dead_fine); |
||
303 | |||
304 | if (mSaplingTrees.count() > mLiving*1.3) |
||
305 | cleanupStorage(); |
||
306 | |||
307 | mRUS->statistics().add(this); |
||
308 | //qDebug() << ru->index() << species->id()<< ": (living/avg.height):" << mLiving << mAvgHeight; |
||
309 | } |
||
310 | |||
311 | /// fill a grid with the maximum height of saplings per pixel (2x2m). |
||
312 | /// this function is used for visualization only |
||
313 | void Sapling::fillMaxHeightGrid(Grid<float> &grid) const |
||
314 | { |
||
315 | QVector<SaplingTree>::const_iterator it; |
||
316 | for (it = mSaplingTrees.begin(); it!=mSaplingTrees.end(); ++it) { |
||
317 | if (it->isValid()) { |
||
318 | QPoint p=GlobalSettings::instance()->model()->grid()->indexOf(it->pixel); |
||
319 | if (grid.valueAtIndex(p)<it->height) |
||
320 | grid.valueAtIndex(p) = it->height; |
||
321 | } |
||
322 | } |
||
323 | |||
324 | } |