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