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