Rev 593 | Rev 600 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 593 | Rev 595 | ||
---|---|---|---|
Line 1... | Line 1... | ||
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/sapling.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/sapling.cpp': |
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()->addRegeneration(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 | }
|
- | |
- | 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 | }
|