Rev 1104 | Rev 1118 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed
Rev 1104 | Rev 1114 | ||
---|---|---|---|
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.cpp': |
1 | Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/tree.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 | #include "global.h"
|
20 | #include "global.h"
|
21 | #include "tree.h"
|
21 | #include "tree.h"
|
22 | 22 | ||
23 | #include "grid.h"
|
23 | #include "grid.h"
|
24 | 24 | ||
25 | #include "stamp.h"
|
25 | #include "stamp.h"
|
26 | #include "species.h"
|
26 | #include "species.h"
|
27 | #include "resourceunit.h"
|
27 | #include "resourceunit.h"
|
28 | #include "model.h"
|
28 | #include "model.h"
|
29 | #include "snag.h"
|
29 | #include "snag.h"
|
30 | 30 | ||
31 | #include "forestmanagementengine.h"
|
31 | #include "forestmanagementengine.h"
|
32 | #include "modules.h"
|
32 | #include "modules.h"
|
33 | 33 | ||
34 | #include "treeout.h"
|
34 | #include "treeout.h"
|
- | 35 | #include "landscapeout.h"
|
|
35 | 36 | ||
36 | // static varaibles
|
37 | // static varaibles
|
37 | FloatGrid *Tree::mGrid = 0; |
38 | FloatGrid *Tree::mGrid = 0; |
38 | HeightGrid *Tree::mHeightGrid = 0; |
39 | HeightGrid *Tree::mHeightGrid = 0; |
39 | TreeRemovedOut *Tree::mRemovalOutput = 0; |
40 | TreeRemovedOut *Tree::mRemovalOutput = 0; |
- | 41 | LandscapeRemovedOut *Tree::mLSRemovalOutput = 0; |
|
40 | int Tree::m_statPrint=0; |
42 | int Tree::m_statPrint=0; |
41 | int Tree::m_statAboveZ=0; |
43 | int Tree::m_statAboveZ=0; |
42 | int Tree::m_statCreated=0; |
44 | int Tree::m_statCreated=0; |
43 | int Tree::m_nextId=0; |
45 | int Tree::m_nextId=0; |
44 | 46 | ||
45 | #ifdef ALT_TREE_MORTALITY
|
47 | #ifdef ALT_TREE_MORTALITY
|
46 | static double _stress_threshold = 0.05; |
48 | static double _stress_threshold = 0.05; |
47 | static int _stress_years = 5; |
49 | static int _stress_years = 5; |
48 | static double _stress_death_prob = 0.33; |
50 | static double _stress_death_prob = 0.33; |
49 | #endif
|
51 | #endif
|
50 | 52 | ||
51 | /** @class Tree
|
53 | /** @class Tree
|
52 | @ingroup core
|
54 | @ingroup core
|
53 | A tree is the basic simulation entity of iLand and represents a single tree.
|
55 | A tree is the basic simulation entity of iLand and represents a single tree.
|
54 | Trees in iLand are designed to be lightweight, thus the list of stored properties is limited. Basic properties
|
56 | Trees in iLand are designed to be lightweight, thus the list of stored properties is limited. Basic properties
|
55 | are dimensions (dbh, height), biomass pools (stem, leaves, roots), the reserve NPP pool. Additionally, the location and species are stored.
|
57 | are dimensions (dbh, height), biomass pools (stem, leaves, roots), the reserve NPP pool. Additionally, the location and species are stored.
|
56 | A Tree has a height of at least 4m; trees below this threshold are covered by the regeneration layer (see Sapling).
|
58 | A Tree has a height of at least 4m; trees below this threshold are covered by the regeneration layer (see Sapling).
|
57 | Trees are stored in lists managed at the resource unit level.
|
59 | Trees are stored in lists managed at the resource unit level.
|
58 | 60 | ||
59 | */
|
61 | */
|
60 | 62 | ||
61 | /** get distance and direction between two points.
|
63 | /** get distance and direction between two points.
|
62 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */
|
64 | returns the distance (m), and the angle between PStart and PEnd (radians) in referenced param rAngle. */
|
63 | double dist_and_direction(const QPointF &PStart, const QPointF &PEnd, double &rAngle) |
65 | double dist_and_direction(const QPointF &PStart, const QPointF &PEnd, double &rAngle) |
64 | {
|
66 | {
|
65 | double dx = PEnd.x() - PStart.x(); |
67 | double dx = PEnd.x() - PStart.x(); |
66 | double dy = PEnd.y() - PStart.y(); |
68 | double dy = PEnd.y() - PStart.y(); |
67 | double d = sqrt(dx*dx + dy*dy); |
69 | double d = sqrt(dx*dx + dy*dy); |
68 | // direction:
|
70 | // direction:
|
69 | rAngle = atan2(dx, dy); |
71 | rAngle = atan2(dx, dy); |
70 | return d; |
72 | return d; |
71 | }
|
73 | }
|
72 | 74 | ||
73 | 75 | ||
74 | // lifecycle
|
76 | // lifecycle
|
75 | Tree::Tree() |
77 | Tree::Tree() |
76 | {
|
78 | {
|
77 | mDbh = mHeight = 0; |
79 | mDbh = mHeight = 0; |
78 | mRU = 0; mSpecies = 0; |
80 | mRU = 0; mSpecies = 0; |
79 | mFlags = mAge = 0; |
81 | mFlags = mAge = 0; |
80 | mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.; |
82 | mOpacity=mFoliageMass=mWoodyMass=mCoarseRootMass=mFineRootMass=mLeafArea=0.; |
81 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
83 | mDbhDelta=mNPPReserve=mLRI=mStressIndex=0.; |
82 | mLightResponse = 0.; |
84 | mLightResponse = 0.; |
83 | mStamp=0; |
85 | mStamp=0; |
84 | mId = m_nextId++; |
86 | mId = m_nextId++; |
85 | m_statCreated++; |
87 | m_statCreated++; |
86 | }
|
88 | }
|
87 | 89 | ||
88 | float Tree::crownRadius() const |
90 | float Tree::crownRadius() const |
89 | {
|
91 | {
|
90 | Q_ASSERT(mStamp!=0); |
92 | Q_ASSERT(mStamp!=0); |
91 | return mStamp->crownRadius(); |
93 | return mStamp->crownRadius(); |
92 | }
|
94 | }
|
93 | 95 | ||
94 | float Tree::biomassBranch() const |
96 | float Tree::biomassBranch() const |
95 | {
|
97 | {
|
96 | return static_cast<float>( mSpecies->biomassBranch(mDbh) ); |
98 | return static_cast<float>( mSpecies->biomassBranch(mDbh) ); |
97 | }
|
99 | }
|
98 | 100 | ||
99 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
101 | void Tree::setGrid(FloatGrid* gridToStamp, Grid<HeightGridValue> *dominanceGrid) |
100 | {
|
102 | {
|
101 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
103 | mGrid = gridToStamp; mHeightGrid = dominanceGrid; |
102 | }
|
104 | }
|
103 | 105 | ||
104 | // calculate the thickness of the bark of the tree
|
106 | // calculate the thickness of the bark of the tree
|
105 | double Tree::barkThickness() const |
107 | double Tree::barkThickness() const |
106 | {
|
108 | {
|
107 | return mSpecies->barkThickness(mDbh); |
109 | return mSpecies->barkThickness(mDbh); |
108 | }
|
110 | }
|
109 | 111 | ||
110 | /// dumps some core variables of a tree to a string.
|
112 | /// dumps some core variables of a tree to a string.
|
111 | QString Tree::dump() |
113 | QString Tree::dump() |
112 | {
|
114 | {
|
113 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
115 | QString result = QString("id %1 species %2 dbh %3 h %4 x/y %5/%6 ru# %7 LRI %8") |
114 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
116 | .arg(mId).arg(species()->id()).arg(mDbh).arg(mHeight) |
115 | .arg(position().x()).arg(position().y()) |
117 | .arg(position().x()).arg(position().y()) |
116 | .arg(mRU->index()).arg(mLRI); |
118 | .arg(mRU->index()).arg(mLRI); |
117 | return result; |
119 | return result; |
118 | }
|
120 | }
|
119 | 121 | ||
120 | void Tree::dumpList(DebugList &rTargetList) |
122 | void Tree::dumpList(DebugList &rTargetList) |
121 | {
|
123 | {
|
122 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
124 | rTargetList << mId << species()->id() << mDbh << mHeight << position().x() << position().y() << mRU->index() << mLRI |
123 | << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea; |
125 | << mWoodyMass << mCoarseRootMass << mFoliageMass << mLeafArea; |
124 | }
|
126 | }
|
125 | 127 | ||
126 | void Tree::setup() |
128 | void Tree::setup() |
127 | {
|
129 | {
|
128 | if (mDbh<=0 || mHeight<=0) { |
130 | if (mDbh<=0 || mHeight<=0) { |
129 | throw IException(QString("Error: trying to set up a tree with invalid dimensions: dbh: %1 height: %2 id: %3 RU-index: %4").arg(mDbh).arg(mHeight).arg(mId).arg(mRU->index())); |
131 | throw IException(QString("Error: trying to set up a tree with invalid dimensions: dbh: %1 height: %2 id: %3 RU-index: %4").arg(mDbh).arg(mHeight).arg(mId).arg(mRU->index())); |
130 | }
|
132 | }
|
131 | // check stamp
|
133 | // check stamp
|
132 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
134 | Q_ASSERT_X(species()!=0, "Tree::setup()", "species is NULL"); |
133 | mStamp = species()->stamp(mDbh, mHeight); |
135 | mStamp = species()->stamp(mDbh, mHeight); |
134 | if (!mStamp) { |
136 | if (!mStamp) { |
135 | throw IException("Tree::setup() with invalid stamp!"); |
137 | throw IException("Tree::setup() with invalid stamp!"); |
136 | }
|
138 | }
|
137 | 139 | ||
138 | mFoliageMass = static_cast<float>( species()->biomassFoliage(mDbh) ); |
140 | mFoliageMass = static_cast<float>( species()->biomassFoliage(mDbh) ); |
139 | mCoarseRootMass = static_cast<float>( species()->biomassRoot(mDbh) ); // coarse root (allometry) |
141 | mCoarseRootMass = static_cast<float>( species()->biomassRoot(mDbh) ); // coarse root (allometry) |
140 | mFineRootMass = static_cast<float>( mFoliageMass * species()->finerootFoliageRatio() ); // fine root (size defined by finerootFoliageRatio) |
142 | mFineRootMass = static_cast<float>( mFoliageMass * species()->finerootFoliageRatio() ); // fine root (size defined by finerootFoliageRatio) |
141 | mWoodyMass = static_cast<float>( species()->biomassWoody(mDbh) ); |
143 | mWoodyMass = static_cast<float>( species()->biomassWoody(mDbh) ); |
142 | 144 | ||
143 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
|
145 | // LeafArea[m2] = LeafMass[kg] * specificLeafArea[m2/kg]
|
144 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); |
146 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); |
145 | mOpacity = static_cast<float>( 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea()) ); |
147 | mOpacity = static_cast<float>( 1. - exp(- Model::settings().lightExtinctionCoefficientOpacity * mLeafArea / mStamp->crownArea()) ); |
146 | mNPPReserve = static_cast<float>( (1+species()->finerootFoliageRatio())*mFoliageMass ); // initial value |
148 | mNPPReserve = static_cast<float>( (1+species()->finerootFoliageRatio())*mFoliageMass ); // initial value |
147 | mDbhDelta = 0.1f; // initial value: used in growth() to estimate diameter increment |
149 | mDbhDelta = 0.1f; // initial value: used in growth() to estimate diameter increment |
148 | 150 | ||
149 | }
|
151 | }
|
150 | 152 | ||
151 | void Tree::setAge(const int age, const float treeheight) |
153 | void Tree::setAge(const int age, const float treeheight) |
152 | {
|
154 | {
|
153 | mAge = age; |
155 | mAge = age; |
154 | if (age==0) { |
156 | if (age==0) { |
155 | // estimate age using the tree height
|
157 | // estimate age using the tree height
|
156 | mAge = mSpecies->estimateAge(treeheight); |
158 | mAge = mSpecies->estimateAge(treeheight); |
157 | }
|
159 | }
|
158 | }
|
160 | }
|
159 | 161 | ||
160 | //////////////////////////////////////////////////
|
162 | //////////////////////////////////////////////////
|
161 | //// Light functions (Pattern-stuff)
|
163 | //// Light functions (Pattern-stuff)
|
162 | //////////////////////////////////////////////////
|
164 | //////////////////////////////////////////////////
|
163 | 165 | ||
164 | #define NOFULLDBG
|
166 | #define NOFULLDBG
|
165 | //#define NOFULLOPT
|
167 | //#define NOFULLOPT
|
166 | 168 | ||
167 | 169 | ||
168 | void Tree::applyLIP() |
170 | void Tree::applyLIP() |
169 | {
|
171 | {
|
170 | if (!mStamp) |
172 | if (!mStamp) |
171 | return; |
173 | return; |
172 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
174 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
173 | QPoint pos = mPositionIndex; |
175 | QPoint pos = mPositionIndex; |
174 | int offset = mStamp->offset(); |
176 | int offset = mStamp->offset(); |
175 | pos-=QPoint(offset, offset); |
177 | pos-=QPoint(offset, offset); |
176 | 178 | ||
177 | float local_dom; // height of Z* on the current position |
179 | float local_dom; // height of Z* on the current position |
178 | int x,y; |
180 | int x,y; |
179 | float value, z, z_zstar; |
181 | float value, z, z_zstar; |
180 | int gr_stamp = mStamp->size(); |
182 | int gr_stamp = mStamp->size(); |
181 | 183 | ||
182 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
184 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
183 | // this should not happen because of the buffer
|
185 | // this should not happen because of the buffer
|
184 | return; |
186 | return; |
185 | }
|
187 | }
|
186 | int grid_y = pos.y(); |
188 | int grid_y = pos.y(); |
187 | for (y=0;y<gr_stamp; ++y) { |
189 | for (y=0;y<gr_stamp; ++y) { |
188 | 190 | ||
189 | float *grid_value_ptr = mGrid->ptr(pos.x(), grid_y); |
191 | float *grid_value_ptr = mGrid->ptr(pos.x(), grid_y); |
190 | int grid_x = pos.x(); |
192 | int grid_x = pos.x(); |
191 | for (x=0;x<gr_stamp;++x, ++grid_x, ++grid_value_ptr) { |
193 | for (x=0;x<gr_stamp;++x, ++grid_x, ++grid_value_ptr) { |
192 | // suppose there is no stamping outside
|
194 | // suppose there is no stamping outside
|
193 | value = (*mStamp)(x,y); // stampvalue |
195 | value = (*mStamp)(x,y); // stampvalue |
194 | //if (value>0.f) {
|
196 | //if (value>0.f) {
|
195 | local_dom = (*mHeightGrid)(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height; |
197 | local_dom = (*mHeightGrid)(grid_x/cPxPerHeight, grid_y/cPxPerHeight).height; |
196 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
198 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
197 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
199 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
198 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
200 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
199 | value = std::max(value, 0.02f); // limit value |
201 | value = std::max(value, 0.02f); // limit value |
200 | 202 | ||
201 | *grid_value_ptr *= value; |
203 | *grid_value_ptr *= value; |
202 | //}
|
204 | //}
|
203 | 205 | ||
204 | }
|
206 | }
|
205 | grid_y++; |
207 | grid_y++; |
206 | }
|
208 | }
|
207 | 209 | ||
208 | m_statPrint++; // count # of stamp applications... |
210 | m_statPrint++; // count # of stamp applications... |
209 | }
|
211 | }
|
210 | 212 | ||
211 | /// helper function for gluing the edges together
|
213 | /// helper function for gluing the edges together
|
212 | /// index: index at grid
|
214 | /// index: index at grid
|
213 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)
|
215 | /// count: number of pixels that are the simulation area (e.g. 100m and 2m pixel -> 50)
|
214 | /// buffer: size of buffer around simulation area (in pixels)
|
216 | /// buffer: size of buffer around simulation area (in pixels)
|
215 | inline int torusIndex(int index, int count, int buffer, int ru_index) |
217 | inline int torusIndex(int index, int count, int buffer, int ru_index) |
216 | {
|
218 | {
|
217 | return buffer + ru_index + (index-buffer+count)%count; |
219 | return buffer + ru_index + (index-buffer+count)%count; |
218 | }
|
220 | }
|
219 | 221 | ||
220 | 222 | ||
221 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.
|
223 | /** Apply LIPs. This "Torus" functions wraps the influence at the edges of a 1ha simulation area.
|
222 | */
|
224 | */
|
223 | void Tree::applyLIP_torus() |
225 | void Tree::applyLIP_torus() |
224 | {
|
226 | {
|
225 | if (!mStamp) |
227 | if (!mStamp) |
226 | return; |
228 | return; |
227 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
229 | Q_ASSERT(mGrid!=0 && mStamp!=0 && mRU!=0); |
228 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
230 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
229 | QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
231 | QPoint pos = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
230 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
232 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
231 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index |
233 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos.x(), mPositionIndex.y() - pos.y()); // offset of the corner of the resource index |
232 | 234 | ||
233 | int offset = mStamp->offset(); |
235 | int offset = mStamp->offset(); |
234 | pos-=QPoint(offset, offset); |
236 | pos-=QPoint(offset, offset); |
235 | 237 | ||
236 | float local_dom; // height of Z* on the current position |
238 | float local_dom; // height of Z* on the current position |
237 | int x,y; |
239 | int x,y; |
238 | float value; |
240 | float value; |
239 | int gr_stamp = mStamp->size(); |
241 | int gr_stamp = mStamp->size(); |
240 | int grid_x, grid_y; |
242 | int grid_x, grid_y; |
241 | float *grid_value; |
243 | float *grid_value; |
242 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
244 | if (!mGrid->isIndexValid(pos) || !mGrid->isIndexValid(pos+QPoint(gr_stamp, gr_stamp))) { |
243 | // todo: in this case we should use another algorithm!!! necessary????
|
245 | // todo: in this case we should use another algorithm!!! necessary????
|
244 | return; |
246 | return; |
245 | }
|
247 | }
|
246 | float z, z_zstar; |
248 | float z, z_zstar; |
247 | int xt, yt; // wraparound coordinates |
249 | int xt, yt; // wraparound coordinates |
248 | for (y=0;y<gr_stamp; ++y) { |
250 | for (y=0;y<gr_stamp; ++y) { |
249 | grid_y = pos.y() + y; |
251 | grid_y = pos.y() + y; |
250 | yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m |
252 | yt = torusIndex(grid_y, cPxPerRU,bufferOffset, ru_offset.y()); // 50 cells per 100m |
251 | for (x=0;x<gr_stamp;++x) { |
253 | for (x=0;x<gr_stamp;++x) { |
252 | // suppose there is no stamping outside
|
254 | // suppose there is no stamping outside
|
253 | grid_x = pos.x() + x; |
255 | grid_x = pos.x() + x; |
254 | xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x()); |
256 | xt = torusIndex(grid_x,cPxPerRU,bufferOffset, ru_offset.x()); |
255 | 257 | ||
256 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height; |
258 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight,yt/cPxPerHeight).height; |
257 | 259 | ||
258 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
260 | z = std::max(mHeight - (*mStamp).distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
259 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
261 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
260 | value = (*mStamp)(x,y); // stampvalue |
262 | value = (*mStamp)(x,y); // stampvalue |
261 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
263 | value = 1.f - value*mOpacity * z_zstar; // calculated value |
262 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
264 | // old: value = 1. - value*mOpacity / local_dom; // calculated value
|
263 | value = qMax(value, 0.02f); // limit value |
265 | value = qMax(value, 0.02f); // limit value |
264 | 266 | ||
265 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
267 | grid_value = mGrid->ptr(xt, yt); // use wraparound coordinates |
266 | *grid_value *= value; |
268 | *grid_value *= value; |
267 | }
|
269 | }
|
268 | }
|
270 | }
|
269 | 271 | ||
270 | m_statPrint++; // count # of stamp applications... |
272 | m_statPrint++; // count # of stamp applications... |
271 | }
|
273 | }
|
272 | 274 | ||
273 | /** heightGrid()
|
275 | /** heightGrid()
|
274 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.
|
276 | This function calculates the "dominant height field". This grid is coarser as the fine-scaled light-grid.
|
275 | */
|
277 | */
|
276 | void Tree::heightGrid() |
278 | void Tree::heightGrid() |
277 | {
|
279 | {
|
278 | 280 | ||
279 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
281 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
280 | 282 | ||
281 | // count trees that are on height-grid cells (used for stockable area)
|
283 | // count trees that are on height-grid cells (used for stockable area)
|
282 | mHeightGrid->valueAtIndex(p).increaseCount(); |
284 | mHeightGrid->valueAtIndex(p).increaseCount(); |
283 | if (mHeight > mHeightGrid->valueAtIndex(p).height) |
285 | if (mHeight > mHeightGrid->valueAtIndex(p).height) |
284 | mHeightGrid->valueAtIndex(p).height=mHeight; |
286 | mHeightGrid->valueAtIndex(p).height=mHeight; |
285 | 287 | ||
286 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
288 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
287 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
289 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
288 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
290 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
289 | if (index_eastwest - r < 0) { // east |
291 | if (index_eastwest - r < 0) { // east |
290 | mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight); |
292 | mHeightGrid->valueAtIndex(p.x()-1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()-1, p.y()).height,mHeight); |
291 | }
|
293 | }
|
292 | if (index_eastwest + r >= cPxPerHeight) { // west |
294 | if (index_eastwest + r >= cPxPerHeight) { // west |
293 | mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight); |
295 | mHeightGrid->valueAtIndex(p.x()+1, p.y()).height=qMax(mHeightGrid->valueAtIndex(p.x()+1, p.y()).height,mHeight); |
294 | }
|
296 | }
|
295 | if (index_northsouth - r < 0) { // south |
297 | if (index_northsouth - r < 0) { // south |
296 | mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight); |
298 | mHeightGrid->valueAtIndex(p.x(), p.y()-1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()-1).height,mHeight); |
297 | }
|
299 | }
|
298 | if (index_northsouth + r >= cPxPerHeight) { // north |
300 | if (index_northsouth + r >= cPxPerHeight) { // north |
299 | mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight); |
301 | mHeightGrid->valueAtIndex(p.x(), p.y()+1).height=qMax(mHeightGrid->valueAtIndex(p.x(), p.y()+1).height,mHeight); |
300 | }
|
302 | }
|
301 | 303 | ||
302 | 304 | ||
303 | // without spread of the height grid
|
305 | // without spread of the height grid
|
304 | 306 | ||
305 | // // height of Z*
|
307 | // // height of Z*
|
306 | // const float cellsize = mHeightGrid->cellsize();
|
308 | // const float cellsize = mHeightGrid->cellsize();
|
307 | //
|
309 | //
|
308 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
310 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
309 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
311 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
310 | // int dist[9];
|
312 | // int dist[9];
|
311 | // dist[3] = index_northsouth * 2 + 1; // south
|
313 | // dist[3] = index_northsouth * 2 + 1; // south
|
312 | // dist[1] = index_eastwest * 2 + 1; // west
|
314 | // dist[1] = index_eastwest * 2 + 1; // west
|
313 | // dist[5] = 10 - dist[3]; // north
|
315 | // dist[5] = 10 - dist[3]; // north
|
314 | // dist[7] = 10 - dist[1]; // east
|
316 | // dist[7] = 10 - dist[1]; // east
|
315 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
317 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
316 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
318 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
317 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
319 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
318 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
320 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
319 | // dist[4] = 0; // center cell
|
321 | // dist[4] = 0; // center cell
|
320 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
322 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
321 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
323 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
322 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
324 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
323 | // */
|
325 | // */
|
324 | //
|
326 | //
|
325 | //
|
327 | //
|
326 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
328 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
327 | // int ix, iy;
|
329 | // int ix, iy;
|
328 | // int ring;
|
330 | // int ring;
|
329 | // float hdom;
|
331 | // float hdom;
|
330 | //
|
332 | //
|
331 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
333 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
332 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
334 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
333 | // ring = qMax(abs(ix), abs(iy));
|
335 | // ring = qMax(abs(ix), abs(iy));
|
334 | // QPoint pos(ix+p.x(), iy+p.y());
|
336 | // QPoint pos(ix+p.x(), iy+p.y());
|
335 | // if (mHeightGrid->isIndexValid(pos)) {
|
337 | // if (mHeightGrid->isIndexValid(pos)) {
|
336 | // float &rHGrid = mHeightGrid->valueAtIndex(pos).height;
|
338 | // float &rHGrid = mHeightGrid->valueAtIndex(pos).height;
|
337 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
339 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
338 | // continue;
|
340 | // continue;
|
339 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
341 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
340 | // hdom = mHeight - dist[direction];
|
342 | // hdom = mHeight - dist[direction];
|
341 | // if (ring>1)
|
343 | // if (ring>1)
|
342 | // hdom -= (ring-1)*10;
|
344 | // hdom -= (ring-1)*10;
|
343 | //
|
345 | //
|
344 | // rHGrid = qMax(rHGrid, hdom); // write value
|
346 | // rHGrid = qMax(rHGrid, hdom); // write value
|
345 | // } // is valid
|
347 | // } // is valid
|
346 | // } // for (y)
|
348 | // } // for (y)
|
347 | }
|
349 | }
|
348 | 350 | ||
349 | void Tree::heightGrid_torus() |
351 | void Tree::heightGrid_torus() |
350 | {
|
352 | {
|
351 | // height of Z*
|
353 | // height of Z*
|
352 | 354 | ||
353 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
355 | QPoint p = QPoint(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight); // pos of tree on height grid |
354 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels) |
356 | int bufferOffset = mHeightGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer (i.e.: size of buffer in height-pixels) |
355 | p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m |
357 | p.setX((p.x()-bufferOffset)%10 + bufferOffset); // 10: 10 x 10m pixeln in 100m |
356 | p.setY((p.y()-bufferOffset)%10 + bufferOffset); |
358 | p.setY((p.y()-bufferOffset)%10 + bufferOffset); |
357 | 359 | ||
358 | 360 | ||
359 | // torus coordinates: ru_offset = coords of lower left corner of 1ha patch
|
361 | // torus coordinates: ru_offset = coords of lower left corner of 1ha patch
|
360 | QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y()); |
362 | QPoint ru_offset =QPoint(mPositionIndex.x()/cPxPerHeight - p.x(), mPositionIndex.y()/cPxPerHeight - p.y()); |
361 | 363 | ||
362 | // count trees that are on height-grid cells (used for stockable area)
|
364 | // count trees that are on height-grid cells (used for stockable area)
|
363 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
365 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
364 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
366 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
365 | v.increaseCount(); |
367 | v.increaseCount(); |
366 | v.height = qMax(v.height, mHeight); |
368 | v.height = qMax(v.height, mHeight); |
367 | 369 | ||
368 | 370 | ||
369 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
371 | int r = mStamp->reader()->offset(); // distance between edge and the center pixel. e.g.: if r = 2 -> stamp=5x5 |
370 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
372 | int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge |
371 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
373 | int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge |
372 | if (index_eastwest - r < 0) { // east |
374 | if (index_eastwest - r < 0) { // east |
373 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()), |
375 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()-1,10,bufferOffset,ru_offset.x()), |
374 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
376 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
375 | v.height = qMax(v.height, mHeight); |
377 | v.height = qMax(v.height, mHeight); |
376 | }
|
378 | }
|
377 | if (index_eastwest + r >= cPxPerHeight) { // west |
379 | if (index_eastwest + r >= cPxPerHeight) { // west |
378 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()), |
380 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x()+1,10,bufferOffset,ru_offset.x()), |
379 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
381 | torusIndex(p.y(),10,bufferOffset,ru_offset.y())); |
380 | v.height = qMax(v.height, mHeight); |
382 | v.height = qMax(v.height, mHeight); |
381 | }
|
383 | }
|
382 | if (index_northsouth - r < 0) { // south |
384 | if (index_northsouth - r < 0) { // south |
383 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
385 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
384 | torusIndex(p.y()-1,10,bufferOffset,ru_offset.y())); |
386 | torusIndex(p.y()-1,10,bufferOffset,ru_offset.y())); |
385 | v.height = qMax(v.height, mHeight); |
387 | v.height = qMax(v.height, mHeight); |
386 | }
|
388 | }
|
387 | if (index_northsouth + r >= cPxPerHeight) { // north |
389 | if (index_northsouth + r >= cPxPerHeight) { // north |
388 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
390 | HeightGridValue &v = mHeightGrid->valueAtIndex(torusIndex(p.x(),10,bufferOffset,ru_offset.x()), |
389 | torusIndex(p.y()+1,10,bufferOffset,ru_offset.y())); |
391 | torusIndex(p.y()+1,10,bufferOffset,ru_offset.y())); |
390 | v.height = qMax(v.height, mHeight); |
392 | v.height = qMax(v.height, mHeight); |
391 | }
|
393 | }
|
392 | 394 | ||
393 | 395 | ||
394 | 396 | ||
395 | 397 | ||
396 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
398 | // int index_eastwest = mPositionIndex.x() % cPxPerHeight; // 4: very west, 0 east edge
|
397 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
399 | // int index_northsouth = mPositionIndex.y() % cPxPerHeight; // 4: northern edge, 0: southern edge
|
398 | // int dist[9];
|
400 | // int dist[9];
|
399 | // dist[3] = index_northsouth * 2 + 1; // south
|
401 | // dist[3] = index_northsouth * 2 + 1; // south
|
400 | // dist[1] = index_eastwest * 2 + 1; // west
|
402 | // dist[1] = index_eastwest * 2 + 1; // west
|
401 | // dist[5] = 10 - dist[3]; // north
|
403 | // dist[5] = 10 - dist[3]; // north
|
402 | // dist[7] = 10 - dist[1]; // east
|
404 | // dist[7] = 10 - dist[1]; // east
|
403 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
405 | // dist[8] = qMax(dist[5], dist[7]); // north-east
|
404 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
406 | // dist[6] = qMax(dist[3], dist[7]); // south-east
|
405 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
407 | // dist[0] = qMax(dist[3], dist[1]); // south-west
|
406 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
408 | // dist[2] = qMax(dist[5], dist[1]); // north-west
|
407 | // dist[4] = 0; // center cell
|
409 | // dist[4] = 0; // center cell
|
408 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
410 | // /* the scheme of indices is as follows: if sign(ix)= -1, if ix<0, 0 for ix=0, 1 for ix>0 (detto iy), then:
|
409 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
411 | // index = 4 + 3*sign(ix) + sign(iy) transforms combinations of directions to unique ids (0..8), which are used above.
|
410 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
412 | // e.g.: sign(ix) = -1, sign(iy) = 1 (=north-west) -> index = 4 + -3 + 1 = 2
|
411 | // */
|
413 | // */
|
412 | //
|
414 | //
|
413 | //
|
415 | //
|
414 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
416 | // int ringcount = int(floor(mHeight / cellsize)) + 1;
|
415 | // int ix, iy;
|
417 | // int ix, iy;
|
416 | // int ring;
|
418 | // int ring;
|
417 | // float hdom;
|
419 | // float hdom;
|
418 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
420 | // for (ix=-ringcount;ix<=ringcount;ix++)
|
419 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
421 | // for (iy=-ringcount; iy<=+ringcount; iy++) {
|
420 | // ring = qMax(abs(ix), abs(iy));
|
422 | // ring = qMax(abs(ix), abs(iy));
|
421 | // QPoint pos(ix+p.x(), iy+p.y());
|
423 | // QPoint pos(ix+p.x(), iy+p.y());
|
422 | // QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()),
|
424 | // QPoint p_torus(torusIndex(pos.x(),10,bufferOffset,ru_offset.x()),
|
423 | // torusIndex(pos.y(),10,bufferOffset,ru_offset.y()));
|
425 | // torusIndex(pos.y(),10,bufferOffset,ru_offset.y()));
|
424 | // if (mHeightGrid->isIndexValid(p_torus)) {
|
426 | // if (mHeightGrid->isIndexValid(p_torus)) {
|
425 | // float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height;
|
427 | // float &rHGrid = mHeightGrid->valueAtIndex(p_torus.x(),p_torus.y()).height;
|
426 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
428 | // if (rHGrid > mHeight) // skip calculation if grid is higher than tree
|
427 | // continue;
|
429 | // continue;
|
428 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
430 | // int direction = 4 + (ix?(ix<0?-3:3):0) + (iy?(iy<0?-1:1):0); // 4 + 3*sgn(x) + sgn(y)
|
429 | // hdom = mHeight - dist[direction];
|
431 | // hdom = mHeight - dist[direction];
|
430 | // if (ring>1)
|
432 | // if (ring>1)
|
431 | // hdom -= (ring-1)*10;
|
433 | // hdom -= (ring-1)*10;
|
432 | //
|
434 | //
|
433 | // rHGrid = qMax(rHGrid, hdom); // write value
|
435 | // rHGrid = qMax(rHGrid, hdom); // write value
|
434 | // } // is valid
|
436 | // } // is valid
|
435 | // } // for (y)
|
437 | // } // for (y)
|
436 | }
|
438 | }
|
437 | 439 | ||
438 | 440 | ||
439 | /** reads the light influence field value for a tree.
|
441 | /** reads the light influence field value for a tree.
|
440 | The LIF field is scanned within the crown area of the focal tree, and the influence of
|
442 | The LIF field is scanned within the crown area of the focal tree, and the influence of
|
441 | the focal tree is "subtracted" from the LIF values.
|
443 | the focal tree is "subtracted" from the LIF values.
|
442 | Finally, the "LRI correction" is applied.
|
444 | Finally, the "LRI correction" is applied.
|
443 | see http://iland.boku.ac.at/competition+for+light for details.
|
445 | see http://iland.boku.ac.at/competition+for+light for details.
|
444 | */
|
446 | */
|
445 | void Tree::readLIF() |
447 | void Tree::readLIF() |
446 | {
|
448 | {
|
447 | if (!mStamp) |
449 | if (!mStamp) |
448 | return; |
450 | return; |
449 | const Stamp *reader = mStamp->reader(); |
451 | const Stamp *reader = mStamp->reader(); |
450 | if (!reader) |
452 | if (!reader) |
451 | return; |
453 | return; |
452 | QPoint pos_reader = mPositionIndex; |
454 | QPoint pos_reader = mPositionIndex; |
453 | const float outside_area_factor = 0.1f; // |
455 | const float outside_area_factor = 0.1f; // |
454 | 456 | ||
455 | int offset_reader = reader->offset(); |
457 | int offset_reader = reader->offset(); |
456 | int offset_writer = mStamp->offset(); |
458 | int offset_writer = mStamp->offset(); |
457 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
459 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
458 | 460 | ||
459 | pos_reader-=QPoint(offset_reader, offset_reader); |
461 | pos_reader-=QPoint(offset_reader, offset_reader); |
460 | 462 | ||
461 | float local_dom; |
463 | float local_dom; |
462 | 464 | ||
463 | int x,y; |
465 | int x,y; |
464 | double sum=0.; |
466 | double sum=0.; |
465 | double value, own_value; |
467 | double value, own_value; |
466 | float *grid_value; |
468 | float *grid_value; |
467 | float z, z_zstar; |
469 | float z, z_zstar; |
468 | int reader_size = reader->size(); |
470 | int reader_size = reader->size(); |
469 | int rx = pos_reader.x(); |
471 | int rx = pos_reader.x(); |
470 | int ry = pos_reader.y(); |
472 | int ry = pos_reader.y(); |
471 | for (y=0;y<reader_size; ++y, ++ry) { |
473 | for (y=0;y<reader_size; ++y, ++ry) { |
472 | grid_value = mGrid->ptr(rx, ry); |
474 | grid_value = mGrid->ptr(rx, ry); |
473 | for (x=0;x<reader_size;++x) { |
475 | for (x=0;x<reader_size;++x) { |
474 | 476 | ||
475 | const HeightGridValue &hgv = mHeightGrid->constValueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight); // the height grid value, ry: gets ++ed in outer loop, rx not |
477 | const HeightGridValue &hgv = mHeightGrid->constValueAtIndex((rx+x)/cPxPerHeight, ry/cPxPerHeight); // the height grid value, ry: gets ++ed in outer loop, rx not |
476 | local_dom = hgv.height; |
478 | local_dom = hgv.height; |
477 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
479 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
478 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
480 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
479 | 481 | ||
480 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
482 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
481 | own_value = qMax(own_value, 0.02); |
483 | own_value = qMax(own_value, 0.02); |
482 | value = *grid_value++ / own_value; // remove impact of focal tree |
484 | value = *grid_value++ / own_value; // remove impact of focal tree |
483 | // additional punishment if pixel is outside:
|
485 | // additional punishment if pixel is outside:
|
484 | if (hgv.isForestOutside()) |
486 | if (hgv.isForestOutside()) |
485 | value *= outside_area_factor; |
487 | value *= outside_area_factor; |
486 | 488 | ||
487 | //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset);
|
489 | //qDebug() << x << y << local_dom << z << z_zstar << own_value << value << *(grid_value-1) << (*reader)(x,y) << mStamp->offsetValue(x,y,d_offset);
|
488 | //if (value>0.)
|
490 | //if (value>0.)
|
489 | sum += value * (*reader)(x,y); |
491 | sum += value * (*reader)(x,y); |
490 | 492 | ||
491 | }
|
493 | }
|
492 | }
|
494 | }
|
493 | mLRI = static_cast<float>( sum ); |
495 | mLRI = static_cast<float>( sum ); |
494 | // LRI correction...
|
496 | // LRI correction...
|
495 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
497 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
496 | if (hrel<1.) |
498 | if (hrel<1.) |
497 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
499 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
498 | 500 | ||
499 | 501 | ||
500 | if (mLRI > 1.) |
502 | if (mLRI > 1.) |
501 | mLRI = 1.; |
503 | mLRI = 1.; |
502 | 504 | ||
503 | // Finally, add LRI of this Tree to the ResourceUnit!
|
505 | // Finally, add LRI of this Tree to the ResourceUnit!
|
504 | mRU->addWLA(mLeafArea, mLRI); |
506 | mRU->addWLA(mLeafArea, mLRI); |
505 | 507 | ||
506 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
508 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
507 | }
|
509 | }
|
508 | 510 | ||
509 | /// Torus version of read stamp (glued edges)
|
511 | /// Torus version of read stamp (glued edges)
|
510 | void Tree::readLIF_torus() |
512 | void Tree::readLIF_torus() |
511 | {
|
513 | {
|
512 | if (!mStamp) |
514 | if (!mStamp) |
513 | return; |
515 | return; |
514 | const Stamp *reader = mStamp->reader(); |
516 | const Stamp *reader = mStamp->reader(); |
515 | if (!reader) |
517 | if (!reader) |
516 | return; |
518 | return; |
517 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
519 | int bufferOffset = mGrid->indexAt(QPointF(0.,0.)).x(); // offset of buffer |
518 | 520 | ||
519 | QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
521 | QPoint pos_reader = QPoint((mPositionIndex.x()-bufferOffset)%cPxPerRU + bufferOffset, |
520 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
522 | (mPositionIndex.y()-bufferOffset)%cPxPerRU + bufferOffset); // offset within the ha |
521 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index |
523 | QPoint ru_offset = QPoint(mPositionIndex.x() - pos_reader.x(), mPositionIndex.y() - pos_reader.y()); // offset of the corner of the resource index |
522 | 524 | ||
523 | int offset_reader = reader->offset(); |
525 | int offset_reader = reader->offset(); |
524 | int offset_writer = mStamp->offset(); |
526 | int offset_writer = mStamp->offset(); |
525 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
527 | int d_offset = offset_writer - offset_reader; // offset on the *stamp* to the crown-cells |
526 | 528 | ||
527 | pos_reader-=QPoint(offset_reader, offset_reader); |
529 | pos_reader-=QPoint(offset_reader, offset_reader); |
528 | 530 | ||
529 | float local_dom; |
531 | float local_dom; |
530 | 532 | ||
531 | int x,y; |
533 | int x,y; |
532 | double sum=0.; |
534 | double sum=0.; |
533 | double value, own_value; |
535 | double value, own_value; |
534 | float *grid_value; |
536 | float *grid_value; |
535 | float z, z_zstar; |
537 | float z, z_zstar; |
536 | int reader_size = reader->size(); |
538 | int reader_size = reader->size(); |
537 | int rx = pos_reader.x(); |
539 | int rx = pos_reader.x(); |
538 | int ry = pos_reader.y(); |
540 | int ry = pos_reader.y(); |
539 | int xt, yt; // wrapped coords |
541 | int xt, yt; // wrapped coords |
540 | 542 | ||
541 | for (y=0;y<reader_size; ++y) { |
543 | for (y=0;y<reader_size; ++y) { |
542 | yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y()); |
544 | yt = torusIndex(ry+y,cPxPerRU, bufferOffset, ru_offset.y()); |
543 | for (x=0;x<reader_size;++x) { |
545 | for (x=0;x<reader_size;++x) { |
544 | xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x()); |
546 | xt = torusIndex(rx+x,cPxPerRU, bufferOffset, ru_offset.x()); |
545 | grid_value = mGrid->ptr(xt,yt); |
547 | grid_value = mGrid->ptr(xt,yt); |
546 | 548 | ||
547 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
549 | local_dom = mHeightGrid->valueAtIndex(xt/cPxPerHeight, yt/cPxPerHeight).height; // ry: gets ++ed in outer loop, rx not |
548 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
550 | z = std::max(mHeight - reader->distanceToCenter(x,y), 0.f); // distance to center = height (45 degree line) |
549 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
551 | z_zstar = (z>=local_dom)?1.f:z/local_dom; |
550 | 552 | ||
551 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
553 | own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity * z_zstar; |
552 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
554 | // old: own_value = 1. - mStamp->offsetValue(x,y,d_offset)*mOpacity / local_dom; // old: dom_height;
|
553 | own_value = qMax(own_value, 0.02); |
555 | own_value = qMax(own_value, 0.02); |
554 | value = *grid_value++ / own_value; // remove impact of focal tree |
556 | value = *grid_value++ / own_value; // remove impact of focal tree |
555 | 557 | ||
556 | // debug for one tree in HJA
|
558 | // debug for one tree in HJA
|
557 | //if (id()==178020)
|
559 | //if (id()==178020)
|
558 | // qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y);
|
560 | // qDebug() << x << y << xt << yt << *grid_value << local_dom << own_value << value << (*reader)(x,y);
|
559 | //if (_isnan(value))
|
561 | //if (_isnan(value))
|
560 | // qDebug() << "isnan" << id();
|
562 | // qDebug() << "isnan" << id();
|
561 | if (value * (*reader)(x,y)>1.) |
563 | if (value * (*reader)(x,y)>1.) |
562 | qDebug() << "LIFTorus: value>1."; |
564 | qDebug() << "LIFTorus: value>1."; |
563 | sum += value * (*reader)(x,y); |
565 | sum += value * (*reader)(x,y); |
564 | 566 | ||
565 | //} // isIndexValid
|
567 | //} // isIndexValid
|
566 | }
|
568 | }
|
567 | }
|
569 | }
|
568 | mLRI = static_cast<float>( sum ); |
570 | mLRI = static_cast<float>( sum ); |
569 | 571 | ||
570 | // LRI correction...
|
572 | // LRI correction...
|
571 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
573 | double hrel = mHeight / mHeightGrid->valueAtIndex(mPositionIndex.x()/cPxPerHeight, mPositionIndex.y()/cPxPerHeight).height; |
572 | if (hrel<1.) |
574 | if (hrel<1.) |
573 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
575 | mLRI = static_cast<float>( species()->speciesSet()->LRIcorrection(mLRI, hrel) ); |
574 | 576 | ||
575 | 577 | ||
576 | if (isnan(mLRI)) { |
578 | if (isnan(mLRI)) { |
577 | qDebug() << "LRI invalid (nan)!" << id(); |
579 | qDebug() << "LRI invalid (nan)!" << id(); |
578 | mLRI=0.; |
580 | mLRI=0.; |
579 | //qDebug() << reader->dump();
|
581 | //qDebug() << reader->dump();
|
580 | }
|
582 | }
|
581 | if (mLRI > 1.) |
583 | if (mLRI > 1.) |
582 | mLRI = 1.; |
584 | mLRI = 1.; |
583 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
585 | //qDebug() << "Tree #"<< id() << "value" << sum << "Impact" << mImpact;
|
584 | 586 | ||
585 | // Finally, add LRI of this Tree to the ResourceUnit!
|
587 | // Finally, add LRI of this Tree to the ResourceUnit!
|
586 | mRU->addWLA(mLeafArea, mLRI); |
588 | mRU->addWLA(mLeafArea, mLRI); |
587 | }
|
589 | }
|
588 | 590 | ||
589 | 591 | ||
590 | void Tree::resetStatistics() |
592 | void Tree::resetStatistics() |
591 | {
|
593 | {
|
592 | m_statPrint=0; |
594 | m_statPrint=0; |
593 | m_statCreated=0; |
595 | m_statCreated=0; |
594 | m_statAboveZ=0; |
596 | m_statAboveZ=0; |
595 | m_nextId=1; |
597 | m_nextId=1; |
596 | }
|
598 | }
|
597 | 599 | ||
598 | #ifdef ALT_TREE_MORTALITY
|
600 | #ifdef ALT_TREE_MORTALITY
|
599 | void Tree::mortalityParams(double dbh_inc_threshold, int stress_years, double stress_mort_prob) |
601 | void Tree::mortalityParams(double dbh_inc_threshold, int stress_years, double stress_mort_prob) |
600 | {
|
602 | {
|
601 | _stress_threshold = dbh_inc_threshold; |
603 | _stress_threshold = dbh_inc_threshold; |
602 | _stress_years = stress_years; |
604 | _stress_years = stress_years; |
603 | _stress_death_prob = stress_mort_prob; |
605 | _stress_death_prob = stress_mort_prob; |
604 | qDebug() << "Alternative Mortality enabled: threshold" << dbh_inc_threshold << ", years:" << _stress_years << ", level:" << _stress_death_prob; |
606 | qDebug() << "Alternative Mortality enabled: threshold" << dbh_inc_threshold << ", years:" << _stress_years << ", level:" << _stress_death_prob; |
605 | }
|
607 | }
|
606 | #endif
|
608 | #endif
|
607 | 609 | ||
608 | void Tree::calcLightResponse() |
610 | void Tree::calcLightResponse() |
609 | {
|
611 | {
|
610 | // calculate a light response from lri:
|
612 | // calculate a light response from lri:
|
611 | // http://iland.boku.ac.at/individual+tree+light+availability
|
613 | // http://iland.boku.ac.at/individual+tree+light+availability
|
612 | double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.); // Eq. (3) |
614 | double lri = limit(mLRI * mRU->LRImodifier(), 0., 1.); // Eq. (3) |
613 | mLightResponse = static_cast<float>( mSpecies->lightResponse(lri) ); // Eq. (4) |
615 | mLightResponse = static_cast<float>( mSpecies->lightResponse(lri) ); // Eq. (4) |
614 | mRU->addLR(mLeafArea, mLightResponse); |
616 | mRU->addLR(mLeafArea, mLightResponse); |
615 | 617 | ||
616 | }
|
618 | }
|
617 | 619 | ||
618 | //////////////////////////////////////////////////
|
620 | //////////////////////////////////////////////////
|
619 | //// Growth Functions
|
621 | //// Growth Functions
|
620 | //////////////////////////////////////////////////
|
622 | //////////////////////////////////////////////////
|
621 | 623 | ||
622 | /** grow() is the main function of the yearly tree growth.
|
624 | /** grow() is the main function of the yearly tree growth.
|
623 | The main steps are:
|
625 | The main steps are:
|
624 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability
|
626 | - Production of GPP/NPP @sa http://iland.boku.ac.at/primary+production http://iland.boku.ac.at/individual+tree+light+availability
|
625 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation
|
627 | - Partitioning of NPP to biomass compartments of the tree @sa http://iland.boku.ac.at/allocation
|
626 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???)
|
628 | - Growth of the stem http://iland.boku.ac.at/stem+growth (???)
|
627 | Further activties: * the age of the tree is increased
|
629 | Further activties: * the age of the tree is increased
|
628 | * the mortality sub routine is executed
|
630 | * the mortality sub routine is executed
|
629 | * seeds are produced */
|
631 | * seeds are produced */
|
630 | void Tree::grow() |
632 | void Tree::grow() |
631 | {
|
633 | {
|
632 | TreeGrowthData d;
|
634 | TreeGrowthData d;
|
633 | mAge++; // increase age |
635 | mAge++; // increase age |
634 | // step 1: get "interception area" of the tree individual [m2]
|
636 | // step 1: get "interception area" of the tree individual [m2]
|
635 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert)
|
637 | // the sum of all area of all trees of a unit equal the total stocked area * interception_factor(Beer-Lambert)
|
636 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
638 | double effective_area = mRU->interceptedArea(mLeafArea, mLightResponse); |
637 | 639 | ||
638 | // step 2: calculate GPP of the tree based
|
640 | // step 2: calculate GPP of the tree based
|
639 | // (1) get the amount of GPP for a "unit area" of the tree species
|
641 | // (1) get the amount of GPP for a "unit area" of the tree species
|
640 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
642 | double raw_gpp_per_area = mRU->resourceUnitSpecies(species()).prod3PG().GPPperArea(); |
641 | // (2) GPP (without aging-effect) in kg Biomass / year
|
643 | // (2) GPP (without aging-effect) in kg Biomass / year
|
642 | double raw_gpp = raw_gpp_per_area * effective_area; |
644 | double raw_gpp = raw_gpp_per_area * effective_area; |
643 | 645 | ||
644 | // apply aging according to the state of the individuum
|
646 | // apply aging according to the state of the individuum
|
645 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
647 | const double aging_factor = mSpecies->aging(mHeight, mAge); |
646 | mRU->addTreeAging(mLeafArea, aging_factor); |
648 | mRU->addTreeAging(mLeafArea, aging_factor); |
647 | double gpp = raw_gpp * aging_factor; // |
649 | double gpp = raw_gpp * aging_factor; // |
648 | d.NPP = gpp * cAutotrophicRespiration; // respiration loss (0.47), cf. Waring et al 1998. |
650 | d.NPP = gpp * cAutotrophicRespiration; // respiration loss (0.47), cf. Waring et al 1998. |
649 | 651 | ||
650 | //DBGMODE(
|
652 | //DBGMODE(
|
651 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
653 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeNPP) && isDebugging()) { |
652 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
654 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeNPP); |
653 | dumpList(out); // add tree headers |
655 | dumpList(out); // add tree headers |
654 | out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
656 | out << mLRI * mRU->LRImodifier() << mLightResponse << effective_area << raw_gpp << gpp << d.NPP << aging_factor; |
655 | }
|
657 | }
|
656 | //); // DBGMODE()
|
658 | //); // DBGMODE()
|
657 | if (d.NPP>0.) |
659 | if (d.NPP>0.) |
658 | partitioning(d); // split npp to compartments and grow (diameter, height) |
660 | partitioning(d); // split npp to compartments and grow (diameter, height) |
659 | 661 | ||
660 | // mortality
|
662 | // mortality
|
661 | #ifdef ALT_TREE_MORTALITY
|
663 | #ifdef ALT_TREE_MORTALITY
|
662 | // alternative variant of tree mortality (note: mStrssIndex used otherwise)
|
664 | // alternative variant of tree mortality (note: mStrssIndex used otherwise)
|
663 | altMortality(d); |
665 | altMortality(d); |
664 | 666 | ||
665 | #else
|
667 | #else
|
666 | if (Model::settings().mortalityEnabled) |
668 | if (Model::settings().mortalityEnabled) |
667 | mortality(d); |
669 | mortality(d); |
668 | 670 | ||
669 | mStressIndex = d.stress_index; |
671 | mStressIndex = d.stress_index; |
670 | #endif
|
672 | #endif
|
671 | 673 | ||
672 | if (!isDead()) |
674 | if (!isDead()) |
673 | mRU->resourceUnitSpecies(species()).statistics().add(this, &d); |
675 | mRU->resourceUnitSpecies(species()).statistics().add(this, &d); |
674 | 676 | ||
675 | // regeneration
|
677 | // regeneration
|
676 | mSpecies->seedProduction(mAge, mHeight, mPositionIndex); |
678 | mSpecies->seedProduction(mAge, mHeight, mPositionIndex); |
677 | 679 | ||
678 | }
|
680 | }
|
679 | 681 | ||
680 | /** partitioning of this years assimilates (NPP) to biomass compartments.
|
682 | /** partitioning of this years assimilates (NPP) to biomass compartments.
|
681 | Conceptionally, the algorithm is based on Duursma, 2007.
|
683 | Conceptionally, the algorithm is based on Duursma, 2007.
|
682 | @sa http://iland.boku.ac.at/allocation */
|
684 | @sa http://iland.boku.ac.at/allocation */
|
683 | inline void Tree::partitioning(TreeGrowthData &d) |
685 | inline void Tree::partitioning(TreeGrowthData &d) |
684 | {
|
686 | {
|
685 | double npp = d.NPP; |
687 | double npp = d.NPP; |
686 | // add content of reserve pool
|
688 | // add content of reserve pool
|
687 | npp += mNPPReserve; |
689 | npp += mNPPReserve; |
688 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
690 | const double foliage_mass_allo = species()->biomassFoliage(mDbh); |
689 | const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio()); |
691 | const double reserve_size = foliage_mass_allo * (1. + mSpecies->finerootFoliageRatio()); |
690 | double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100% |
692 | double refill_reserve = qMin(reserve_size, (1. + mSpecies->finerootFoliageRatio())*mFoliageMass); // not always try to refill reserve 100% |
691 | 693 | ||
692 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
694 | double apct_wood, apct_root, apct_foliage; // allocation percentages (sum=1) (eta) |
693 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
695 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
694 | // turnover rates
|
696 | // turnover rates
|
695 | const double &to_fol = species()->turnoverLeaf(); |
697 | const double &to_fol = species()->turnoverLeaf(); |
696 | const double &to_root = species()->turnoverRoot(); |
698 | const double &to_root = species()->turnoverRoot(); |
697 | // the turnover rate of wood depends on the size of the reserve pool:
|
699 | // the turnover rate of wood depends on the size of the reserve pool:
|
698 | 700 | ||
699 | 701 | ||
700 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
702 | double to_wood = refill_reserve / (mWoodyMass + refill_reserve); |
701 | 703 | ||
702 | apct_root = rus.prod3PG().rootFraction(); |
704 | apct_root = rus.prod3PG().rootFraction(); |
703 | d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots |
705 | d.NPP_above = d.NPP * (1. - apct_root); // aboveground: total NPP - fraction to roots |
704 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage) |
706 | double b_wf = species()->allometricRatio_wf(); // ratio of allometric exponents (b_woody / b_foliage) |
705 | 707 | ||
706 | // Duursma 2007, Eq. (20)
|
708 | // Duursma 2007, Eq. (20)
|
707 | apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf ); |
709 | apct_wood = (foliage_mass_allo*to_wood/npp + b_wf*(1.-apct_root) - b_wf*foliage_mass_allo*to_fol/npp) / ( foliage_mass_allo/mWoodyMass + b_wf ); |
708 | 710 | ||
709 | apct_wood = limit(apct_wood, 0., 1.-apct_root); |
711 | apct_wood = limit(apct_wood, 0., 1.-apct_root); |
710 | 712 | ||
711 | apct_foliage = 1. - apct_root - apct_wood; |
713 | apct_foliage = 1. - apct_root - apct_wood; |
712 | 714 | ||
713 | 715 | ||
714 | DBGMODE(
|
716 | DBGMODE(
|
715 | if (apct_foliage<0 || apct_wood<0) |
717 | if (apct_foliage<0 || apct_wood<0) |
716 | qDebug() << "transfer to foliage or wood < 0"; |
718 | qDebug() << "transfer to foliage or wood < 0"; |
717 | if (npp<0) |
719 | if (npp<0) |
718 | qDebug() << "NPP < 0"; |
720 | qDebug() << "NPP < 0"; |
719 | ); |
721 | ); |
720 | 722 | ||
721 | // Change of biomass compartments
|
723 | // Change of biomass compartments
|
722 | double sen_root = mFineRootMass * to_root; |
724 | double sen_root = mFineRootMass * to_root; |
723 | double sen_foliage = mFoliageMass * to_fol; |
725 | double sen_foliage = mFoliageMass * to_fol; |
724 | if (ru()->snag()) |
726 | if (ru()->snag()) |
725 | ru()->snag()->addTurnoverLitter(this->species(), sen_foliage, sen_root); |
727 | ru()->snag()->addTurnoverLitter(this->species(), sen_foliage, sen_root); |
726 | 728 | ||
727 | // Roots
|
729 | // Roots
|
728 | // http://iland.boku.ac.at/allocation#belowground_NPP
|
730 | // http://iland.boku.ac.at/allocation#belowground_NPP
|
729 | mFineRootMass -= sen_root; // reduce only fine root pool |
731 | mFineRootMass -= sen_root; // reduce only fine root pool |
730 | double delta_root = apct_root * npp; |
732 | double delta_root = apct_root * npp; |
731 | // 1st, refill the fine root pool
|
733 | // 1st, refill the fine root pool
|
732 | double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass; |
734 | double fineroot_miss = mFoliageMass * mSpecies->finerootFoliageRatio() - mFineRootMass; |
733 | if (fineroot_miss>0.){ |
735 | if (fineroot_miss>0.){ |
734 | double delta_fineroot = qMin(fineroot_miss, delta_root); |
736 | double delta_fineroot = qMin(fineroot_miss, delta_root); |
735 | mFineRootMass += delta_fineroot; |
737 | mFineRootMass += delta_fineroot; |
736 | delta_root -= delta_fineroot; |
738 | delta_root -= delta_fineroot; |
737 | }
|
739 | }
|
738 | // 2nd, the rest of NPP allocated to roots go to coarse roots
|
740 | // 2nd, the rest of NPP allocated to roots go to coarse roots
|
739 | double max_coarse_root = species()->biomassRoot(mDbh); |
741 | double max_coarse_root = species()->biomassRoot(mDbh); |
740 | mCoarseRootMass += delta_root; |
742 | mCoarseRootMass += delta_root; |
741 | if (mCoarseRootMass > max_coarse_root) { |
743 | if (mCoarseRootMass > max_coarse_root) { |
742 | // if the coarse root pool exceeds the value given by the allometry, then the
|
744 | // if the coarse root pool exceeds the value given by the allometry, then the
|
743 | // surplus is accounted as turnover
|
745 | // surplus is accounted as turnover
|
744 | if (ru()->snag()) |
746 | if (ru()->snag()) |
745 | ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root); |
747 | ru()->snag()->addTurnoverWood(species(), mCoarseRootMass-max_coarse_root); |
746 | 748 | ||
747 | mCoarseRootMass = static_cast<float>( max_coarse_root ); |
749 | mCoarseRootMass = static_cast<float>( max_coarse_root ); |
748 | }
|
750 | }
|
749 | 751 | ||
750 | // Foliage
|
752 | // Foliage
|
751 | double delta_foliage = apct_foliage * npp - sen_foliage; |
753 | double delta_foliage = apct_foliage * npp - sen_foliage; |
752 | mFoliageMass += delta_foliage; |
754 | mFoliageMass += delta_foliage; |
753 | if (isnan(mFoliageMass)) |
755 | if (isnan(mFoliageMass)) |
754 | qDebug() << "foliage mass invalid!"; |
756 | qDebug() << "foliage mass invalid!"; |
755 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
757 | if (mFoliageMass<0.) mFoliageMass=0.; // limit to zero |
756 | 758 | ||
757 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); // update leaf area |
759 | mLeafArea = static_cast<float>( mFoliageMass * species()->specificLeafArea() ); // update leaf area |
758 | 760 | ||
759 | // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild,
|
761 | // stress index: different varaints at denominator: to_fol*foliage_mass = leafmass to rebuild,
|
760 | // foliage_mass_allo: simply higher chance for stress
|
762 | // foliage_mass_allo: simply higher chance for stress
|
761 | // note: npp = NPP + reserve (see above)
|
763 | // note: npp = NPP + reserve (see above)
|
762 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.); |
764 | d.stress_index =qMax(1. - (npp) / ( to_fol*foliage_mass_allo + to_root*foliage_mass_allo*species()->finerootFoliageRatio() + reserve_size), 0.); |
763 | 765 | ||
764 | // Woody compartments
|
766 | // Woody compartments
|
765 | // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth
|
767 | // see also: http://iland.boku.ac.at/allocation#reserve_and_allocation_to_stem_growth
|
766 | // (1) transfer to reserve pool
|
768 | // (1) transfer to reserve pool
|
767 | double gross_woody = apct_wood * npp; |
769 | double gross_woody = apct_wood * npp; |
768 | double to_reserve = qMin(reserve_size, gross_woody); |
770 | double to_reserve = qMin(reserve_size, gross_woody); |
769 | mNPPReserve = static_cast<float>( to_reserve ); |
771 | mNPPReserve = static_cast<float>( to_reserve ); |
770 | double net_woody = gross_woody - to_reserve; |
772 | double net_woody = gross_woody - to_reserve; |
771 | double net_stem = 0.; |
773 | double net_stem = 0.; |
772 | mDbhDelta = 0.; |
774 | mDbhDelta = 0.; |
773 | 775 | ||
774 | 776 | ||
775 | if (net_woody > 0.) { |
777 | if (net_woody > 0.) { |
776 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter)
|
778 | // (2) calculate part of increment that is dedicated to the stem (which is a function of diameter)
|
777 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
779 | net_stem = net_woody * species()->allometricFractionStem(mDbh); |
778 | d.NPP_stem = net_stem; |
780 | d.NPP_stem = net_stem; |
779 | mWoodyMass += net_woody; |
781 | mWoodyMass += net_woody; |
780 | // (3) growth of diameter and height baseed on net stem increment
|
782 | // (3) growth of diameter and height baseed on net stem increment
|
781 | grow_diameter(d); |
783 | grow_diameter(d); |
782 | }
|
784 | }
|
783 | 785 | ||
784 | //DBGMODE(
|
786 | //DBGMODE(
|
785 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
787 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreePartition) |
786 | && isDebugging() ) { |
788 | && isDebugging() ) { |
787 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
789 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreePartition); |
788 | dumpList(out); // add tree headers |
790 | dumpList(out); // add tree headers |
789 | out << npp << apct_foliage << apct_wood << apct_root |
791 | out << npp << apct_foliage << apct_wood << apct_root |
790 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index; |
792 | << delta_foliage << net_woody << delta_root << mNPPReserve << net_stem << d.stress_index; |
791 | }
|
793 | }
|
792 | 794 | ||
793 | //); // DBGMODE()
|
795 | //); // DBGMODE()
|
794 | DBGMODE(
|
796 | DBGMODE(
|
795 | if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000 |
797 | if (mWoodyMass<0. || mWoodyMass>50000 || mFoliageMass<0. || mFoliageMass>2000. || mCoarseRootMass<0. || mCoarseRootMass>30000 |
796 | || mNPPReserve>4000.) { |
798 | || mNPPReserve>4000.) { |
797 | qDebug() << "Tree:partitioning: invalid or unlikely pools."; |
799 | qDebug() << "Tree:partitioning: invalid or unlikely pools."; |
798 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
800 | qDebug() << GlobalSettings::instance()->debugListCaptions(GlobalSettings::DebugOutputs(0)); |
799 | DebugList dbg; dumpList(dbg); |
801 | DebugList dbg; dumpList(dbg); |
800 | qDebug() << dbg; |
802 | qDebug() << dbg; |
801 | } ); |
803 | } ); |
802 | 804 | ||
803 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump()
|
805 | /*DBG_IF_X(mId == 1 , "Tree::partitioning", "dump", dump()
|
804 | + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8")
|
806 | + QString("npp %1 npp_reserve %9 sen_fol %2 sen_stem %3 sen_root %4 net_fol %5 net_stem %6 net_root %7 to_reserve %8")
|
805 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root)
|
807 | .arg(npp).arg(senescence_foliage).arg(senescence_stem).arg(senescence_root)
|
806 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/
|
808 | .arg(net_foliage).arg(net_stem).arg(net_root).arg(to_reserve).arg(mNPPReserve) );*/
|
807 | 809 | ||
808 | }
|
810 | }
|
809 | 811 | ||
810 | 812 | ||
811 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp).
|
813 | /** Determination of diamter and height growth based on increment of the stem mass (@p net_stem_npp).
|
812 | Refer to XXX for equations and variables.
|
814 | Refer to XXX for equations and variables.
|
813 | This function updates the dbh and height of the tree.
|
815 | This function updates the dbh and height of the tree.
|
814 | The equations are based on dbh in meters! */
|
816 | The equations are based on dbh in meters! */
|
815 | inline void Tree::grow_diameter(TreeGrowthData &d) |
817 | inline void Tree::grow_diameter(TreeGrowthData &d) |
816 | {
|
818 | {
|
817 | // determine dh-ratio of increment
|
819 | // determine dh-ratio of increment
|
818 | // height increment is a function of light competition:
|
820 | // height increment is a function of light competition:
|
819 | double hd_growth = relative_height_growth(); // hd of height growth |
821 | double hd_growth = relative_height_growth(); // hd of height growth |
820 | double d_m = mDbh / 100.; // current diameter in [m] |
822 | double d_m = mDbh / 100.; // current diameter in [m] |
821 | double net_stem_npp = d.NPP_stem; |
823 | double net_stem_npp = d.NPP_stem; |
822 | 824 | ||
823 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
825 | const double d_delta_m = mDbhDelta / 100.; // increment of last year in [m] |
824 | 826 | ||
825 | const double mass_factor = species()->volumeFactor() * species()->density(); |
827 | const double mass_factor = species()->volumeFactor() * species()->density(); |
826 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
828 | double stem_mass = mass_factor * d_m*d_m * mHeight; // result: kg, dbh[cm], h[meter] |
827 | 829 | ||
828 | // factor is in diameter increment per NPP [m/kg]
|
830 | // factor is in diameter increment per NPP [m/kg]
|
829 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
831 | double factor_diameter = 1. / ( mass_factor * (d_m + d_delta_m)*(d_m + d_delta_m) * ( 2. * mHeight/d_m + hd_growth) ); |
830 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
832 | double delta_d_estimate = factor_diameter * net_stem_npp; // estimated dbh-inc using last years increment |
831 | 833 | ||
832 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9)
|
834 | // using that dbh-increment we estimate a stem-mass-increment and the residual (Eq. 9)
|
833 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
835 | double stem_estimate = mass_factor * (d_m + delta_d_estimate)*(d_m + delta_d_estimate)*(mHeight + delta_d_estimate*hd_growth); |
834 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
836 | double stem_residual = stem_estimate - (stem_mass + net_stem_npp); |
835 | 837 | ||
836 | // the final increment is then:
|
838 | // the final increment is then:
|
837 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
839 | double d_increment = factor_diameter * (net_stem_npp - stem_residual); // Eq. (11) |
838 | double res_final = 0.; |
840 | double res_final = 0.; |
839 | if (fabs(stem_residual) > 1.) { |
841 | if (fabs(stem_residual) > 1.) { |
840 | 842 | ||
841 | // calculate final residual in stem
|
843 | // calculate final residual in stem
|
842 | res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
844 | res_final = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)); |
843 | if (fabs(res_final)>1.) { |
845 | if (fabs(res_final)>1.) { |
844 | // for large errors in stem biomass due to errors in diameter increment (> 1kg), we solve the increment iteratively.
|
846 | // for large errors in stem biomass due to errors in diameter increment (> 1kg), we solve the increment iteratively.
|
845 | // first, increase increment with constant step until we overestimate the first time
|
847 | // first, increase increment with constant step until we overestimate the first time
|
846 | // then,
|
848 | // then,
|
847 | d_increment = 0.02; // start with 2cm increment |
849 | d_increment = 0.02; // start with 2cm increment |
848 | bool reached_error = false; |
850 | bool reached_error = false; |
849 | double step=0.01; // step-width 1cm |
851 | double step=0.01; // step-width 1cm |
850 | double est_stem; |
852 | double est_stem; |
851 | do { |
853 | do { |
852 | est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment |
854 | est_stem = mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth); // estimate with current increment |
853 | stem_residual = est_stem - (stem_mass + net_stem_npp); |
855 | stem_residual = est_stem - (stem_mass + net_stem_npp); |
854 | 856 | ||
855 | if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg |
857 | if (fabs(stem_residual) <1.) // finished, if stem residual below 1kg |
856 | break; |
858 | break; |
857 | if (stem_residual > 0.) { |
859 | if (stem_residual > 0.) { |
858 | d_increment -= step; |
860 | d_increment -= step; |
859 | reached_error=true; |
861 | reached_error=true; |
860 | } else { |
862 | } else { |
861 | d_increment += step; |
863 | d_increment += step; |
862 | }
|
864 | }
|
863 | if (reached_error) |
865 | if (reached_error) |
864 | step /= 2.; |
866 | step /= 2.; |
865 | } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm |
867 | } while (step>0.00001); // continue until diameter "accuracy" falls below 1/100mm |
866 | }
|
868 | }
|
867 | }
|
869 | }
|
868 | 870 | ||
869 | if (d_increment<0.f) |
871 | if (d_increment<0.f) |
870 | qDebug() << "Tree::grow_diameter: d_inc < 0."; |
872 | qDebug() << "Tree::grow_diameter: d_inc < 0."; |
871 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
873 | DBG_IF_X(d_increment<0. || d_increment>0.1, "Tree::grow_dimater", "increment out of range.", dump() |
872 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
874 | + QString("\nhdz %1 factor_diameter %2 stem_residual %3 delta_d_estimate %4 d_increment %5 final residual(kg) %6") |
873 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
875 | .arg(hd_growth).arg(factor_diameter).arg(stem_residual).arg(delta_d_estimate).arg(d_increment) |
874 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
876 | .arg( mass_factor * (mDbh + d_increment)*(mDbh + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp)) )); |
875 | 877 | ||
876 | //DBGMODE(
|
878 | //DBGMODE(
|
877 | // do not calculate res_final twice if already done
|
879 | // do not calculate res_final twice if already done
|
878 | DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
880 | DBG_IF_X( (res_final==0.?fabs(mass_factor * (d_m + d_increment)*(d_m + d_increment)*(mHeight + d_increment*hd_growth)-((stem_mass + net_stem_npp))):res_final) > 1, "Tree::grow_diameter", "final residual stem estimate > 1kg", dump()); |
879 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump()); |
881 | DBG_IF_X(d_increment > 10. || d_increment*hd_growth >10., "Tree::grow_diameter", "growth out of bound:",QString("d-increment %1 h-increment %2 ").arg(d_increment).arg(d_increment*hd_growth/100.) + dump()); |
880 | 882 | ||
881 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
883 | if (GlobalSettings::instance()->isDebugEnabled(GlobalSettings::dTreeGrowth) && isDebugging() ) { |
882 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
884 | DebugList &out = GlobalSettings::instance()->debugList(mId, GlobalSettings::dTreeGrowth); |
883 | dumpList(out); // add tree headers |
885 | dumpList(out); // add tree headers |
884 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
886 | out << net_stem_npp << stem_mass << hd_growth << factor_diameter << delta_d_estimate*100 << d_increment*100; |
885 | }
|
887 | }
|
886 | 888 | ||
887 | //); // DBGMODE()
|
889 | //); // DBGMODE()
|
888 | 890 | ||
889 | d_increment = qMax(d_increment, 0.); |
891 | d_increment = qMax(d_increment, 0.); |
890 | 892 | ||
891 | // update state variables
|
893 | // update state variables
|
892 | mDbh += d_increment*100.f; // convert from [m] to [cm] |
894 | mDbh += d_increment*100.f; // convert from [m] to [cm] |
893 | mDbhDelta = static_cast<float>( d_increment*100. ); // save for next year's growth |
895 | mDbhDelta = static_cast<float>( d_increment*100. ); // save for next year's growth |
894 | mHeight += d_increment * hd_growth; |
896 | mHeight += d_increment * hd_growth; |
895 | 897 | ||
896 | // update state of LIP stamp and opacity
|
898 | // update state of LIP stamp and opacity
|
897 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
899 | mStamp = species()->stamp(mDbh, mHeight); // get new stamp for updated dimensions |
898 | // calculate the CrownFactor which reflects the opacity of the crown
|
900 | // calculate the CrownFactor which reflects the opacity of the crown
|
899 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
901 | const double k=Model::settings().lightExtinctionCoefficientOpacity; |
900 | mOpacity = static_cast<float>( 1. - exp(-k * mLeafArea / mStamp->crownArea()) ); |
902 | mOpacity = static_cast<float>( 1. - exp(-k * mLeafArea / mStamp->crownArea()) ); |
901 | 903 | ||
902 | }
|
904 | }
|
903 | 905 | ||
904 | 906 | ||
905 | /// return the HD ratio of this year's increment based on the light status.
|
907 | /// return the HD ratio of this year's increment based on the light status.
|
906 | inline double Tree::relative_height_growth() |
908 | inline double Tree::relative_height_growth() |
907 | {
|
909 | {
|
908 | double hd_low, hd_high; |
910 | double hd_low, hd_high; |
909 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
911 | mSpecies->hdRange(mDbh, hd_low, hd_high); |
910 | 912 | ||
911 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
913 | DBG_IF_X(hd_low>hd_high, "Tree::relative_height_growth", "hd low higher dann hd_high for ", dump()); |
912 | DBG_IF_X(hd_low < 10 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high)); |
914 | DBG_IF_X(hd_low < 10 || hd_high>250, "Tree::relative_height_growth", "hd out of range ", dump() + QString(" hd-low: %1 hd-high: %2").arg(hd_low).arg(hd_high)); |
913 | 915 | ||
914 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees)
|
916 | // scale according to LRI: if receiving much light (LRI=1), the result is hd_low (for open grown trees)
|
915 | // use the corrected LRI (see tracker#11)
|
917 | // use the corrected LRI (see tracker#11)
|
916 | double lri = limit(mLRI * mRU->LRImodifier(),0.,1.); |
918 | double lri = limit(mLRI * mRU->LRImodifier(),0.,1.); |
917 | double hd_ratio = hd_high - (hd_high-hd_low)*lri; |
919 | double hd_ratio = hd_high - (hd_high-hd_low)*lri; |
918 | return hd_ratio; |
920 | return hd_ratio; |
919 | }
|
921 | }
|
920 | 922 | ||
921 | /** This function is called if a tree dies.
|
923 | /** This function is called if a tree dies.
|
922 | @sa ResourceUnit::cleanTreeList(), remove() */
|
924 | @sa ResourceUnit::cleanTreeList(), remove() */
|
923 | void Tree::die(TreeGrowthData *d) |
925 | void Tree::die(TreeGrowthData *d) |
924 | {
|
926 | {
|
925 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
927 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
926 | mRU->treeDied(); |
928 | mRU->treeDied(); |
927 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
929 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
928 | rus.statisticsDead().add(this, d); // add tree to statistics |
930 | rus.statisticsDead().add(this, d); // add tree to statistics |
929 | notifyTreeRemoved(TreeDeath); |
931 | notifyTreeRemoved(TreeDeath); |
930 | if (ru()->snag()) |
932 | if (ru()->snag()) |
931 | ru()->snag()->addMortality(this); |
933 | ru()->snag()->addMortality(this); |
932 | }
|
934 | }
|
933 | 935 | ||
934 | /// remove a tree (most likely due to harvest) from the system.
|
936 | /// remove a tree (most likely due to harvest) from the system.
|
935 | void Tree::remove(double removeFoliage, double removeBranch, double removeStem ) |
937 | void Tree::remove(double removeFoliage, double removeBranch, double removeStem ) |
936 | {
|
938 | {
|
937 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
939 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
938 | setIsHarvested(); |
940 | setIsHarvested(); |
939 | mRU->treeDied(); |
941 | mRU->treeDied(); |
940 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
942 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
941 | rus.statisticsMgmt().add(this, 0); |
943 | rus.statisticsMgmt().add(this, 0); |
942 | notifyTreeRemoved(TreeHarvest); |
944 | notifyTreeRemoved(TreeHarvest); |
943 | 945 | ||
944 | if (ru()->snag()) |
946 | if (ru()->snag()) |
945 | ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage); |
947 | ru()->snag()->addHarvest(this, removeStem, removeBranch, removeFoliage); |
946 | }
|
948 | }
|
947 | 949 | ||
948 | /// remove the tree due to an special event (disturbance)
|
950 | /// remove the tree due to an special event (disturbance)
|
949 | /// this is +- the same as die().
|
951 | /// this is +- the same as die().
|
950 | void Tree::removeDisturbance(const double stem_to_soil_fraction, const double stem_to_snag_fraction, const double branch_to_soil_fraction, const double branch_to_snag_fraction, const double foliage_to_soil_fraction) |
952 | void Tree::removeDisturbance(const double stem_to_soil_fraction, const double stem_to_snag_fraction, const double branch_to_soil_fraction, const double branch_to_snag_fraction, const double foliage_to_soil_fraction) |
951 | {
|
953 | {
|
952 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
954 | setFlag(Tree::TreeDead, true); // set flag that tree is dead |
953 | mRU->treeDied(); |
955 | mRU->treeDied(); |
954 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
956 | ResourceUnitSpecies &rus = mRU->resourceUnitSpecies(species()); |
955 | rus.statisticsDead().add(this, 0); |
957 | rus.statisticsDead().add(this, 0); |
956 | notifyTreeRemoved(TreeDisturbance); |
958 | notifyTreeRemoved(TreeDisturbance); |
957 | 959 | ||
958 | 960 | ||
959 | if (ru()->snag()) { |
961 | if (ru()->snag()) { |
960 | if (isHarvested()) { // if the tree is harvested, do the same as in normal tree harvest (but with default values) |
962 | if (isHarvested()) { // if the tree is harvested, do the same as in normal tree harvest (but with default values) |
961 | ru()->snag()->addHarvest(this, 1., 0., 0.); |
963 | ru()->snag()->addHarvest(this, 1., 0., 0.); |
962 | } else { |
964 | } else { |
963 | ru()->snag()->addDisturbance(this, stem_to_snag_fraction, stem_to_soil_fraction, branch_to_snag_fraction, branch_to_soil_fraction, foliage_to_soil_fraction); |
965 | ru()->snag()->addDisturbance(this, stem_to_snag_fraction, stem_to_soil_fraction, branch_to_snag_fraction, branch_to_soil_fraction, foliage_to_soil_fraction); |
964 | }
|
966 | }
|
965 | }
|
967 | }
|
966 | }
|
968 | }
|
967 | 969 | ||
968 | /// remove a part of the biomass of the tree, e.g. due to fire.
|
970 | /// remove a part of the biomass of the tree, e.g. due to fire.
|
969 | void Tree::removeBiomassOfTree(const double removeFoliageFraction, const double removeBranchFraction, const double removeStemFraction) |
971 | void Tree::removeBiomassOfTree(const double removeFoliageFraction, const double removeBranchFraction, const double removeStemFraction) |
970 | {
|
972 | {
|
971 | mFoliageMass *= 1. - removeFoliageFraction; |
973 | mFoliageMass *= 1. - removeFoliageFraction; |
972 | mWoodyMass *= (1. - removeStemFraction); |
974 | mWoodyMass *= (1. - removeStemFraction); |
973 | // we have a problem with the branches: this currently cannot be done properly!
|
975 | // we have a problem with the branches: this currently cannot be done properly!
|
974 | (void) removeBranchFraction; // silence warning |
976 | (void) removeBranchFraction; // silence warning |
975 | }
|
977 | }
|
976 | 978 | ||
977 | void Tree::setHeight(const float height) |
979 | void Tree::setHeight(const float height) |
978 | {
|
980 | {
|
979 | if (height<=0. || height>150.) |
981 | if (height<=0. || height>150.) |
980 | qWarning() << "trying to set tree height to invalid value:" << height << " for tree on RU:" << (mRU?mRU->boundingBox():QRect()); |
982 | qWarning() << "trying to set tree height to invalid value:" << height << " for tree on RU:" << (mRU?mRU->boundingBox():QRect()); |
981 | mHeight=height; |
983 | mHeight=height; |
982 | }
|
984 | }
|
983 | 985 | ||
984 | void Tree::mortality(TreeGrowthData &d) |
986 | void Tree::mortality(TreeGrowthData &d) |
985 | {
|
987 | {
|
986 | // death if leaf area is 0
|
988 | // death if leaf area is 0
|
987 | if (mFoliageMass<0.00001) |
989 | if (mFoliageMass<0.00001) |
988 | die(); |
990 | die(); |
989 | 991 | ||
990 | double p_death, p_stress, p_intrinsic; |
992 | double p_death, p_stress, p_intrinsic; |
991 | p_intrinsic = species()->deathProb_intrinsic(); |
993 | p_intrinsic = species()->deathProb_intrinsic(); |
992 | p_stress = species()->deathProb_stress(d.stress_index); |
994 | p_stress = species()->deathProb_stress(d.stress_index); |
993 | p_death = p_intrinsic + p_stress; |
995 | p_death = p_intrinsic + p_stress; |
994 | double p = drandom(); //0..1 |
996 | double p = drandom(); //0..1 |
995 | if (p<p_death) { |
997 | if (p<p_death) { |
996 | // die...
|
998 | // die...
|
997 | die(); |
999 | die(); |
998 | }
|
1000 | }
|
999 | }
|
1001 | }
|
1000 | 1002 | ||
1001 | #ifdef ALT_TREE_MORTALITY
|
1003 | #ifdef ALT_TREE_MORTALITY
|
1002 | void Tree::altMortality(TreeGrowthData &d) |
1004 | void Tree::altMortality(TreeGrowthData &d) |
1003 | {
|
1005 | {
|
1004 | // death if leaf area is 0
|
1006 | // death if leaf area is 0
|
1005 | if (mFoliageMass<0.00001) |
1007 | if (mFoliageMass<0.00001) |
1006 | die(); |
1008 | die(); |
1007 | 1009 | ||
1008 | double p_intrinsic, p_stress=0.; |
1010 | double p_intrinsic, p_stress=0.; |
1009 | p_intrinsic = species()->deathProb_intrinsic(); |
1011 | p_intrinsic = species()->deathProb_intrinsic(); |
1010 | 1012 | ||
1011 | if (mDbhDelta < _stress_threshold) { |
1013 | if (mDbhDelta < _stress_threshold) { |
1012 | mStressIndex++; |
1014 | mStressIndex++; |
1013 | if (mStressIndex> _stress_years) |
1015 | if (mStressIndex> _stress_years) |
1014 | p_stress = _stress_death_prob; |
1016 | p_stress = _stress_death_prob; |
1015 | } else |
1017 | } else |
1016 | mStressIndex = 0; |
1018 | mStressIndex = 0; |
1017 | 1019 | ||
1018 | double p = drandom(); //0..1 |
1020 | double p = drandom(); //0..1 |
1019 | if (p<p_intrinsic + p_stress) { |
1021 | if (p<p_intrinsic + p_stress) { |
1020 | // die...
|
1022 | // die...
|
1021 | die(); |
1023 | die(); |
1022 | }
|
1024 | }
|
1023 | }
|
1025 | }
|
1024 | #endif
|
1026 | #endif
|
1025 | 1027 | ||
1026 | void Tree::notifyTreeRemoved(TreeRemovalType reason) |
1028 | void Tree::notifyTreeRemoved(TreeRemovalType reason) |
1027 | {
|
1029 | {
|
1028 | // add the volume of the current tree to the height grid
|
1030 | // add the volume of the current tree to the height grid
|
1029 | // this information is used to track the removed volume for stands based on grids (and for salvaging operations)
|
1031 | // this information is used to track the removed volume for stands based on grids (and for salvaging operations)
|
1030 | ABE::ForestManagementEngine *abe = GlobalSettings::instance()->model()->ABEngine(); |
1032 | ABE::ForestManagementEngine *abe = GlobalSettings::instance()->model()->ABEngine(); |
1031 | if (abe) |
1033 | if (abe) |
1032 | abe->notifyTreeRemoval(this, static_cast<int>(reason)); |
1034 | abe->notifyTreeRemoval(this, static_cast<int>(reason)); |
1033 | 1035 | ||
1034 | // tell disturbance modules that a tree died
|
1036 | // tell disturbance modules that a tree died
|
1035 | GlobalSettings::instance()->model()->modules()->treeDeath(this, static_cast<int>(reason) ); |
1037 | GlobalSettings::instance()->model()->modules()->treeDeath(this, static_cast<int>(reason) ); |
1036 | 1038 | ||
1037 | // create output for tree removals
|
1039 | // create output for tree removals
|
1038 | if (mRemovalOutput && mRemovalOutput->isEnabled()) |
1040 | if (mRemovalOutput && mRemovalOutput->isEnabled()) |
1039 | mRemovalOutput->execRemovedTree(this, static_cast<int>(reason)); |
1041 | mRemovalOutput->execRemovedTree(this, static_cast<int>(reason)); |
- | 1042 | if (mLSRemovalOutput && mLSRemovalOutput->isEnabled()) |
|
- | 1043 | mLSRemovalOutput->execRemovedTree(this, static_cast<int>(reason)); |
|
1040 | }
|
1044 | }
|
1041 | 1045 | ||
1042 | //////////////////////////////////////////////////
|
1046 | //////////////////////////////////////////////////
|
1043 | //// value functions
|
1047 | //// value functions
|
1044 | //////////////////////////////////////////////////
|
1048 | //////////////////////////////////////////////////
|
1045 | 1049 | ||
1046 | double Tree::volume() const |
1050 | double Tree::volume() const |
1047 | {
|
1051 | {
|
1048 | /// @see Species::volumeFactor() for details
|
1052 | /// @see Species::volumeFactor() for details
|
1049 | const double volume_factor = species()->volumeFactor(); |
1053 | const double volume_factor = species()->volumeFactor(); |
1050 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
1054 | const double volume = volume_factor * mDbh*mDbh*mHeight * 0.0001; // dbh in cm: cm/100 * cm/100 = cm*cm * 0.0001 = m2 |
1051 | return volume; |
1055 | return volume; |
1052 | }
|
1056 | }
|
1053 | 1057 | ||
1054 | /// return the basal area in m2
|
1058 | /// return the basal area in m2
|
1055 | double Tree::basalArea() const |
1059 | double Tree::basalArea() const |
1056 | {
|
1060 | {
|
1057 | // A = r^2 * pi = d/2*pi; from cm->m: d/200
|
1061 | // A = r^2 * pi = d/2*pi; from cm->m: d/200
|
1058 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
1062 | const double b = (mDbh/200.)*(mDbh/200.)*M_PI; |
1059 | return b; |
1063 | return b; |
1060 | }
|
1064 | }
|
1061 | 1065 | ||
1062 | 1066 |