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