Rev 1204 | Rev 1217 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed
Rev | Author | Line No. | Line |
---|---|---|---|
1 | |||
671 | werner | 2 | /******************************************************************************************** |
3 | ** iLand - an individual based forest landscape and disturbance model |
||
4 | ** http://iland.boku.ac.at |
||
5 | ** Copyright (C) 2009- Werner Rammer, Rupert Seidl |
||
6 | ** |
||
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 |
||
9 | ** the Free Software Foundation, either version 3 of the License, or |
||
10 | ** (at your option) any later version. |
||
11 | ** |
||
12 | ** This program is distributed in the hope that it will be useful, |
||
13 | ** but WITHOUT ANY WARRANTY; without even the implied warranty of |
||
14 | ** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
||
15 | ** GNU General Public License for more details. |
||
16 | ** |
||
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/>. |
||
19 | ********************************************************************************************/ |
||
20 | |||
92 | Werner | 21 | /** @class Model |
247 | werner | 22 | Main object of the iLand model composited of various sub models / sub components. |
697 | werner | 23 | @ingroup core |
24 | The class Model is the top level container of iLand. The Model holds a collection of ResourceUnits, links to SpeciesSet and Climate. |
||
25 | ResourceUnit are grid cells with (currently) a size of 1 ha (100x100m). Many stand level processes (NPP produciton, WaterCycle) operate on this |
||
26 | level. |
||
27 | The Model also contain the landscape-wide 2m LIF-grid (http://iland.boku.ac.at/competition+for+light). |
||
28 | |||
92 | Werner | 29 | */ |
30 | #include "global.h" |
||
31 | #include "model.h" |
||
286 | werner | 32 | #include "sqlhelper.h" |
1182 | werner | 33 | #include "version.h" |
92 | Werner | 34 | |
105 | Werner | 35 | #include "xmlhelper.h" |
808 | werner | 36 | #include "debugtimer.h" |
281 | werner | 37 | #include "environment.h" |
340 | werner | 38 | #include "timeevents.h" |
106 | Werner | 39 | #include "helper.h" |
189 | iland | 40 | #include "resourceunit.h" |
208 | werner | 41 | #include "climate.h" |
92 | Werner | 42 | #include "speciesset.h" |
106 | Werner | 43 | #include "standloader.h" |
44 | #include "tree.h" |
||
185 | werner | 45 | #include "management.h" |
1111 | werner | 46 | #include "saplings.h" |
200 | werner | 47 | #include "modelsettings.h" |
615 | werner | 48 | #include "standstatistics.h" |
543 | werner | 49 | #include "mapgrid.h" |
632 | werner | 50 | #include "modelcontroller.h" |
641 | werner | 51 | #include "modules.h" |
654 | werner | 52 | #include "dem.h" |
1062 | werner | 53 | #include "grasscover.h" |
92 | Werner | 54 | |
202 | werner | 55 | #include "outputmanager.h" |
176 | werner | 56 | |
890 | werner | 57 | #include "forestmanagementengine.h" |
58 | |||
105 | Werner | 59 | #include <QtCore> |
60 | #include <QtXml> |
||
61 | |||
107 | Werner | 62 | /** iterate over all trees of the model. return NULL if all trees processed. |
63 | Usage: |
||
64 | @code |
||
65 | AllTreeIterator trees(model); |
||
66 | while (Tree *tree = trees.next()) { // returns NULL when finished. |
||
67 | tree->something(); // do something |
||
68 | } |
||
281 | werner | 69 | @endcode */ |
107 | Werner | 70 | Tree *AllTreeIterator::next() |
71 | { |
||
143 | Werner | 72 | |
107 | Werner | 73 | if (!mTreeEnd) { |
74 | // initialize to first ressource unit |
||
75 | mRUIterator = mModel->ruList().constBegin(); |
||
314 | werner | 76 | // fast forward to the first RU with trees |
77 | while (mRUIterator!=mModel->ruList().constEnd()) { |
||
78 | if ((*mRUIterator)->trees().count()>0) |
||
79 | break; |
||
753 | werner | 80 | ++mRUIterator; |
314 | werner | 81 | } |
82 | // finished if all RU processed |
||
83 | if (mRUIterator == mModel->ruList().constEnd()) |
||
84 | return NULL; |
||
143 | Werner | 85 | mTreeEnd = &((*mRUIterator)->trees().back()) + 1; // let end point to "1 after end" (STL-style) |
107 | Werner | 86 | mCurrent = &((*mRUIterator)->trees().front()); |
87 | } |
||
88 | if (mCurrent==mTreeEnd) { |
||
753 | werner | 89 | ++mRUIterator; // switch to next RU (loop until RU with trees is found) |
314 | werner | 90 | while (mRUIterator!=mModel->ruList().constEnd()) { |
91 | if ((*mRUIterator)->trees().count()>0) { |
||
92 | break; |
||
93 | } |
||
753 | werner | 94 | ++mRUIterator; |
314 | werner | 95 | } |
107 | Werner | 96 | if (mRUIterator == mModel->ruList().constEnd()) { |
143 | Werner | 97 | mCurrent = NULL; |
98 | return NULL; // finished!! |
||
107 | Werner | 99 | }else { |
143 | Werner | 100 | mTreeEnd = &((*mRUIterator)->trees().back()) + 1; |
107 | Werner | 101 | mCurrent = &((*mRUIterator)->trees().front()); |
102 | } |
||
103 | } |
||
143 | Werner | 104 | |
105 | return mCurrent++; |
||
107 | Werner | 106 | } |
157 | werner | 107 | Tree *AllTreeIterator::nextLiving() |
108 | { |
||
109 | while (Tree *t = next()) |
||
158 | werner | 110 | if (!t->isDead()) return t; |
157 | werner | 111 | return NULL; |
112 | } |
||
113 | Tree *AllTreeIterator::current() const |
||
114 | { |
||
115 | return mCurrent?mCurrent-1:NULL; |
||
116 | } |
||
107 | Werner | 117 | |
118 | |||
200 | werner | 119 | ModelSettings Model::mSettings; |
92 | Werner | 120 | Model::Model() |
121 | { |
||
122 | initialize(); |
||
137 | Werner | 123 | GlobalSettings::instance()->setModel(this); |
767 | werner | 124 | GlobalSettings::instance()->resetScriptEngine(); // clear the script |
1204 | werner | 125 | QString dbg="extended debug checks disabled."; |
126 | DBGMODE( dbg="extended debug checks enabled."; ); |
||
130 | Werner | 127 | qDebug() << dbg; |
92 | Werner | 128 | } |
129 | |||
130 | Model::~Model() |
||
131 | { |
||
132 | clear(); |
||
137 | Werner | 133 | GlobalSettings::instance()->setModel(NULL); |
92 | Werner | 134 | } |
135 | |||
136 | /** Initial setup of the Model. |
||
137 | */ |
||
138 | void Model::initialize() |
||
139 | { |
||
151 | iland | 140 | mSetup = false; |
162 | werner | 141 | GlobalSettings::instance()->setCurrentYear(0); |
151 | iland | 142 | mGrid = 0; |
143 | mHeightGrid = 0; |
||
185 | werner | 144 | mManagement = 0; |
909 | werner | 145 | mABEManagement = 0; |
281 | werner | 146 | mEnvironment = 0; |
340 | werner | 147 | mTimeEvents = 0; |
549 | werner | 148 | mStandGrid = 0; |
641 | werner | 149 | mModules = 0; |
654 | werner | 150 | mDEM = 0; |
1062 | werner | 151 | mGrassCover = 0; |
1111 | werner | 152 | mSaplings=0; |
92 | Werner | 153 | } |
154 | |||
261 | werner | 155 | /** sets up the simulation space. |
156 | */ |
||
103 | Werner | 157 | void Model::setupSpace() |
158 | { |
||
194 | werner | 159 | XmlHelper xml(GlobalSettings::instance()->settings().node("model.world")); |
192 | werner | 160 | double cellSize = xml.value("cellSize", "2").toDouble(); |
161 | double width = xml.value("width", "100").toDouble(); |
||
162 | double height = xml.value("height", "100").toDouble(); |
||
549 | werner | 163 | double buffer = xml.value("buffer", "60").toDouble(); |
164 | mModelRect = QRectF(0., 0., width, height); |
||
165 | |||
103 | Werner | 166 | qDebug() << QString("setup of the world: %1x%2m with cell-size=%3m and %4m buffer").arg(width).arg(height).arg(cellSize).arg(buffer); |
167 | |||
168 | QRectF total_grid(QPointF(-buffer, -buffer), QPointF(width+buffer, height+buffer)); |
||
169 | qDebug() << "setup grid rectangle:" << total_grid; |
||
170 | |||
151 | iland | 171 | if (mGrid) |
172 | delete mGrid; |
||
103 | Werner | 173 | mGrid = new FloatGrid(total_grid, cellSize); |
156 | werner | 174 | mGrid->initialize(1.f); |
151 | iland | 175 | if (mHeightGrid) |
176 | delete mHeightGrid; |
||
1165 | werner | 177 | mHeightGrid = new HeightGrid(total_grid, cellSize*cPxPerHeight); |
285 | werner | 178 | mHeightGrid->wipe(); // set all to zero |
156 | werner | 179 | Tree::setGrid(mGrid, mHeightGrid); |
105 | Werner | 180 | |
569 | werner | 181 | // setup the spatial location of the project area |
182 | if (xml.hasNode("location")) { |
||
183 | // setup of spatial location |
||
184 | double loc_x = xml.valueDouble("location.x"); |
||
185 | double loc_y = xml.valueDouble("location.y"); |
||
186 | double loc_z = xml.valueDouble("location.z"); |
||
187 | double loc_rot = xml.valueDouble("location.rotation"); |
||
188 | setupGISTransformation(loc_x, loc_y, loc_z, loc_rot); |
||
189 | qDebug() << "setup of spatial location: x/y/z" << loc_x << loc_y << loc_z << "rotation:" << loc_rot; |
||
190 | } else { |
||
191 | setupGISTransformation(0., 0., 0., 0.); |
||
192 | } |
||
193 | |||
567 | werner | 194 | // load environment (multiple climates, speciesSets, ... |
195 | if (mEnvironment) |
||
196 | delete mEnvironment; |
||
197 | mEnvironment = new Environment(); |
||
281 | werner | 198 | |
567 | werner | 199 | if (xml.valueBool("environmentEnabled", false)) { |
200 | QString env_file = GlobalSettings::instance()->path(xml.value("environmentFile")); |
||
201 | bool grid_mode = (xml.value("environmentMode")=="grid"); |
||
202 | QString grid_file = GlobalSettings::instance()->path(xml.value("environmentGrid")); |
||
893 | werner | 203 | if (grid_mode) { |
1210 | werner | 204 | if (QFile::exists(grid_file) && !xml.value("environmentGrid").isEmpty()) |
893 | werner | 205 | mEnvironment->setGridMode(grid_file); |
206 | else |
||
207 | throw IException(QString("File '%1' specified in key 'environmentGrid' does not exit ('environmentMode' is 'grid').").arg(grid_file) ); |
||
208 | } |
||
567 | werner | 209 | |
210 | if (!mEnvironment->loadFromFile(env_file)) |
||
211 | return; |
||
212 | } else { |
||
213 | // load and prepare default values |
||
214 | // (2) SpeciesSets: currently only one a global species set. |
||
215 | SpeciesSet *speciesSet = new SpeciesSet(); |
||
216 | mSpeciesSets.push_back(speciesSet); |
||
217 | speciesSet->setup(); |
||
218 | // Climate... |
||
219 | Climate *c = new Climate(); |
||
220 | mClimates.push_back(c); |
||
221 | mEnvironment->setDefaultValues(c, speciesSet); |
||
222 | } // environment? |
||
223 | |||
224 | // time series data |
||
225 | if (xml.valueBool(".timeEventsEnabled", false)) { |
||
226 | mTimeEvents = new TimeEvents(); |
||
584 | werner | 227 | mTimeEvents->loadFromFile(GlobalSettings::instance()->path(xml.value("timeEventsFile"), "script")); |
567 | werner | 228 | } |
229 | |||
646 | werner | 230 | |
105 | Werner | 231 | // simple case: create ressource units in a regular grid. |
194 | werner | 232 | if (xml.valueBool("resourceUnitsAsGrid")) { |
881 | werner | 233 | |
281 | werner | 234 | mRUmap.setup(QRectF(0., 0., width, height),100.); // Grid, that holds positions of resource units |
881 | werner | 235 | mRUmap.wipe(); |
646 | werner | 236 | |
539 | werner | 237 | bool mask_is_setup = false; |
238 | if (xml.valueBool("standGrid.enabled")) { |
||
239 | QString fileName = GlobalSettings::instance()->path(xml.value("standGrid.fileName")); |
||
882 | werner | 240 | mStandGrid = new MapGrid(fileName,false); // create stand grid index later |
543 | werner | 241 | |
549 | werner | 242 | if (mStandGrid->isValid()) { |
714 | werner | 243 | for (int i=0;i<mStandGrid->grid().count();i++) { |
244 | const int &grid_value = mStandGrid->grid().constValueAtIndex(i); |
||
245 | mHeightGrid->valueAtIndex(i).setValid( grid_value > -1 ); |
||
881 | werner | 246 | if (grid_value>-1) |
247 | mRUmap.valueAt(mStandGrid->grid().cellCenterPoint(i)) = (ResourceUnit*)1; |
||
714 | werner | 248 | if (grid_value < -1) |
718 | werner | 249 | mHeightGrid->valueAtIndex(i).setForestOutside(true); |
714 | werner | 250 | } |
539 | werner | 251 | } |
252 | mask_is_setup = true; |
||
802 | werner | 253 | } else { |
1182 | werner | 254 | if (!settings().torusMode) { |
802 | werner | 255 | // in the case we have no stand grid but only a large rectangle (without the torus option) |
256 | // we assume a forest outside |
||
257 | for (int i=0;i<mHeightGrid->count();++i) { |
||
258 | const QPointF &p = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i)); |
||
259 | if (p.x() < 0. || p.x()>width || p.y()<0. || p.y()>height) { |
||
260 | mHeightGrid->valueAtIndex(i).setForestOutside(true); |
||
261 | mHeightGrid->valueAtIndex(i).setValid(false); |
||
262 | } |
||
263 | } |
||
264 | |||
265 | } |
||
539 | werner | 266 | } |
267 | |||
881 | werner | 268 | ResourceUnit **p; // ptr to ptr! |
269 | ResourceUnit *new_ru; |
||
270 | |||
271 | int ru_index = 0; |
||
272 | for (p=mRUmap.begin(); p!=mRUmap.end(); ++p) { |
||
273 | QRectF r = mRUmap.cellRect(mRUmap.indexOf(p)); |
||
1032 | werner | 274 | if (!mStandGrid || !mStandGrid->isValid() || *p>NULL) { |
917 | werner | 275 | mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used. |
881 | werner | 276 | // create resource units for valid positions only |
277 | new_ru = new ResourceUnit(ru_index++); // create resource unit |
||
278 | new_ru->setClimate( mEnvironment->climate() ); |
||
279 | new_ru->setSpeciesSet( mEnvironment->speciesSet() ); |
||
280 | new_ru->setup(); |
||
281 | new_ru->setID( mEnvironment->currentID() ); // set id of resource unit in grid mode |
||
282 | new_ru->setBoundingBox(r); |
||
283 | mRU.append(new_ru); |
||
889 | werner | 284 | *p = new_ru; // save in the RUmap grid |
881 | werner | 285 | } |
286 | } |
||
1013 | werner | 287 | if (mEnvironment) { |
288 | // retrieve species sets and climates (that were really used) |
||
289 | mSpeciesSets << mEnvironment->speciesSetList(); |
||
290 | mClimates << mEnvironment->climateList(); |
||
1064 | werner | 291 | QString climate_file_list; |
292 | for (int i=0, c=0;i<mClimates.count();++i) { |
||
293 | climate_file_list += mClimates[i]->name() + ", "; |
||
294 | if (++c>5) { |
||
295 | climate_file_list += "..."; |
||
296 | break; |
||
297 | } |
||
881 | werner | 298 | |
1064 | werner | 299 | } |
300 | qDebug() << "Setup of climates: #loaded:" << mClimates.count() << "tables:" << climate_file_list; |
||
301 | |||
302 | |||
1013 | werner | 303 | } |
304 | |||
1011 | werner | 305 | qDebug() << "setup of" << mEnvironment->climateList().size() << "climates performed."; |
306 | |||
889 | werner | 307 | if (mStandGrid && mStandGrid->isValid()) |
882 | werner | 308 | mStandGrid->createIndex(); |
881 | werner | 309 | // now store the pointers in the grid. |
310 | // Important: This has to be done after the mRU-QList is complete - otherwise pointers would |
||
311 | // point to invalid memory when QList's memory is reorganized (expanding) |
||
312 | // ru_index = 0; |
||
313 | // for (p=mRUmap.begin();p!=mRUmap.end(); ++p) { |
||
314 | // *p = mRU.value(ru_index++); |
||
315 | // } |
||
316 | qDebug() << "created a grid of ResourceUnits: count=" << mRU.count() << "number of RU-map-cells:" << mRUmap.count(); |
||
317 | |||
318 | |||
574 | werner | 319 | calculateStockableArea(); |
320 | |||
285 | werner | 321 | // setup of the project area mask |
539 | werner | 322 | if (!mask_is_setup && xml.valueBool("areaMask.enabled", false) && xml.hasNode("areaMask.imageFile")) { |
285 | werner | 323 | // to be extended!!! e.g. to load ESRI-style text files.... |
324 | // setup a grid with the same size as the height grid... |
||
325 | FloatGrid tempgrid((int)mHeightGrid->cellsize(), mHeightGrid->sizeX(), mHeightGrid->sizeY()); |
||
326 | QString fileName = GlobalSettings::instance()->path(xml.value("areaMask.imageFile")); |
||
327 | loadGridFromImage(fileName, tempgrid); // fetch from image |
||
328 | for (int i=0;i<tempgrid.count(); i++) |
||
329 | mHeightGrid->valueAtIndex(i).setValid( tempgrid.valueAtIndex(i)>0.99 ); |
||
330 | qDebug() << "loaded project area mask from" << fileName; |
||
331 | } |
||
332 | |||
590 | werner | 333 | // list of "valid" resource units |
334 | QList<ResourceUnit*> valid_rus; |
||
335 | foreach(ResourceUnit* ru, mRU) |
||
336 | if (ru->id()!=-1) |
||
337 | valid_rus.append(ru); |
||
338 | |||
654 | werner | 339 | // setup of the digital elevation map (if present) |
340 | QString dem_file = xml.value("DEM"); |
||
341 | if (!dem_file.isEmpty()) { |
||
656 | werner | 342 | mDEM = new DEM(GlobalSettings::instance()->path(dem_file)); |
343 | // add them to the visuals... |
||
344 | GlobalSettings::instance()->controller()->addGrid(mDEM, "DEM height", GridViewRainbow, 0, 1000); |
||
345 | GlobalSettings::instance()->controller()->addGrid(mDEM->slopeGrid(), "DEM slope", GridViewRainbow, 0, 3); |
||
346 | GlobalSettings::instance()->controller()->addGrid(mDEM->aspectGrid(), "DEM aspect", GridViewRainbow, 0, 360); |
||
347 | GlobalSettings::instance()->controller()->addGrid(mDEM->viewGrid(), "DEM view", GridViewGray, 0, 1); |
||
348 | |||
654 | werner | 349 | } |
350 | |||
1111 | werner | 351 | // setup of saplings |
352 | if (mSaplings) { |
||
353 | delete mSaplings; mSaplings=0; |
||
354 | } |
||
355 | if (settings().regenerationEnabled) { |
||
356 | mSaplings = new Saplings(); |
||
357 | mSaplings->setup(); |
||
358 | } |
||
359 | |||
360 | |||
1062 | werner | 361 | // setup of the grass cover |
362 | if (!mGrassCover) |
||
363 | mGrassCover = new GrassCover(); |
||
364 | mGrassCover->setup(); |
||
365 | |||
646 | werner | 366 | // setup of external modules |
367 | mModules->setup(); |
||
368 | if (mModules->hasSetupResourceUnits()) { |
||
369 | for (ResourceUnit **p=mRUmap.begin(); p!=mRUmap.end(); ++p) { |
||
1087 | werner | 370 | if (*p) { |
371 | QRectF r = mRUmap.cellRect(mRUmap.indexOf(p)); |
||
372 | mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used. |
||
373 | mModules->setupResourceUnit( *p ); |
||
374 | } |
||
646 | werner | 375 | } |
376 | } |
||
377 | |||
767 | werner | 378 | // setup of scripting environment |
379 | ScriptGlobal::setupGlobalScripting(); |
||
380 | |||
123 | Werner | 381 | // setup the helper that does the multithreading |
590 | werner | 382 | threadRunner.setup(valid_rus); |
194 | werner | 383 | threadRunner.setMultithreading(GlobalSettings::instance()->settings().valueBool("system.settings.multithreading")); |
123 | Werner | 384 | threadRunner.print(); |
105 | Werner | 385 | |
1071 | werner | 386 | |
194 | werner | 387 | } else { |
388 | throw IException("resourceUnitsAsGrid MUST be set to true - at least currently :)"); |
||
105 | Werner | 389 | } |
120 | Werner | 390 | mSetup = true; |
103 | Werner | 391 | } |
392 | |||
393 | |||
92 | Werner | 394 | /** clear() frees all ressources allocated with the run of a simulation. |
395 | |||
396 | */ |
||
397 | void Model::clear() |
||
398 | { |
||
143 | Werner | 399 | mSetup = false; |
151 | iland | 400 | qDebug() << "Model clear: attempting to clear" << mRU.count() << "RU, " << mSpeciesSets.count() << "SpeciesSets."; |
92 | Werner | 401 | // clear ressource units |
103 | Werner | 402 | qDeleteAll(mRU); // delete ressource units (and trees) |
92 | Werner | 403 | mRU.clear(); |
103 | Werner | 404 | |
405 | qDeleteAll(mSpeciesSets); // delete species sets |
||
92 | Werner | 406 | mSpeciesSets.clear(); |
103 | Werner | 407 | |
208 | werner | 408 | // delete climate data |
409 | qDeleteAll(mClimates); |
||
410 | |||
151 | iland | 411 | // delete the grids |
412 | if (mGrid) |
||
413 | delete mGrid; |
||
414 | if (mHeightGrid) |
||
415 | delete mHeightGrid; |
||
1111 | werner | 416 | if (mSaplings) |
417 | delete mSaplings; |
||
185 | werner | 418 | if (mManagement) |
419 | delete mManagement; |
||
281 | werner | 420 | if (mEnvironment) |
421 | delete mEnvironment; |
||
340 | werner | 422 | if (mTimeEvents) |
423 | delete mTimeEvents; |
||
549 | werner | 424 | if (mStandGrid) |
425 | delete mStandGrid; |
||
641 | werner | 426 | if (mModules) |
427 | delete mModules; |
||
654 | werner | 428 | if (mDEM) |
429 | delete mDEM; |
||
1062 | werner | 430 | if (mGrassCover) |
431 | delete mGrassCover; |
||
909 | werner | 432 | if (mABEManagement) |
433 | delete mABEManagement; |
||
103 | Werner | 434 | |
185 | werner | 435 | mGrid = 0; |
436 | mHeightGrid = 0; |
||
437 | mManagement = 0; |
||
281 | werner | 438 | mEnvironment = 0; |
340 | werner | 439 | mTimeEvents = 0; |
549 | werner | 440 | mStandGrid = 0; |
641 | werner | 441 | mModules = 0; |
654 | werner | 442 | mDEM = 0; |
1062 | werner | 443 | mGrassCover = 0; |
909 | werner | 444 | mABEManagement = 0; |
185 | werner | 445 | |
583 | werner | 446 | GlobalSettings::instance()->outputManager()->close(); |
447 | |||
92 | Werner | 448 | qDebug() << "Model ressources freed."; |
449 | } |
||
450 | |||
451 | /** Setup of the Simulation. |
||
452 | This really creates the simulation environment and does the setup of various aspects. |
||
453 | */ |
||
102 | Werner | 454 | void Model::loadProject() |
92 | Werner | 455 | { |
109 | Werner | 456 | DebugTimer dt("load project"); |
93 | Werner | 457 | GlobalSettings *g = GlobalSettings::instance(); |
1157 | werner | 458 | g->printDirectories(); |
102 | Werner | 459 | const XmlHelper &xml = g->settings(); |
189 | iland | 460 | |
93 | Werner | 461 | g->clearDatabaseConnections(); |
92 | Werner | 462 | // database connections: reset |
94 | Werner | 463 | GlobalSettings::instance()->clearDatabaseConnections(); |
286 | werner | 464 | // input and climate connection |
465 | // see initOutputDatabase() for output database |
||
191 | werner | 466 | QString dbPath = g->path( xml.value("system.database.in"), "database"); |
194 | werner | 467 | GlobalSettings::instance()->setupDatabaseConnection("in", dbPath, true); |
193 | werner | 468 | dbPath = g->path( xml.value("system.database.climate"), "database"); |
194 | werner | 469 | GlobalSettings::instance()->setupDatabaseConnection("climate", dbPath, true); |
92 | Werner | 470 | |
200 | werner | 471 | mSettings.loadModelSettings(); |
472 | mSettings.print(); |
||
416 | werner | 473 | // random seed: if stored value is <> 0, use this as the random seed (and produce hence always an equal sequence of random numbers) |
474 | uint seed = xml.value("system.settings.randomSeed","0").toUInt(); |
||
708 | werner | 475 | RandomGenerator::setup(RandomGenerator::ergMersenneTwister, seed); // use the MersenneTwister as default |
428 | werner | 476 | // linearization of expressions: if true *and* linearize() is explicitely called, then |
477 | // function results will be cached over a defined range of values. |
||
478 | bool do_linearization = xml.valueBool("system.settings.expressionLinearizationEnabled", false); |
||
479 | Expression::setLinearizationEnabled(do_linearization); |
||
431 | werner | 480 | if (do_linearization) |
1157 | werner | 481 | qDebug() << "The linearization of expressions is enabled (performance optimization)."; |
431 | werner | 482 | |
483 | // log level |
||
484 | QString log_level = xml.value("system.settings.logLevel", "debug").toLower(); |
||
485 | if (log_level=="debug") setLogLevel(0); |
||
486 | if (log_level=="info") setLogLevel(1); |
||
487 | if (log_level=="warning") setLogLevel(2); |
||
488 | if (log_level=="error") setLogLevel(3); |
||
489 | |||
475 | werner | 490 | // snag dynamics / soil model enabled? (info used during setup of world) |
491 | changeSettings().carbonCycleEnabled = xml.valueBool("model.settings.carbonCycleEnabled", false); |
||
528 | werner | 492 | // class size of snag classes |
493 | Snag::setupThresholds(xml.valueDouble("model.settings.soil.swdDBHClass12"), |
||
494 | xml.valueDouble("model.settings.soil.swdDBHClass23")); |
||
475 | werner | 495 | |
646 | werner | 496 | // setup of modules |
497 | if (mModules) |
||
498 | delete mModules; |
||
499 | mModules = new Modules(); |
||
500 | |||
1111 | werner | 501 | changeSettings().regenerationEnabled = xml.valueBool("model.settings.regenerationEnabled", false); |
502 | |||
503 | |||
103 | Werner | 504 | setupSpace(); |
281 | werner | 505 | if (mRU.isEmpty()) |
506 | throw IException("Setup of Model: no resource units present!"); |
||
185 | werner | 507 | |
508 | // (3) additional issues |
||
837 | werner | 509 | // (3.1) load javascript code into the engine |
1065 | werner | 510 | QString script_file = xml.value("system.javascript.fileName"); |
1064 | werner | 511 | if (!script_file.isEmpty()) { |
1065 | werner | 512 | script_file = g->path(script_file, "script"); |
1064 | werner | 513 | ScriptGlobal::loadScript(script_file); |
514 | g->controller()->setLoadedJavascriptFile(script_file); |
||
515 | } |
||
837 | werner | 516 | |
517 | // (3.2) setup of regeneration |
||
391 | werner | 518 | if (settings().regenerationEnabled) { |
519 | foreach(SpeciesSet *ss, mSpeciesSets) |
||
520 | ss->setupRegeneration(); |
||
387 | werner | 521 | } |
1115 | werner | 522 | Saplings::setRecruitmentVariation(xml.valueDouble("model.settings.seedDispersal.recruitmentDimensionVariation",0.1)); |
468 | werner | 523 | |
837 | werner | 524 | // (3.3) management |
890 | werner | 525 | bool use_abe = xml.valueBool("model.management.abeEnabled"); |
526 | if (use_abe) { |
||
527 | // use the agent based forest management engine |
||
909 | werner | 528 | mABEManagement = new ABE::ForestManagementEngine(); |
890 | werner | 529 | // setup of ABE after loading of trees. |
530 | |||
185 | werner | 531 | } |
910 | werner | 532 | // use the standard management |
533 | QString mgmtFile = xml.value("model.management.file"); |
||
1194 | werner | 534 | if (xml.valueBool("model.management.enabled")) { |
910 | werner | 535 | mManagement = new Management(); |
536 | QString path = GlobalSettings::instance()->path(mgmtFile, "script"); |
||
537 | mManagement->loadScript(path); |
||
538 | qDebug() << "setup management using script" << path; |
||
539 | } |
||
468 | werner | 540 | |
541 | |||
910 | werner | 542 | |
92 | Werner | 543 | } |
105 | Werner | 544 | |
1058 | werner | 545 | void Model::reloadABE() |
546 | { |
||
547 | // delete firest |
||
548 | if (mABEManagement) |
||
549 | delete mABEManagement; |
||
550 | mABEManagement = new ABE::ForestManagementEngine(); |
||
551 | // and setup |
||
552 | mABEManagement->setup(); |
||
1089 | werner | 553 | mABEManagement->runOnInit(true); |
105 | Werner | 554 | |
1058 | werner | 555 | mABEManagement->initialize(); |
556 | |||
1089 | werner | 557 | mABEManagement->runOnInit(false); |
558 | |||
1058 | werner | 559 | } |
560 | |||
561 | |||
543 | werner | 562 | ResourceUnit *Model::ru(QPointF coord) |
105 | Werner | 563 | { |
106 | Werner | 564 | if (!mRUmap.isEmpty() && mRUmap.coordValid(coord)) |
565 | return mRUmap.valueAt(coord); |
||
1159 | werner | 566 | if (mRUmap.isEmpty()) |
567 | return ru(); // default RU if there is only one |
||
568 | else |
||
569 | return 0; // in this case, no valid coords were provided |
||
105 | Werner | 570 | } |
571 | |||
286 | werner | 572 | void Model::initOutputDatabase() |
573 | { |
||
574 | GlobalSettings *g = GlobalSettings::instance(); |
||
575 | QString dbPath = g->path(g->settings().value("system.database.out"), "output"); |
||
576 | // create run-metadata |
||
577 | int maxid = SqlHelper::queryValue("select max(id) from runs", g->dbin()).toInt(); |
||
578 | |||
579 | maxid++; |
||
580 | QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss"); |
||
287 | werner | 581 | SqlHelper::executeSql(QString("insert into runs (id, timestamp) values (%1, '%2')").arg(maxid).arg(timestamp), g->dbin()); |
286 | werner | 582 | // replace path information |
583 | dbPath.replace("$id$", QString::number(maxid)); |
||
584 | dbPath.replace("$date$", timestamp); |
||
585 | // setup final path |
||
586 | g->setupDatabaseConnection("out", dbPath, false); |
||
587 | |||
588 | } |
||
589 | |||
1115 | werner | 590 | /// multithreaded running function for the resource unit level establishment |
1157 | werner | 591 | void nc_establishment(ResourceUnit *unit) |
1115 | werner | 592 | { |
593 | Saplings *s = GlobalSettings::instance()->model()->saplings(); |
||
594 | try { |
||
595 | s->establishment(unit); |
||
903 | werner | 596 | |
1115 | werner | 597 | } catch (const IException& e) { |
598 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
599 | } |
||
600 | |||
601 | } |
||
602 | |||
824 | werner | 603 | /// multithreaded running function for the resource unit level establishment |
1157 | werner | 604 | void nc_sapling_growth(ResourceUnit *unit) |
440 | werner | 605 | { |
1115 | werner | 606 | Saplings *s = GlobalSettings::instance()->model()->saplings(); |
632 | werner | 607 | try { |
1115 | werner | 608 | s->saplingGrowth(unit); |
609 | |||
610 | } catch (const IException& e) { |
||
611 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
612 | } |
||
613 | } |
||
614 | |||
444 | werner | 615 | |
450 | werner | 616 | |
475 | werner | 617 | /// multithreaded execution of the carbon cycle routine |
1157 | werner | 618 | void nc_carbonCycle(ResourceUnit *unit) |
475 | werner | 619 | { |
611 | werner | 620 | try { |
632 | werner | 621 | // (1) do calculations on snag dynamics for the resource unit |
622 | unit->calculateCarbonCycle(); |
||
623 | // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N) |
||
624 | } catch (const IException& e) { |
||
625 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
611 | werner | 626 | } |
627 | |||
475 | werner | 628 | } |
629 | |||
440 | werner | 630 | /// beforeRun performs several steps before the models starts running. |
631 | /// inter alia: * setup of the stands |
||
632 | /// * setup of the climates |
||
105 | Werner | 633 | void Model::beforeRun() |
634 | { |
||
395 | werner | 635 | // setup outputs |
636 | // setup output database |
||
539 | werner | 637 | if (GlobalSettings::instance()->dbout().isOpen()) |
638 | GlobalSettings::instance()->dbout().close(); |
||
395 | werner | 639 | initOutputDatabase(); |
640 | GlobalSettings::instance()->outputManager()->setup(); |
||
641 | GlobalSettings::instance()->clearDebugLists(); |
||
642 | |||
105 | Werner | 643 | // initialize stands |
967 | werner | 644 | StandLoader loader(this); |
251 | werner | 645 | { |
106 | Werner | 646 | DebugTimer loadtrees("load trees"); |
647 | loader.processInit(); |
||
251 | werner | 648 | } |
934 | werner | 649 | // initalization of ABE |
650 | if (mABEManagement) { |
||
651 | mABEManagement->setup(); |
||
1089 | werner | 652 | mABEManagement->runOnInit(true); |
934 | werner | 653 | } |
106 | Werner | 654 | |
214 | werner | 655 | // load climate |
251 | werner | 656 | { |
1024 | werner | 657 | if (logLevelDebug()) qDebug() << "attempting to load climate..." ; |
658 | DebugTimer loadclim("load climate"); |
||
659 | foreach(Climate *c, mClimates) { |
||
660 | if (!c->isSetup()) |
||
661 | c->setup(); |
||
662 | } |
||
663 | // load the first year of the climate database |
||
664 | foreach(Climate *c, mClimates) |
||
665 | c->nextYear(); |
||
214 | werner | 666 | |
251 | werner | 667 | } |
668 | |||
934 | werner | 669 | |
251 | werner | 670 | { DebugTimer loadinit("load standstatistics"); |
734 | werner | 671 | if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ; |
135 | Werner | 672 | Tree::setGrid(mGrid, mHeightGrid); |
737 | werner | 673 | // debugCheckAllTrees(); // introduced for debugging session (2012-04-06) |
135 | Werner | 674 | applyPattern(); |
675 | readPattern(); |
||
967 | werner | 676 | loader.processAfterInit(); // e.g. initialization of saplings |
106 | Werner | 677 | |
376 | werner | 678 | // force the compilation of initial stand statistics |
241 | werner | 679 | createStandStatistics(); |
251 | werner | 680 | } |
286 | werner | 681 | |
934 | werner | 682 | // initalization of ABE (now all stands are properly set up) |
909 | werner | 683 | if (mABEManagement) { |
934 | werner | 684 | mABEManagement->initialize(); |
1089 | werner | 685 | mABEManagement->runOnInit(false); |
890 | werner | 686 | } |
687 | |||
688 | |||
264 | werner | 689 | // outputs to create with inital state (without any growth) are called here: |
936 | werner | 690 | GlobalSettings::instance()->setCurrentYear(0); // set clock to "0" (for outputs with initial state) |
691 | |||
257 | werner | 692 | GlobalSettings::instance()->outputManager()->execute("stand"); // year=0 |
837 | werner | 693 | GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0 |
504 | werner | 694 | GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0 |
1182 | werner | 695 | GlobalSettings::instance()->outputManager()->execute("saplingdetail"); // year=0 |
264 | werner | 696 | GlobalSettings::instance()->outputManager()->execute("tree"); // year=0 |
395 | werner | 697 | GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0 |
264 | werner | 698 | |
257 | werner | 699 | GlobalSettings::instance()->setCurrentYear(1); // set to first year |
1024 | werner | 700 | |
105 | Werner | 701 | } |
702 | |||
331 | werner | 703 | /** Main model runner. |
704 | The sequence of actions is as follows: |
||
705 | (1) Load the climate of the new year |
||
706 | (2) Reset statistics for resource unit as well as for dead/managed trees |
||
714 | werner | 707 | (3) Invoke Management. |
331 | werner | 708 | (4) *after* that, calculate Light patterns |
709 | (5) 3PG on stand level, tree growth. Clear stand-statistcs before they are filled by single-tree-growth. calculate water cycle (with LAIs before management) |
||
440 | werner | 710 | (6) execute Regeneration |
714 | werner | 711 | (7) invoke disturbance modules |
1202 | werner | 712 | (8) calculate carbon cycle |
713 | (9) calculate statistics for the year |
||
714 | (10) write database outputs |
||
331 | werner | 715 | */ |
105 | Werner | 716 | void Model::runYear() |
717 | { |
||
223 | werner | 718 | DebugTimer t("Model::runYear()"); |
615 | werner | 719 | GlobalSettings::instance()->systemStatistics()->reset(); |
707 | werner | 720 | RandomGenerator::checkGenerator(); // see if we need to generate new numbers... |
649 | werner | 721 | // initalization at start of year for external modules |
722 | mModules->yearBegin(); |
||
903 | werner | 723 | |
341 | werner | 724 | // execute scheduled events for the current year |
725 | if (mTimeEvents) |
||
726 | mTimeEvents->run(); |
||
727 | |||
1024 | werner | 728 | // load the next year of the climate database (except for the first year - the first climate year is loaded immediately |
729 | if (GlobalSettings::instance()->currentYear()>1) { |
||
730 | foreach(Climate *c, mClimates) |
||
731 | c->nextYear(); |
||
732 | } |
||
214 | werner | 733 | |
278 | werner | 734 | // reset statistics |
735 | foreach(ResourceUnit *ru, mRU) |
||
736 | ru->newYear(); |
||
737 | |||
391 | werner | 738 | foreach(SpeciesSet *set, mSpeciesSets) |
739 | set->newYear(); |
||
890 | werner | 740 | // management classic |
615 | werner | 741 | if (mManagement) { |
742 | DebugTimer t("management"); |
||
278 | werner | 743 | mManagement->run(); |
615 | werner | 744 | GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed(); |
745 | } |
||
909 | werner | 746 | // ... or ABE (the agent based variant) |
747 | if (mABEManagement) { |
||
748 | DebugTimer t("ABE:run"); |
||
749 | mABEManagement->run(); |
||
890 | werner | 750 | GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed(); |
751 | } |
||
278 | werner | 752 | |
937 | werner | 753 | // if trees are dead/removed because of management, the tree lists |
754 | // need to be cleaned (and the statistics need to be recreated) |
||
1157 | werner | 755 | cleanTreeLists(true); // recalculate statistics (LAIs per species needed later in production) |
937 | werner | 756 | |
214 | werner | 757 | // process a cycle of individual growth |
277 | werner | 758 | applyPattern(); // create Light Influence Patterns |
759 | readPattern(); // readout light state of individual trees |
||
760 | grow(); // let the trees grow (growth on stand-level, tree-level, mortality) |
||
1062 | werner | 761 | mGrassCover->execute(); // evaluate the grass / herb cover (and its effect on regeneration) |
106 | Werner | 762 | |
391 | werner | 763 | // regeneration |
764 | if (settings().regenerationEnabled) { |
||
440 | werner | 765 | // seed dispersal |
1117 | werner | 766 | DebugTimer tseed("Seed dispersal, establishment, sapling growth"); |
391 | werner | 767 | foreach(SpeciesSet *set, mSpeciesSets) |
475 | werner | 768 | set->regeneration(); // parallel execution for each species set |
440 | werner | 769 | |
615 | werner | 770 | GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed(); |
440 | werner | 771 | // establishment |
1115 | werner | 772 | Saplings::updateBrowsingPressure(); |
1063 | werner | 773 | |
824 | werner | 774 | |
1115 | werner | 775 | { DebugTimer t("establishment"); |
1158 | werner | 776 | executePerResourceUnit( nc_establishment, false /* true: force single threaded operation */); |
777 | GlobalSettings::instance()->systemStatistics()->tEstablishment+=t.elapsed(); |
||
1115 | werner | 778 | } |
779 | { DebugTimer t("sapling growth"); |
||
1158 | werner | 780 | executePerResourceUnit( nc_sapling_growth, false /* true: force single threaded operation */); |
781 | GlobalSettings::instance()->systemStatistics()->tSapling+=t.elapsed(); |
||
1115 | werner | 782 | } |
783 | |||
1176 | werner | 784 | // Establishment::debugInfo(); // debug test |
1068 | werner | 785 | |
615 | werner | 786 | } |
440 | werner | 787 | |
1202 | werner | 788 | // external modules/disturbances |
789 | mModules->run(); |
||
790 | // cleanup of tree lists if external modules removed trees. |
||
791 | cleanTreeLists(false); // do not recalculate statistics - this is done in ru->yearEnd() |
||
792 | |||
793 | |||
468 | werner | 794 | // calculate soil / snag dynamics |
475 | werner | 795 | if (settings().carbonCycleEnabled) { |
796 | DebugTimer ccycle("carbon cylce"); |
||
1063 | werner | 797 | executePerResourceUnit( nc_carbonCycle, false /* true: force single threaded operation */); |
615 | werner | 798 | GlobalSettings::instance()->systemStatistics()->tCarbonCycle+=ccycle.elapsed(); |
799 | |||
475 | werner | 800 | } |
649 | werner | 801 | |
802 | |||
615 | werner | 803 | DebugTimer toutput("outputs"); |
278 | werner | 804 | // calculate statistics |
805 | foreach(ResourceUnit *ru, mRU) |
||
806 | ru->yearEnd(); |
||
124 | Werner | 807 | |
277 | werner | 808 | // create outputs |
184 | werner | 809 | OutputManager *om = GlobalSettings::instance()->outputManager(); |
285 | werner | 810 | om->execute("tree"); // single tree output |
1102 | werner | 811 | om->execute("treeremoval"); // single removed tree output |
278 | werner | 812 | om->execute("stand"); //resource unit level x species |
837 | werner | 813 | om->execute("landscape"); //landscape x species |
1114 | werner | 814 | om->execute("landscape_removed"); //removed trees on landscape x species |
504 | werner | 815 | om->execute("sapling"); // sapling layer per RU x species |
1182 | werner | 816 | om->execute("saplingdetail"); // individual sapling cohorts (per RU) |
278 | werner | 817 | om->execute("production_month"); // 3pg responses growth per species x RU x month |
285 | werner | 818 | om->execute("dynamicstand"); // output with user-defined columns (based on species x RU) |
278 | werner | 819 | om->execute("standdead"); // resource unit level x species |
820 | om->execute("management"); // resource unit level x species |
||
587 | werner | 821 | om->execute("carbon"); // resource unit level, carbon pools above and belowground |
609 | werner | 822 | om->execute("carbonflow"); // resource unit level, GPP, NPP and total carbon flows (atmosphere, harvest, ...) |
1157 | werner | 823 | om->execute("water"); // resource unit/landscape level water output (ET, rad, snow cover, ...) |
185 | werner | 824 | |
615 | werner | 825 | GlobalSettings::instance()->systemStatistics()->tWriteOutput+=toutput.elapsed(); |
826 | GlobalSettings::instance()->systemStatistics()->tTotalYear+=t.elapsed(); |
||
827 | GlobalSettings::instance()->systemStatistics()->writeOutput(); |
||
828 | |||
1157 | werner | 829 | // global javascript event |
830 | GlobalSettings::instance()->executeJSFunction("onYearEnd"); |
||
831 | |||
162 | werner | 832 | GlobalSettings::instance()->setCurrentYear(GlobalSettings::instance()->currentYear()+1); |
1085 | werner | 833 | |
834 | // try to clean up a bit of memory (useful if many large JS objects (e.g., grids) are used) |
||
835 | GlobalSettings::instance()->scriptEngine()->collectGarbage(); |
||
105 | Werner | 836 | } |
837 | |||
278 | werner | 838 | |
839 | |||
105 | Werner | 840 | void Model::afterStop() |
841 | { |
||
842 | // do some cleanup |
||
843 | } |
||
106 | Werner | 844 | |
369 | werner | 845 | /// multithreaded running function for LIP printing |
1157 | werner | 846 | void nc_applyPattern(ResourceUnit *unit) |
106 | Werner | 847 | { |
848 | |||
849 | QVector<Tree>::iterator tit; |
||
118 | Werner | 850 | QVector<Tree>::iterator tend = unit->trees().end(); |
107 | Werner | 851 | |
632 | werner | 852 | try { |
107 | Werner | 853 | |
632 | werner | 854 | // light concurrence influence |
1182 | werner | 855 | if (!GlobalSettings::instance()->model()->settings().torusMode) { |
632 | werner | 856 | // height dominance grid |
857 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
||
858 | (*tit).heightGrid(); // just do it ;) |
||
155 | werner | 859 | |
632 | werner | 860 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
861 | (*tit).applyLIP(); // just do it ;) |
||
155 | werner | 862 | |
632 | werner | 863 | } else { |
864 | // height dominance grid |
||
865 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
||
866 | (*tit).heightGrid_torus(); // just do it ;) |
||
155 | werner | 867 | |
632 | werner | 868 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
869 | (*tit).applyLIP_torus(); // do it the wraparound way |
||
870 | } |
||
1157 | werner | 871 | |
632 | werner | 872 | } catch (const IException &e) { |
873 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
118 | Werner | 874 | } |
875 | } |
||
106 | Werner | 876 | |
369 | werner | 877 | /// multithreaded running function for LIP value extraction |
1157 | werner | 878 | void nc_readPattern(ResourceUnit *unit) |
118 | Werner | 879 | { |
880 | QVector<Tree>::iterator tit; |
||
881 | QVector<Tree>::iterator tend = unit->trees().end(); |
||
632 | werner | 882 | try { |
1182 | werner | 883 | if (!GlobalSettings::instance()->model()->settings().torusMode) { |
632 | werner | 884 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
885 | (*tit).readLIF(); // multipliactive approach |
||
886 | } else { |
||
887 | for (tit=unit->trees().begin(); tit!=tend; ++tit) |
||
888 | (*tit).readLIF_torus(); // do it the wraparound way |
||
889 | } |
||
890 | } catch (const IException &e) { |
||
891 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
118 | Werner | 892 | } |
106 | Werner | 893 | } |
894 | |||
369 | werner | 895 | /// multithreaded running function for the growth of individual trees |
1157 | werner | 896 | void nc_grow(ResourceUnit *unit) |
118 | Werner | 897 | { |
898 | QVector<Tree>::iterator tit; |
||
899 | QVector<Tree>::iterator tend = unit->trees().end(); |
||
632 | werner | 900 | try { |
901 | unit->beforeGrow(); // reset statistics |
||
902 | // calculate light responses |
||
903 | // responses are based on *modified* values for LightResourceIndex |
||
904 | for (tit=unit->trees().begin(); tit!=tend; ++tit) { |
||
905 | (*tit).calcLightResponse(); |
||
906 | } |
||
118 | Werner | 907 | |
632 | werner | 908 | unit->calculateInterceptedArea(); |
251 | werner | 909 | |
632 | werner | 910 | for (tit=unit->trees().begin(); tit!=tend; ++tit) { |
911 | (*tit).grow(); // actual growth of individual trees |
||
912 | } |
||
913 | } catch (const IException &e) { |
||
914 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
118 | Werner | 915 | } |
615 | werner | 916 | |
917 | GlobalSettings::instance()->systemStatistics()->treeCount+=unit->trees().count(); |
||
118 | Werner | 918 | } |
919 | |||
369 | werner | 920 | /// multithreaded running function for the resource level production |
1157 | werner | 921 | void nc_production(ResourceUnit *unit) |
369 | werner | 922 | { |
632 | werner | 923 | try { |
924 | unit->production(); |
||
925 | } catch (const IException &e) { |
||
926 | GlobalSettings::instance()->controller()->throwError(e.message()); |
||
927 | } |
||
369 | werner | 928 | } |
929 | |||
440 | werner | 930 | |
124 | Werner | 931 | void Model::test() |
932 | { |
||
933 | // Test-funktion: braucht 1/3 time von readGrid() |
||
934 | DebugTimer t("test"); |
||
935 | FloatGrid averaged = mGrid->averaged(10); |
||
936 | int count = 0; |
||
937 | float *end = averaged.end(); |
||
938 | for (float *p=averaged.begin(); p!=end; ++p) |
||
939 | if (*p > 0.9) |
||
940 | count++; |
||
941 | qDebug() << count << "LIF>0.9 of " << averaged.count(); |
||
942 | } |
||
943 | |||
734 | werner | 944 | void Model::debugCheckAllTrees() |
945 | { |
||
946 | AllTreeIterator at(this); |
||
947 | bool has_errors = false; double dummy=0.; |
||
948 | while (Tree *t = at.next()) { |
||
949 | // plausibility |
||
950 | if (t->dbh()<0 || t->dbh()>10000. || t->biomassFoliage()<0. || t->height()>1000. || t->height() < 0. |
||
951 | || t->biomassFoliage() <0.) |
||
952 | has_errors = true; |
||
953 | // check for objects.... |
||
954 | dummy = t->stamp()->offset() + t->ru()->ruSpecies()[1]->statistics().count(); |
||
955 | } |
||
956 | if (has_errors) |
||
957 | qDebug() << "model: debugCheckAllTrees found problems" << dummy; |
||
958 | } |
||
959 | |||
118 | Werner | 960 | void Model::applyPattern() |
961 | { |
||
962 | |||
963 | DebugTimer t("applyPattern()"); |
||
964 | // intialize grids... |
||
720 | werner | 965 | initializeGrid(); |
551 | werner | 966 | |
406 | werner | 967 | // initialize height grid with a value of 4m. This is the height of the regeneration layer |
551 | werner | 968 | for (HeightGridValue *h=mHeightGrid->begin();h!=mHeightGrid->end();++h) { |
969 | h->resetCount(); // set count = 0, but do not touch the flags |
||
970 | h->height = 4.f; |
||
971 | } |
||
118 | Werner | 972 | |
123 | Werner | 973 | threadRunner.run(nc_applyPattern); |
615 | werner | 974 | GlobalSettings::instance()->systemStatistics()->tApplyPattern+=t.elapsed(); |
118 | Werner | 975 | } |
976 | |||
106 | Werner | 977 | void Model::readPattern() |
978 | { |
||
979 | DebugTimer t("readPattern()"); |
||
123 | Werner | 980 | threadRunner.run(nc_readPattern); |
615 | werner | 981 | GlobalSettings::instance()->systemStatistics()->tReadPattern+=t.elapsed(); |
982 | |||
106 | Werner | 983 | } |
984 | |||
331 | werner | 985 | /** Main function for the growth of stands and trees. |
986 | This includes several steps. |
||
987 | (1) calculate the stocked area (i.e. count pixels in height grid) |
||
988 | (2) 3PG production (including response calculation, water cycle) |
||
989 | (3) single tree growth (including mortality) |
||
990 | (4) cleanup of tree lists (remove dead trees) |
||
991 | */ |
||
106 | Werner | 992 | void Model::grow() |
993 | { |
||
151 | iland | 994 | |
1196 | werner | 995 | |
369 | werner | 996 | { DebugTimer t("growRU()"); |
997 | calculateStockedArea(); |
||
113 | Werner | 998 | |
1161 | werner | 999 | // Production of biomass (stand level, 3PG) |
374 | werner | 1000 | threadRunner.run(nc_production); |
370 | werner | 1001 | } |
1002 | |||
1003 | DebugTimer t("growTrees()"); |
||
251 | werner | 1004 | threadRunner.run(nc_grow); // actual growth of individual trees |
159 | werner | 1005 | |
187 | iland | 1006 | foreach(ResourceUnit *ru, mRU) { |
159 | werner | 1007 | ru->cleanTreeList(); |
376 | werner | 1008 | ru->afterGrow(); |
168 | werner | 1009 | //qDebug() << (b-n) << "trees died (of" << b << ")."; |
159 | werner | 1010 | } |
615 | werner | 1011 | GlobalSettings::instance()->systemStatistics()->tTreeGrowth+=t.elapsed(); |
123 | Werner | 1012 | } |
151 | iland | 1013 | |
240 | werner | 1014 | /** calculate for each resource unit the fraction of area which is stocked. |
1015 | This is done by checking the pixels of the global height grid. |
||
151 | iland | 1016 | */ |
1017 | void Model::calculateStockedArea() |
||
1018 | { |
||
1019 | // iterate over the whole heightgrid and count pixels for each ressource unit |
||
1020 | HeightGridValue *end = mHeightGrid->end(); |
||
1021 | QPointF cp; |
||
187 | iland | 1022 | ResourceUnit *ru; |
151 | iland | 1023 | for (HeightGridValue *i=mHeightGrid->begin(); i!=end; ++i) { |
1024 | cp = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i)); |
||
1025 | if (mRUmap.coordValid(cp)) { |
||
1026 | ru = mRUmap.valueAt(cp); |
||
1027 | if (ru) { |
||
285 | werner | 1028 | ru->countStockedPixel( (*i).count()>0 ); |
151 | iland | 1029 | } |
1030 | } |
||
1031 | |||
1032 | } |
||
1033 | } |
||
240 | werner | 1034 | |
664 | werner | 1035 | /** calculate for each resource unit the stockable area. |
574 | werner | 1036 | "stockability" is determined by the isValid flag of resource units which in turn |
1037 | is derived from stand grid values. |
||
1038 | */ |
||
1039 | void Model::calculateStockableArea() |
||
1040 | { |
||
1041 | |||
1114 | werner | 1042 | mTotalStockableArea = 0.; |
574 | werner | 1043 | foreach(ResourceUnit *ru, mRU) { |
720 | werner | 1044 | // // |
1045 | // if (ru->id()==-1) { |
||
1046 | // ru->setStockableArea(0.); |
||
1047 | // continue; |
||
1048 | // } |
||
574 | werner | 1049 | GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox()); |
1050 | int valid=0, total=0; |
||
1051 | while (runner.next()) { |
||
1052 | if ( runner.current()->isValid() ) |
||
1053 | valid++; |
||
1054 | total++; |
||
1055 | } |
||
575 | werner | 1056 | if (total) { |
1114 | werner | 1057 | ru->setStockableArea( cHeightPixelArea * valid); // in m2 |
1157 | werner | 1058 | if (ru->snag()) |
1059 | ru->snag()->scaleInitialState(); |
||
1060 | mTotalStockableArea += cHeightPixelArea * valid / cRUArea; // in ha |
||
734 | werner | 1061 | if (valid==0 && ru->id()>-1) { |
1062 | // invalidate this resource unit |
||
1063 | ru->setID(-1); |
||
1064 | } |
||
575 | werner | 1065 | if (valid>0 && ru->id()==-1) { |
1066 | qDebug() << "Warning: a resource unit has id=-1 but stockable area (id was set to 0)!!! ru: " << ru->boundingBox() << "with index" << ru->index(); |
||
1067 | ru->setID(0); |
||
664 | werner | 1068 | // test-code |
1069 | //GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox()); |
||
1070 | //while (runner.next()) { |
||
1071 | // qDebug() << mHeightGrid->cellCenterPoint(mHeightGrid->indexOf( runner.current() )) << ": " << runner.current()->isValid(); |
||
1072 | //} |
||
585 | werner | 1073 | |
575 | werner | 1074 | } |
1075 | } else |
||
574 | werner | 1076 | throw IException("calculateStockableArea: resource unit without pixels!"); |
1077 | |||
1078 | } |
||
720 | werner | 1079 | // mark those pixels that are at the edge of a "forest-out-of-area" |
1080 | GridRunner<HeightGridValue> runner(*mHeightGrid, mHeightGrid->metricRect()); |
||
1081 | HeightGridValue* neighbors[8]; |
||
1082 | while (runner.next()) { |
||
1083 | if (runner.current()->isForestOutside()) { |
||
1084 | // if the current pixel is a "radiating" border pixel, |
||
1085 | // then check the neighbors and set a flag if the pixel is a neighbor of a in-project-area pixel. |
||
1086 | runner.neighbors8(neighbors); |
||
1087 | for (int i=0;i<8;++i) |
||
1088 | if (neighbors[i] && neighbors[i]->isValid()) |
||
1089 | runner.current()->setIsRadiating(); |
||
1090 | |||
1091 | } |
||
1092 | } |
||
1093 | |||
1114 | werner | 1094 | qDebug() << "Total stockable area of the landscape is" << mTotalStockableArea << "ha."; |
1095 | |||
574 | werner | 1096 | } |
1097 | |||
720 | werner | 1098 | void Model::initializeGrid() |
1099 | { |
||
1100 | // fill the whole grid with a value of "1." |
||
1101 | mGrid->initialize(1.f); |
||
574 | werner | 1102 | |
720 | werner | 1103 | // apply special values for grid cells border regions where out-of-area cells |
1104 | // radiate into the main LIF grid. |
||
1105 | QPoint p; |
||
1106 | int ix_min, ix_max, iy_min, iy_max, ix_center, iy_center; |
||
1107 | const int px_offset = cPxPerHeight / 2; // for 5 px per height grid cell, the offset is 2 |
||
1108 | const int max_radiate_distance = 7; |
||
1109 | const float step_width = 1.f / (float)max_radiate_distance; |
||
1110 | int c_rad = 0; |
||
1111 | for (HeightGridValue *hgv=mHeightGrid->begin(); hgv!=mHeightGrid->end(); ++hgv) { |
||
1112 | if (hgv->isRadiating()) { |
||
1113 | p=mHeightGrid->indexOf(hgv); |
||
1114 | ix_min = p.x() * cPxPerHeight - max_radiate_distance + px_offset; |
||
1115 | ix_max = ix_min + 2*max_radiate_distance + 1; |
||
1116 | ix_center = ix_min + max_radiate_distance; |
||
1117 | iy_min = p.y() * cPxPerHeight - max_radiate_distance + px_offset; |
||
1118 | iy_max = iy_min + 2*max_radiate_distance + 1; |
||
1119 | iy_center = iy_min + max_radiate_distance; |
||
1120 | for (int y=iy_min; y<=iy_max; ++y) { |
||
1121 | for (int x=ix_min; x<=ix_max; ++x) { |
||
802 | werner | 1122 | if (!mGrid->isIndexValid(x,y) || !(*mHeightGrid)(x/cPxPerHeight, y/cPxPerHeight).isValid()) |
720 | werner | 1123 | continue; |
1124 | float value = qMax(qAbs(x-ix_center), qAbs(y-iy_center)) * step_width; |
||
1125 | float &v = mGrid->valueAtIndex(x, y); |
||
1126 | if (value>=0.f && v>value) |
||
1127 | v = value; |
||
1128 | } |
||
1129 | } |
||
1130 | c_rad++; |
||
1131 | } |
||
1132 | } |
||
721 | werner | 1133 | if (logLevelDebug()) |
1134 | qDebug() << "initialize grid:" << c_rad << "radiating pixels..."; |
||
720 | werner | 1135 | |
1136 | } |
||
1137 | |||
1138 | |||
241 | werner | 1139 | /// Force the creation of stand statistics. |
1140 | /// - stocked area (for resourceunit-areas) |
||
1141 | /// - ru - statistics |
||
240 | werner | 1142 | void Model::createStandStatistics() |
1143 | { |
||
241 | werner | 1144 | calculateStockedArea(); |
482 | werner | 1145 | foreach(ResourceUnit *ru, mRU) { |
1146 | ru->addTreeAgingForAllTrees(); |
||
240 | werner | 1147 | ru->createStandStatistics(); |
482 | werner | 1148 | } |
240 | werner | 1149 | } |
584 | werner | 1150 | |
1157 | werner | 1151 | void Model::cleanTreeLists(bool recalculate_stats) |
936 | werner | 1152 | { |
937 | werner | 1153 | foreach(ResourceUnit *ru, GlobalSettings::instance()->model()->ruList()) { |
1154 | if (ru->hasDiedTrees()) { |
||
1155 | ru->cleanTreeList(); |
||
1157 | werner | 1156 | ru->recreateStandStatistics(recalculate_stats); |
937 | werner | 1157 | } |
1158 | } |
||
936 | werner | 1159 | } |
584 | werner | 1160 | |
936 | werner | 1161 |