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