Subversion Repositories public iLand

Rev

Rev 1111 | Rev 1115 | Go to most recent revision | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 1111 Rev 1114
1
Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/model.cpp':
1
Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/model.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
20
21
/** @class Model
21
/** @class Model
22
  Main object of the iLand model composited of various sub models / sub components.
22
  Main object of the iLand model composited of various sub models / sub components.
23
  @ingroup core
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.
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
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.
26
  level.
27
  The Model also contain the landscape-wide 2m LIF-grid (http://iland.boku.ac.at/competition+for+light).
27
  The Model also contain the landscape-wide 2m LIF-grid (http://iland.boku.ac.at/competition+for+light).
28

28

29
  */
29
  */
30
#include "global.h"
30
#include "global.h"
31
#include "model.h"
31
#include "model.h"
32
#include "sqlhelper.h"
32
#include "sqlhelper.h"
33
33
34
#include "xmlhelper.h"
34
#include "xmlhelper.h"
35
#include "debugtimer.h"
35
#include "debugtimer.h"
36
#include "environment.h"
36
#include "environment.h"
37
#include "timeevents.h"
37
#include "timeevents.h"
38
#include "helper.h"
38
#include "helper.h"
39
#include "resourceunit.h"
39
#include "resourceunit.h"
40
#include "climate.h"
40
#include "climate.h"
41
#include "speciesset.h"
41
#include "speciesset.h"
42
#include "standloader.h"
42
#include "standloader.h"
43
#include "tree.h"
43
#include "tree.h"
44
#include "management.h"
44
#include "management.h"
45
#include "saplings.h"
45
#include "saplings.h"
46
#include "modelsettings.h"
46
#include "modelsettings.h"
47
#include "standstatistics.h"
47
#include "standstatistics.h"
48
#include "mapgrid.h"
48
#include "mapgrid.h"
49
#include "modelcontroller.h"
49
#include "modelcontroller.h"
50
#include "modules.h"
50
#include "modules.h"
51
#include "dem.h"
51
#include "dem.h"
52
#include "grasscover.h"
52
#include "grasscover.h"
53
53
54
#include "outputmanager.h"
54
#include "outputmanager.h"
55
55
56
#include "forestmanagementengine.h"
56
#include "forestmanagementengine.h"
57
57
58
#include <QtCore>
58
#include <QtCore>
59
#include <QtXml>
59
#include <QtXml>
60
60
61
/** iterate over all trees of the model. return NULL if all trees processed.
61
/** iterate over all trees of the model. return NULL if all trees processed.
62
  Usage:
62
  Usage:
63
  @code
63
  @code
64
  AllTreeIterator trees(model);
64
  AllTreeIterator trees(model);
65
  while (Tree *tree = trees.next()) { // returns NULL when finished.
65
  while (Tree *tree = trees.next()) { // returns NULL when finished.
66
     tree->something(); // do something
66
     tree->something(); // do something
67
  }
67
  }
68
  @endcode  */
68
  @endcode  */
69
Tree *AllTreeIterator::next()
69
Tree *AllTreeIterator::next()
70
{
70
{
71
71
72
    if (!mTreeEnd) {
72
    if (!mTreeEnd) {
73
        // initialize to first ressource unit
73
        // initialize to first ressource unit
74
        mRUIterator = mModel->ruList().constBegin();
74
        mRUIterator = mModel->ruList().constBegin();
75
        // fast forward to the first RU with trees
75
        // fast forward to the first RU with trees
76
        while (mRUIterator!=mModel->ruList().constEnd()) {
76
        while (mRUIterator!=mModel->ruList().constEnd()) {
77
            if ((*mRUIterator)->trees().count()>0)
77
            if ((*mRUIterator)->trees().count()>0)
78
                break;
78
                break;
79
            ++mRUIterator;
79
            ++mRUIterator;
80
        }
80
        }
81
            // finished if all RU processed
81
            // finished if all RU processed
82
        if (mRUIterator == mModel->ruList().constEnd())
82
        if (mRUIterator == mModel->ruList().constEnd())
83
            return NULL;
83
            return NULL;
84
        mTreeEnd = &((*mRUIterator)->trees().back()) + 1; // let end point to "1 after end" (STL-style)
84
        mTreeEnd = &((*mRUIterator)->trees().back()) + 1; // let end point to "1 after end" (STL-style)
85
        mCurrent = &((*mRUIterator)->trees().front());
85
        mCurrent = &((*mRUIterator)->trees().front());
86
    }
86
    }
87
    if (mCurrent==mTreeEnd) {
87
    if (mCurrent==mTreeEnd) {
88
        ++mRUIterator; // switch to next RU (loop until RU with trees is found)
88
        ++mRUIterator; // switch to next RU (loop until RU with trees is found)
89
        while (mRUIterator!=mModel->ruList().constEnd()) {
89
        while (mRUIterator!=mModel->ruList().constEnd()) {
90
            if ((*mRUIterator)->trees().count()>0) {
90
            if ((*mRUIterator)->trees().count()>0) {
91
                break;
91
                break;
92
            }
92
            }
93
            ++mRUIterator;
93
            ++mRUIterator;
94
        }
94
        }
95
        if (mRUIterator == mModel->ruList().constEnd()) {
95
        if (mRUIterator == mModel->ruList().constEnd()) {
96
            mCurrent = NULL;
96
            mCurrent = NULL;
97
            return NULL; // finished!!
97
            return NULL; // finished!!
98
        }else {
98
        }else {
99
            mTreeEnd = &((*mRUIterator)->trees().back()) + 1;
99
            mTreeEnd = &((*mRUIterator)->trees().back()) + 1;
100
            mCurrent = &((*mRUIterator)->trees().front());
100
            mCurrent = &((*mRUIterator)->trees().front());
101
        }
101
        }
102
    }
102
    }
103
103
104
    return mCurrent++;
104
    return mCurrent++;
105
}
105
}
106
Tree *AllTreeIterator::nextLiving()
106
Tree *AllTreeIterator::nextLiving()
107
{
107
{
108
    while (Tree *t = next())
108
    while (Tree *t = next())
109
        if (!t->isDead()) return t;
109
        if (!t->isDead()) return t;
110
    return NULL;
110
    return NULL;
111
}
111
}
112
Tree *AllTreeIterator::current() const
112
Tree *AllTreeIterator::current() const
113
{
113
{
114
    return mCurrent?mCurrent-1:NULL;
114
    return mCurrent?mCurrent-1:NULL;
115
}
115
}
116
116
117
117
118
ModelSettings Model::mSettings;
118
ModelSettings Model::mSettings;
119
Model::Model()
119
Model::Model()
120
{
120
{
121
    initialize();
121
    initialize();
122
    GlobalSettings::instance()->setModel(this);
122
    GlobalSettings::instance()->setModel(this);
123
    GlobalSettings::instance()->resetScriptEngine(); // clear the script
123
    GlobalSettings::instance()->resetScriptEngine(); // clear the script
124
    QString dbg="running in release mode.";
124
    QString dbg="running in release mode.";
125
    DBGMODE( dbg="running in debug mode."; );
125
    DBGMODE( dbg="running in debug mode."; );
126
    qDebug() << dbg;
126
    qDebug() << dbg;
127
}
127
}
128
128
129
Model::~Model()
129
Model::~Model()
130
{
130
{
131
    clear();
131
    clear();
132
    GlobalSettings::instance()->setModel(NULL);
132
    GlobalSettings::instance()->setModel(NULL);
133
}
133
}
134
134
135
/** Initial setup of the Model.
135
/** Initial setup of the Model.
136
  */
136
  */
137
void Model::initialize()
137
void Model::initialize()
138
{
138
{
139
   mSetup = false;
139
   mSetup = false;
140
   GlobalSettings::instance()->setCurrentYear(0);
140
   GlobalSettings::instance()->setCurrentYear(0);
141
   mGrid = 0;
141
   mGrid = 0;
142
   mHeightGrid = 0;
142
   mHeightGrid = 0;
143
   mManagement = 0;
143
   mManagement = 0;
144
   mABEManagement = 0;
144
   mABEManagement = 0;
145
   mEnvironment = 0;
145
   mEnvironment = 0;
146
   mTimeEvents = 0;
146
   mTimeEvents = 0;
147
   mStandGrid = 0;
147
   mStandGrid = 0;
148
   mModules = 0;
148
   mModules = 0;
149
   mDEM = 0;
149
   mDEM = 0;
150
   mGrassCover = 0;
150
   mGrassCover = 0;
151
   mSaplings=0;
151
   mSaplings=0;
152
}
152
}
153
153
154
/** sets up the simulation space.
154
/** sets up the simulation space.
155
*/
155
*/
156
void Model::setupSpace()
156
void Model::setupSpace()
157
{
157
{
158
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.world"));
158
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.world"));
159
    double cellSize = xml.value("cellSize", "2").toDouble();
159
    double cellSize = xml.value("cellSize", "2").toDouble();
160
    double width = xml.value("width", "100").toDouble();
160
    double width = xml.value("width", "100").toDouble();
161
    double height = xml.value("height", "100").toDouble();
161
    double height = xml.value("height", "100").toDouble();
162
    double buffer = xml.value("buffer", "60").toDouble();
162
    double buffer = xml.value("buffer", "60").toDouble();
163
    mModelRect = QRectF(0., 0., width, height);
163
    mModelRect = QRectF(0., 0., width, height);
164
164
165
    qDebug() << QString("setup of the world: %1x%2m with cell-size=%3m and %4m buffer").arg(width).arg(height).arg(cellSize).arg(buffer);
165
    qDebug() << QString("setup of the world: %1x%2m with cell-size=%3m and %4m buffer").arg(width).arg(height).arg(cellSize).arg(buffer);
166
166
167
    QRectF total_grid(QPointF(-buffer, -buffer), QPointF(width+buffer, height+buffer));
167
    QRectF total_grid(QPointF(-buffer, -buffer), QPointF(width+buffer, height+buffer));
168
    qDebug() << "setup grid rectangle:" << total_grid;
168
    qDebug() << "setup grid rectangle:" << total_grid;
169
169
170
    if (mGrid)
170
    if (mGrid)
171
        delete mGrid;
171
        delete mGrid;
172
    mGrid = new FloatGrid(total_grid, cellSize);
172
    mGrid = new FloatGrid(total_grid, cellSize);
173
    mGrid->initialize(1.f);
173
    mGrid->initialize(1.f);
174
    if (mHeightGrid)
174
    if (mHeightGrid)
175
        delete mHeightGrid;
175
        delete mHeightGrid;
176
    mHeightGrid = new HeightGrid(total_grid, cellSize*5);
176
    mHeightGrid = new HeightGrid(total_grid, cellSize*5);
177
    mHeightGrid->wipe(); // set all to zero
177
    mHeightGrid->wipe(); // set all to zero
178
    Tree::setGrid(mGrid, mHeightGrid);
178
    Tree::setGrid(mGrid, mHeightGrid);
179
179
180
    // setup the spatial location of the project area
180
    // setup the spatial location of the project area
181
    if (xml.hasNode("location")) {
181
    if (xml.hasNode("location")) {
182
        // setup of spatial location
182
        // setup of spatial location
183
        double loc_x = xml.valueDouble("location.x");
183
        double loc_x = xml.valueDouble("location.x");
184
        double loc_y = xml.valueDouble("location.y");
184
        double loc_y = xml.valueDouble("location.y");
185
        double loc_z = xml.valueDouble("location.z");
185
        double loc_z = xml.valueDouble("location.z");
186
        double loc_rot = xml.valueDouble("location.rotation");
186
        double loc_rot = xml.valueDouble("location.rotation");
187
        setupGISTransformation(loc_x, loc_y, loc_z, loc_rot);
187
        setupGISTransformation(loc_x, loc_y, loc_z, loc_rot);
188
        qDebug() << "setup of spatial location: x/y/z" << loc_x << loc_y << loc_z << "rotation:" << loc_rot;
188
        qDebug() << "setup of spatial location: x/y/z" << loc_x << loc_y << loc_z << "rotation:" << loc_rot;
189
    } else {
189
    } else {
190
        setupGISTransformation(0., 0., 0., 0.);
190
        setupGISTransformation(0., 0., 0., 0.);
191
    }
191
    }
192
192
193
    // load environment (multiple climates, speciesSets, ...
193
    // load environment (multiple climates, speciesSets, ...
194
    if (mEnvironment)
194
    if (mEnvironment)
195
        delete mEnvironment;
195
        delete mEnvironment;
196
    mEnvironment = new Environment();
196
    mEnvironment = new Environment();
197
197
198
    if (xml.valueBool("environmentEnabled", false)) {
198
    if (xml.valueBool("environmentEnabled", false)) {
199
        QString env_file = GlobalSettings::instance()->path(xml.value("environmentFile"));
199
        QString env_file = GlobalSettings::instance()->path(xml.value("environmentFile"));
200
        bool grid_mode = (xml.value("environmentMode")=="grid");
200
        bool grid_mode = (xml.value("environmentMode")=="grid");
201
        QString grid_file = GlobalSettings::instance()->path(xml.value("environmentGrid"));
201
        QString grid_file = GlobalSettings::instance()->path(xml.value("environmentGrid"));
202
        if (grid_mode) {
202
        if (grid_mode) {
203
            if (QFile::exists(grid_file))
203
            if (QFile::exists(grid_file))
204
                mEnvironment->setGridMode(grid_file);
204
                mEnvironment->setGridMode(grid_file);
205
            else
205
            else
206
                throw IException(QString("File '%1' specified in key 'environmentGrid' does not exit ('environmentMode' is 'grid').").arg(grid_file) );
206
                throw IException(QString("File '%1' specified in key 'environmentGrid' does not exit ('environmentMode' is 'grid').").arg(grid_file) );
207
        }
207
        }
208
208
209
        if (!mEnvironment->loadFromFile(env_file))
209
        if (!mEnvironment->loadFromFile(env_file))
210
            return;
210
            return;
211
    } else {
211
    } else {
212
        // load and prepare default values
212
        // load and prepare default values
213
        // (2) SpeciesSets: currently only one a global species set.
213
        // (2) SpeciesSets: currently only one a global species set.
214
        SpeciesSet *speciesSet = new SpeciesSet();
214
        SpeciesSet *speciesSet = new SpeciesSet();
215
        mSpeciesSets.push_back(speciesSet);
215
        mSpeciesSets.push_back(speciesSet);
216
        speciesSet->setup();
216
        speciesSet->setup();
217
        // Climate...
217
        // Climate...
218
        Climate *c = new Climate();
218
        Climate *c = new Climate();
219
        mClimates.push_back(c);
219
        mClimates.push_back(c);
220
        mEnvironment->setDefaultValues(c, speciesSet);
220
        mEnvironment->setDefaultValues(c, speciesSet);
221
    } // environment?
221
    } // environment?
222
222
223
    // time series data
223
    // time series data
224
    if (xml.valueBool(".timeEventsEnabled", false)) {
224
    if (xml.valueBool(".timeEventsEnabled", false)) {
225
        mTimeEvents = new TimeEvents();
225
        mTimeEvents = new TimeEvents();
226
        mTimeEvents->loadFromFile(GlobalSettings::instance()->path(xml.value("timeEventsFile"), "script"));
226
        mTimeEvents->loadFromFile(GlobalSettings::instance()->path(xml.value("timeEventsFile"), "script"));
227
    }
227
    }
228
228
229
229
230
    // simple case: create ressource units in a regular grid.
230
    // simple case: create ressource units in a regular grid.
231
    if (xml.valueBool("resourceUnitsAsGrid")) {
231
    if (xml.valueBool("resourceUnitsAsGrid")) {
232
232
233
        mRUmap.setup(QRectF(0., 0., width, height),100.); // Grid, that holds positions of resource units
233
        mRUmap.setup(QRectF(0., 0., width, height),100.); // Grid, that holds positions of resource units
234
        mRUmap.wipe();
234
        mRUmap.wipe();
235
235
236
        bool mask_is_setup = false;
236
        bool mask_is_setup = false;
237
        if (xml.valueBool("standGrid.enabled")) {
237
        if (xml.valueBool("standGrid.enabled")) {
238
            QString fileName = GlobalSettings::instance()->path(xml.value("standGrid.fileName"));
238
            QString fileName = GlobalSettings::instance()->path(xml.value("standGrid.fileName"));
239
            mStandGrid = new MapGrid(fileName,false); // create stand grid index later
239
            mStandGrid = new MapGrid(fileName,false); // create stand grid index later
240
240
241
            if (mStandGrid->isValid()) {
241
            if (mStandGrid->isValid()) {
242
                for (int i=0;i<mStandGrid->grid().count();i++) {
242
                for (int i=0;i<mStandGrid->grid().count();i++) {
243
                    const int &grid_value = mStandGrid->grid().constValueAtIndex(i);
243
                    const int &grid_value = mStandGrid->grid().constValueAtIndex(i);
244
                    mHeightGrid->valueAtIndex(i).setValid( grid_value > -1 );
244
                    mHeightGrid->valueAtIndex(i).setValid( grid_value > -1 );
245
                    if (grid_value>-1)
245
                    if (grid_value>-1)
246
                        mRUmap.valueAt(mStandGrid->grid().cellCenterPoint(i)) = (ResourceUnit*)1;
246
                        mRUmap.valueAt(mStandGrid->grid().cellCenterPoint(i)) = (ResourceUnit*)1;
247
                    if (grid_value < -1)
247
                    if (grid_value < -1)
248
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
248
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
249
                }
249
                }
250
            }
250
            }
251
            mask_is_setup = true;
251
            mask_is_setup = true;
252
        } else {
252
        } else {
253
            if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
253
            if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
254
                // in the case we have no stand grid but only a large rectangle (without the torus option)
254
                // in the case we have no stand grid but only a large rectangle (without the torus option)
255
                // we assume a forest outside
255
                // we assume a forest outside
256
                for (int i=0;i<mHeightGrid->count();++i) {
256
                for (int i=0;i<mHeightGrid->count();++i) {
257
                    const QPointF &p = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
257
                    const QPointF &p = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
258
                    if (p.x() < 0. || p.x()>width || p.y()<0. || p.y()>height) {
258
                    if (p.x() < 0. || p.x()>width || p.y()<0. || p.y()>height) {
259
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
259
                        mHeightGrid->valueAtIndex(i).setForestOutside(true);
260
                        mHeightGrid->valueAtIndex(i).setValid(false);
260
                        mHeightGrid->valueAtIndex(i).setValid(false);
261
                    }
261
                    }
262
                }
262
                }
263
263
264
            }
264
            }
265
        }
265
        }
266
266
267
        ResourceUnit **p; // ptr to ptr!
267
        ResourceUnit **p; // ptr to ptr!
268
        ResourceUnit *new_ru;
268
        ResourceUnit *new_ru;
269
269
270
        int ru_index = 0;
270
        int ru_index = 0;
271
        for (p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
271
        for (p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
272
            QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
272
            QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
273
            if (!mStandGrid || !mStandGrid->isValid() || *p>NULL) {
273
            if (!mStandGrid || !mStandGrid->isValid() || *p>NULL) {
274
                mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
274
                mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
275
                // create resource units for valid positions only
275
                // create resource units for valid positions only
276
                new_ru = new ResourceUnit(ru_index++); // create resource unit
276
                new_ru = new ResourceUnit(ru_index++); // create resource unit
277
                new_ru->setClimate( mEnvironment->climate() );
277
                new_ru->setClimate( mEnvironment->climate() );
278
                new_ru->setSpeciesSet( mEnvironment->speciesSet() );
278
                new_ru->setSpeciesSet( mEnvironment->speciesSet() );
279
                new_ru->setup();
279
                new_ru->setup();
280
                new_ru->setID( mEnvironment->currentID() ); // set id of resource unit in grid mode
280
                new_ru->setID( mEnvironment->currentID() ); // set id of resource unit in grid mode
281
                new_ru->setBoundingBox(r);
281
                new_ru->setBoundingBox(r);
282
                mRU.append(new_ru);
282
                mRU.append(new_ru);
283
                *p = new_ru; // save in the RUmap grid
283
                *p = new_ru; // save in the RUmap grid
284
            }
284
            }
285
        }
285
        }
286
        if (mEnvironment) {
286
        if (mEnvironment) {
287
            // retrieve species sets and climates (that were really used)
287
            // retrieve species sets and climates (that were really used)
288
            mSpeciesSets << mEnvironment->speciesSetList();
288
            mSpeciesSets << mEnvironment->speciesSetList();
289
            mClimates << mEnvironment->climateList();
289
            mClimates << mEnvironment->climateList();
290
            QString climate_file_list;
290
            QString climate_file_list;
291
            for (int i=0, c=0;i<mClimates.count();++i) {
291
            for (int i=0, c=0;i<mClimates.count();++i) {
292
                climate_file_list += mClimates[i]->name() + ", ";
292
                climate_file_list += mClimates[i]->name() + ", ";
293
                if (++c>5) {
293
                if (++c>5) {
294
                    climate_file_list += "...";
294
                    climate_file_list += "...";
295
                    break;
295
                    break;
296
                }
296
                }
297
297
298
            }
298
            }
299
            qDebug() << "Setup of climates: #loaded:" << mClimates.count() << "tables:" << climate_file_list;
299
            qDebug() << "Setup of climates: #loaded:" << mClimates.count() << "tables:" << climate_file_list;
300
300
301
301
302
        }
302
        }
303
303
304
        qDebug() << "setup of" << mEnvironment->climateList().size() << "climates performed.";
304
        qDebug() << "setup of" << mEnvironment->climateList().size() << "climates performed.";
305
305
306
        if (mStandGrid && mStandGrid->isValid())
306
        if (mStandGrid && mStandGrid->isValid())
307
            mStandGrid->createIndex();
307
            mStandGrid->createIndex();
308
        // now store the pointers in the grid.
308
        // now store the pointers in the grid.
309
        // Important: This has to be done after the mRU-QList is complete - otherwise pointers would
309
        // Important: This has to be done after the mRU-QList is complete - otherwise pointers would
310
        // point to invalid memory when QList's memory is reorganized (expanding)
310
        // point to invalid memory when QList's memory is reorganized (expanding)
311
//        ru_index = 0;
311
//        ru_index = 0;
312
//        for (p=mRUmap.begin();p!=mRUmap.end(); ++p) {
312
//        for (p=mRUmap.begin();p!=mRUmap.end(); ++p) {
313
//            *p = mRU.value(ru_index++);
313
//            *p = mRU.value(ru_index++);
314
//        }
314
//        }
315
        qDebug() << "created a grid of ResourceUnits: count=" << mRU.count() << "number of RU-map-cells:" << mRUmap.count();
315
        qDebug() << "created a grid of ResourceUnits: count=" << mRU.count() << "number of RU-map-cells:" << mRUmap.count();
316
316
317
317
318
        calculateStockableArea();
318
        calculateStockableArea();
319
319
320
        // setup of the project area mask
320
        // setup of the project area mask
321
        if (!mask_is_setup && xml.valueBool("areaMask.enabled", false) && xml.hasNode("areaMask.imageFile")) {
321
        if (!mask_is_setup && xml.valueBool("areaMask.enabled", false) && xml.hasNode("areaMask.imageFile")) {
322
            // to be extended!!! e.g. to load ESRI-style text files....
322
            // to be extended!!! e.g. to load ESRI-style text files....
323
            // setup a grid with the same size as the height grid...
323
            // setup a grid with the same size as the height grid...
324
            FloatGrid tempgrid((int)mHeightGrid->cellsize(), mHeightGrid->sizeX(), mHeightGrid->sizeY());
324
            FloatGrid tempgrid((int)mHeightGrid->cellsize(), mHeightGrid->sizeX(), mHeightGrid->sizeY());
325
            QString fileName = GlobalSettings::instance()->path(xml.value("areaMask.imageFile"));
325
            QString fileName = GlobalSettings::instance()->path(xml.value("areaMask.imageFile"));
326
            loadGridFromImage(fileName, tempgrid); // fetch from image
326
            loadGridFromImage(fileName, tempgrid); // fetch from image
327
            for (int i=0;i<tempgrid.count(); i++)
327
            for (int i=0;i<tempgrid.count(); i++)
328
                mHeightGrid->valueAtIndex(i).setValid( tempgrid.valueAtIndex(i)>0.99 );
328
                mHeightGrid->valueAtIndex(i).setValid( tempgrid.valueAtIndex(i)>0.99 );
329
            qDebug() << "loaded project area mask from" << fileName;
329
            qDebug() << "loaded project area mask from" << fileName;
330
        }
330
        }
331
331
332
        // list of "valid" resource units
332
        // list of "valid" resource units
333
        QList<ResourceUnit*> valid_rus;
333
        QList<ResourceUnit*> valid_rus;
334
        foreach(ResourceUnit* ru, mRU)
334
        foreach(ResourceUnit* ru, mRU)
335
            if (ru->id()!=-1)
335
            if (ru->id()!=-1)
336
                valid_rus.append(ru);
336
                valid_rus.append(ru);
337
337
338
        // setup of the digital elevation map (if present)
338
        // setup of the digital elevation map (if present)
339
        QString dem_file = xml.value("DEM");
339
        QString dem_file = xml.value("DEM");
340
        if (!dem_file.isEmpty()) {
340
        if (!dem_file.isEmpty()) {
341
            mDEM = new DEM(GlobalSettings::instance()->path(dem_file));
341
            mDEM = new DEM(GlobalSettings::instance()->path(dem_file));
342
            // add them to the visuals...
342
            // add them to the visuals...
343
            GlobalSettings::instance()->controller()->addGrid(mDEM, "DEM height", GridViewRainbow, 0, 1000);
343
            GlobalSettings::instance()->controller()->addGrid(mDEM, "DEM height", GridViewRainbow, 0, 1000);
344
            GlobalSettings::instance()->controller()->addGrid(mDEM->slopeGrid(), "DEM slope", GridViewRainbow, 0, 3);
344
            GlobalSettings::instance()->controller()->addGrid(mDEM->slopeGrid(), "DEM slope", GridViewRainbow, 0, 3);
345
            GlobalSettings::instance()->controller()->addGrid(mDEM->aspectGrid(), "DEM aspect", GridViewRainbow, 0, 360);
345
            GlobalSettings::instance()->controller()->addGrid(mDEM->aspectGrid(), "DEM aspect", GridViewRainbow, 0, 360);
346
            GlobalSettings::instance()->controller()->addGrid(mDEM->viewGrid(), "DEM view", GridViewGray, 0, 1);
346
            GlobalSettings::instance()->controller()->addGrid(mDEM->viewGrid(), "DEM view", GridViewGray, 0, 1);
347
347
348
        }
348
        }
349
349
350
        // setup of saplings
350
        // setup of saplings
351
        if (mSaplings) {
351
        if (mSaplings) {
352
            delete mSaplings; mSaplings=0;
352
            delete mSaplings; mSaplings=0;
353
        }
353
        }
354
        if (settings().regenerationEnabled) {
354
        if (settings().regenerationEnabled) {
355
            mSaplings = new Saplings();
355
            mSaplings = new Saplings();
356
            mSaplings->setup();
356
            mSaplings->setup();
357
        }
357
        }
358
358
359
359
360
        // setup of the grass cover
360
        // setup of the grass cover
361
        if (!mGrassCover)
361
        if (!mGrassCover)
362
            mGrassCover = new GrassCover();
362
            mGrassCover = new GrassCover();
363
        mGrassCover->setup();
363
        mGrassCover->setup();
364
364
365
        // setup of external modules
365
        // setup of external modules
366
        mModules->setup();
366
        mModules->setup();
367
        if (mModules->hasSetupResourceUnits()) {
367
        if (mModules->hasSetupResourceUnits()) {
368
            for (ResourceUnit **p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
368
            for (ResourceUnit **p=mRUmap.begin(); p!=mRUmap.end(); ++p) {
369
                if (*p) {
369
                if (*p) {
370
                    QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
370
                    QRectF r = mRUmap.cellRect(mRUmap.indexOf(p));
371
                    mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
371
                    mEnvironment->setPosition( r.center() ); // if environment is 'disabled' default values from the project file are used.
372
                    mModules->setupResourceUnit( *p );
372
                    mModules->setupResourceUnit( *p );
373
                }
373
                }
374
            }
374
            }
375
        }
375
        }
376
376
377
        // setup of scripting environment
377
        // setup of scripting environment
378
        ScriptGlobal::setupGlobalScripting();
378
        ScriptGlobal::setupGlobalScripting();
379
379
380
        // setup the helper that does the multithreading
380
        // setup the helper that does the multithreading
381
        threadRunner.setup(valid_rus);
381
        threadRunner.setup(valid_rus);
382
        threadRunner.setMultithreading(GlobalSettings::instance()->settings().valueBool("system.settings.multithreading"));
382
        threadRunner.setMultithreading(GlobalSettings::instance()->settings().valueBool("system.settings.multithreading"));
383
        threadRunner.print();
383
        threadRunner.print();
384
384
385
385
386
    } else  {
386
    } else  {
387
        throw IException("resourceUnitsAsGrid MUST be set to true - at least currently :)");
387
        throw IException("resourceUnitsAsGrid MUST be set to true - at least currently :)");
388
    }
388
    }
389
    mSetup = true;
389
    mSetup = true;
390
}
390
}
391
391
392
392
393
/** clear() frees all ressources allocated with the run of a simulation.
393
/** clear() frees all ressources allocated with the run of a simulation.
394

394

395
  */
395
  */
396
void Model::clear()
396
void Model::clear()
397
{
397
{
398
    mSetup = false;
398
    mSetup = false;
399
    qDebug() << "Model clear: attempting to clear" << mRU.count() << "RU, " << mSpeciesSets.count() << "SpeciesSets.";
399
    qDebug() << "Model clear: attempting to clear" << mRU.count() << "RU, " << mSpeciesSets.count() << "SpeciesSets.";
400
    // clear ressource units
400
    // clear ressource units
401
    qDeleteAll(mRU); // delete ressource units (and trees)
401
    qDeleteAll(mRU); // delete ressource units (and trees)
402
    mRU.clear();
402
    mRU.clear();
403
403
404
    qDeleteAll(mSpeciesSets); // delete species sets
404
    qDeleteAll(mSpeciesSets); // delete species sets
405
    mSpeciesSets.clear();
405
    mSpeciesSets.clear();
406
406
407
    // delete climate data
407
    // delete climate data
408
    qDeleteAll(mClimates);
408
    qDeleteAll(mClimates);
409
409
410
    // delete the grids
410
    // delete the grids
411
    if (mGrid)
411
    if (mGrid)
412
        delete mGrid;
412
        delete mGrid;
413
    if (mHeightGrid)
413
    if (mHeightGrid)
414
        delete mHeightGrid;
414
        delete mHeightGrid;
415
    if (mSaplings)
415
    if (mSaplings)
416
        delete mSaplings;
416
        delete mSaplings;
417
    if (mManagement)
417
    if (mManagement)
418
        delete mManagement;
418
        delete mManagement;
419
    if (mEnvironment)
419
    if (mEnvironment)
420
        delete mEnvironment;
420
        delete mEnvironment;
421
    if (mTimeEvents)
421
    if (mTimeEvents)
422
        delete mTimeEvents;
422
        delete mTimeEvents;
423
    if (mStandGrid)
423
    if (mStandGrid)
424
        delete mStandGrid;
424
        delete mStandGrid;
425
    if (mModules)
425
    if (mModules)
426
        delete mModules;
426
        delete mModules;
427
    if (mDEM)
427
    if (mDEM)
428
        delete mDEM;
428
        delete mDEM;
429
    if (mGrassCover)
429
    if (mGrassCover)
430
        delete mGrassCover;
430
        delete mGrassCover;
431
    if (mABEManagement)
431
    if (mABEManagement)
432
        delete mABEManagement;
432
        delete mABEManagement;
433
433
434
    mGrid = 0;
434
    mGrid = 0;
435
    mHeightGrid = 0;
435
    mHeightGrid = 0;
436
    mManagement = 0;
436
    mManagement = 0;
437
    mEnvironment = 0;
437
    mEnvironment = 0;
438
    mTimeEvents = 0;
438
    mTimeEvents = 0;
439
    mStandGrid  = 0;
439
    mStandGrid  = 0;
440
    mModules = 0;
440
    mModules = 0;
441
    mDEM = 0;
441
    mDEM = 0;
442
    mGrassCover = 0;
442
    mGrassCover = 0;
443
    mABEManagement = 0;
443
    mABEManagement = 0;
444
444
445
    GlobalSettings::instance()->outputManager()->close();
445
    GlobalSettings::instance()->outputManager()->close();
446
446
447
    qDebug() << "Model ressources freed.";
447
    qDebug() << "Model ressources freed.";
448
}
448
}
449
449
450
/** Setup of the Simulation.
450
/** Setup of the Simulation.
451
  This really creates the simulation environment and does the setup of various aspects.
451
  This really creates the simulation environment and does the setup of various aspects.
452
  */
452
  */
453
void Model::loadProject()
453
void Model::loadProject()
454
{
454
{
455
    DebugTimer dt("load project");
455
    DebugTimer dt("load project");
456
    GlobalSettings *g = GlobalSettings::instance();
456
    GlobalSettings *g = GlobalSettings::instance();
457
    g->printDirecories();
457
    g->printDirecories();
458
    const XmlHelper &xml = g->settings();
458
    const XmlHelper &xml = g->settings();
459
459
460
    g->clearDatabaseConnections();
460
    g->clearDatabaseConnections();
461
    // database connections: reset
461
    // database connections: reset
462
    GlobalSettings::instance()->clearDatabaseConnections();
462
    GlobalSettings::instance()->clearDatabaseConnections();
463
    // input and climate connection
463
    // input and climate connection
464
    // see initOutputDatabase() for output database
464
    // see initOutputDatabase() for output database
465
    QString dbPath = g->path( xml.value("system.database.in"), "database");
465
    QString dbPath = g->path( xml.value("system.database.in"), "database");
466
    GlobalSettings::instance()->setupDatabaseConnection("in", dbPath, true);
466
    GlobalSettings::instance()->setupDatabaseConnection("in", dbPath, true);
467
    dbPath = g->path( xml.value("system.database.climate"), "database");
467
    dbPath = g->path( xml.value("system.database.climate"), "database");
468
    GlobalSettings::instance()->setupDatabaseConnection("climate", dbPath, true);
468
    GlobalSettings::instance()->setupDatabaseConnection("climate", dbPath, true);
469
469
470
    mSettings.loadModelSettings();
470
    mSettings.loadModelSettings();
471
    mSettings.print();
471
    mSettings.print();
472
    // random seed: if stored value is <> 0, use this as the random seed (and produce hence always an equal sequence of random numbers)
472
    // random seed: if stored value is <> 0, use this as the random seed (and produce hence always an equal sequence of random numbers)
473
    uint seed = xml.value("system.settings.randomSeed","0").toUInt();
473
    uint seed = xml.value("system.settings.randomSeed","0").toUInt();
474
    RandomGenerator::setup(RandomGenerator::ergMersenneTwister, seed); // use the MersenneTwister as default
474
    RandomGenerator::setup(RandomGenerator::ergMersenneTwister, seed); // use the MersenneTwister as default
475
    // linearization of expressions: if true *and* linearize() is explicitely called, then
475
    // linearization of expressions: if true *and* linearize() is explicitely called, then
476
    // function results will be cached over a defined range of values.
476
    // function results will be cached over a defined range of values.
477
    bool do_linearization = xml.valueBool("system.settings.expressionLinearizationEnabled", false);
477
    bool do_linearization = xml.valueBool("system.settings.expressionLinearizationEnabled", false);
478
    Expression::setLinearizationEnabled(do_linearization);
478
    Expression::setLinearizationEnabled(do_linearization);
479
    if (do_linearization)
479
    if (do_linearization)
480
        qDebug() << "The linearization of certains expressions is enabled (performance optimization).";
480
        qDebug() << "The linearization of certains expressions is enabled (performance optimization).";
481
481
482
    // log level
482
    // log level
483
    QString log_level = xml.value("system.settings.logLevel", "debug").toLower();
483
    QString log_level = xml.value("system.settings.logLevel", "debug").toLower();
484
    if (log_level=="debug") setLogLevel(0);
484
    if (log_level=="debug") setLogLevel(0);
485
    if (log_level=="info") setLogLevel(1);
485
    if (log_level=="info") setLogLevel(1);
486
    if (log_level=="warning") setLogLevel(2);
486
    if (log_level=="warning") setLogLevel(2);
487
    if (log_level=="error") setLogLevel(3);
487
    if (log_level=="error") setLogLevel(3);
488
488
489
    // snag dynamics / soil model enabled? (info used during setup of world)
489
    // snag dynamics / soil model enabled? (info used during setup of world)
490
    changeSettings().carbonCycleEnabled = xml.valueBool("model.settings.carbonCycleEnabled", false);
490
    changeSettings().carbonCycleEnabled = xml.valueBool("model.settings.carbonCycleEnabled", false);
491
    // class size of snag classes
491
    // class size of snag classes
492
    Snag::setupThresholds(xml.valueDouble("model.settings.soil.swdDBHClass12"),
492
    Snag::setupThresholds(xml.valueDouble("model.settings.soil.swdDBHClass12"),
493
                          xml.valueDouble("model.settings.soil.swdDBHClass23"));
493
                          xml.valueDouble("model.settings.soil.swdDBHClass23"));
494
494
495
    // setup of modules
495
    // setup of modules
496
    if (mModules)
496
    if (mModules)
497
        delete mModules;
497
        delete mModules;
498
    mModules = new Modules();
498
    mModules = new Modules();
499
499
500
    changeSettings().regenerationEnabled = xml.valueBool("model.settings.regenerationEnabled", false);
500
    changeSettings().regenerationEnabled = xml.valueBool("model.settings.regenerationEnabled", false);
501
501
502
502
503
    setupSpace();
503
    setupSpace();
504
    if (mRU.isEmpty())
504
    if (mRU.isEmpty())
505
        throw IException("Setup of Model: no resource units present!");
505
        throw IException("Setup of Model: no resource units present!");
506
506
507
    // (3) additional issues
507
    // (3) additional issues
508
    // (3.1) load javascript code into the engine
508
    // (3.1) load javascript code into the engine
509
    QString script_file = xml.value("system.javascript.fileName");
509
    QString script_file = xml.value("system.javascript.fileName");
510
    if (!script_file.isEmpty()) {
510
    if (!script_file.isEmpty()) {
511
        script_file = g->path(script_file, "script");
511
        script_file = g->path(script_file, "script");
512
        ScriptGlobal::loadScript(script_file);
512
        ScriptGlobal::loadScript(script_file);
513
        g->controller()->setLoadedJavascriptFile(script_file);
513
        g->controller()->setLoadedJavascriptFile(script_file);
514
    }
514
    }
515
515
516
    // (3.2) setup of regeneration
516
    // (3.2) setup of regeneration
517
    if (settings().regenerationEnabled) {
517
    if (settings().regenerationEnabled) {
518
        foreach(SpeciesSet *ss, mSpeciesSets)
518
        foreach(SpeciesSet *ss, mSpeciesSets)
519
            ss->setupRegeneration();
519
            ss->setupRegeneration();
520
    }
520
    }
521
    Sapling::setRecruitmentVariation(xml.valueDouble("model.settings.seedDispersal.recruitmentDimensionVariation",0.1));
521
    Sapling::setRecruitmentVariation(xml.valueDouble("model.settings.seedDispersal.recruitmentDimensionVariation",0.1));
522
522
523
    // (3.3) management
523
    // (3.3) management
524
    bool use_abe = xml.valueBool("model.management.abeEnabled");
524
    bool use_abe = xml.valueBool("model.management.abeEnabled");
525
    if (use_abe) {
525
    if (use_abe) {
526
        // use the agent based forest management engine
526
        // use the agent based forest management engine
527
        mABEManagement = new ABE::ForestManagementEngine();
527
        mABEManagement = new ABE::ForestManagementEngine();
528
        // setup of ABE after loading of trees.
528
        // setup of ABE after loading of trees.
529
529
530
    }
530
    }
531
    // use the standard management
531
    // use the standard management
532
    QString mgmtFile = xml.value("model.management.file");
532
    QString mgmtFile = xml.value("model.management.file");
533
    if (!mgmtFile.isEmpty() && xml.valueBool("model.management.enabled")) {
533
    if (!mgmtFile.isEmpty() && xml.valueBool("model.management.enabled")) {
534
        mManagement = new Management();
534
        mManagement = new Management();
535
        QString path = GlobalSettings::instance()->path(mgmtFile, "script");
535
        QString path = GlobalSettings::instance()->path(mgmtFile, "script");
536
        mManagement->loadScript(path);
536
        mManagement->loadScript(path);
537
        qDebug() << "setup management using script" << path;
537
        qDebug() << "setup management using script" << path;
538
    }
538
    }
539
539
540
540
541
541
542
}
542
}
543
543
544
void Model::reloadABE()
544
void Model::reloadABE()
545
{
545
{
546
    // delete firest
546
    // delete firest
547
    if (mABEManagement)
547
    if (mABEManagement)
548
        delete mABEManagement;
548
        delete mABEManagement;
549
    mABEManagement = new ABE::ForestManagementEngine();
549
    mABEManagement = new ABE::ForestManagementEngine();
550
    // and setup
550
    // and setup
551
    mABEManagement->setup();
551
    mABEManagement->setup();
552
    mABEManagement->runOnInit(true);
552
    mABEManagement->runOnInit(true);
553
553
554
    mABEManagement->initialize();
554
    mABEManagement->initialize();
555
555
556
    mABEManagement->runOnInit(false);
556
    mABEManagement->runOnInit(false);
557
557
558
}
558
}
559
559
560
560
561
ResourceUnit *Model::ru(QPointF coord)
561
ResourceUnit *Model::ru(QPointF coord)
562
{
562
{
563
    if (!mRUmap.isEmpty() && mRUmap.coordValid(coord))
563
    if (!mRUmap.isEmpty() && mRUmap.coordValid(coord))
564
        return mRUmap.valueAt(coord);
564
        return mRUmap.valueAt(coord);
565
    return ru(); // default RU if there is only one
565
    return ru(); // default RU if there is only one
566
}
566
}
567
567
568
void Model::initOutputDatabase()
568
void Model::initOutputDatabase()
569
{
569
{
570
    GlobalSettings *g = GlobalSettings::instance();
570
    GlobalSettings *g = GlobalSettings::instance();
571
    QString dbPath = g->path(g->settings().value("system.database.out"), "output");
571
    QString dbPath = g->path(g->settings().value("system.database.out"), "output");
572
    // create run-metadata
572
    // create run-metadata
573
    int maxid = SqlHelper::queryValue("select max(id) from runs", g->dbin()).toInt();
573
    int maxid = SqlHelper::queryValue("select max(id) from runs", g->dbin()).toInt();
574
574
575
    maxid++;
575
    maxid++;
576
    QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss");
576
    QString timestamp = QDateTime::currentDateTime().toString("yyyyMMdd_hhmmss");
577
    SqlHelper::executeSql(QString("insert into runs (id, timestamp) values (%1, '%2')").arg(maxid).arg(timestamp), g->dbin());
577
    SqlHelper::executeSql(QString("insert into runs (id, timestamp) values (%1, '%2')").arg(maxid).arg(timestamp), g->dbin());
578
    // replace path information
578
    // replace path information
579
    dbPath.replace("$id$", QString::number(maxid));
579
    dbPath.replace("$id$", QString::number(maxid));
580
    dbPath.replace("$date$", timestamp);
580
    dbPath.replace("$date$", timestamp);
581
    // setup final path
581
    // setup final path
582
   g->setupDatabaseConnection("out", dbPath, false);
582
   g->setupDatabaseConnection("out", dbPath, false);
583
583
584
}
584
}
585
585
586
586
587
/// multithreaded running function for the resource unit level establishment
587
/// multithreaded running function for the resource unit level establishment
588
ResourceUnit *nc_sapling_growth_establishment(ResourceUnit *unit)
588
ResourceUnit *nc_sapling_growth_establishment(ResourceUnit *unit)
589
{
589
{
590
    try {
590
    try {
591
        { // DebugTimer t("nc_saplingGrowth"); t.setSilent();
591
        { // DebugTimer t("nc_saplingGrowth"); t.setSilent();
592
        // define a height map for the current resource unit on the stack
592
        // define a height map for the current resource unit on the stack
593
        float sapling_map[cPxPerRU*cPxPerRU];
593
        float sapling_map[cPxPerRU*cPxPerRU];
594
        // set the map and initialize it:
594
        // set the map and initialize it:
595
        unit->setSaplingHeightMap(sapling_map);
595
        unit->setSaplingHeightMap(sapling_map);
596
596
597
597
598
        // (1) calculate the growth of (already established) saplings (populate sapling map)
598
        // (1) calculate the growth of (already established) saplings (populate sapling map)
599
        QList<ResourceUnitSpecies*>::const_iterator rus;
599
        QList<ResourceUnitSpecies*>::const_iterator rus;
600
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
600
        for (rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
601
            (*rus)->calclulateSaplingGrowth();
601
            (*rus)->calclulateSaplingGrowth();
602
602
603
        } { // DebugTimer t("nc_Establishment"); t.setSilent();
603
        } { // DebugTimer t("nc_Establishment"); t.setSilent();
604
604
605
        // (2) calculate the establishment probabilities of new saplings
605
        // (2) calculate the establishment probabilities of new saplings
606
        for (QList<ResourceUnitSpecies*>::const_iterator rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
606
        for (QList<ResourceUnitSpecies*>::const_iterator rus=unit->ruSpecies().cbegin(); rus!=unit->ruSpecies().cend(); ++rus)
607
            (*rus)->calculateEstablishment();
607
            (*rus)->calculateEstablishment();
608
608
609
        }
609
        }
610
610
611
    } catch (const IException& e) {
611
    } catch (const IException& e) {
612
        GlobalSettings::instance()->controller()->throwError(e.message());
612
        GlobalSettings::instance()->controller()->throwError(e.message());
613
    }
613
    }
614
614
615
    unit->setSaplingHeightMap(0); // invalidate again
615
    unit->setSaplingHeightMap(0); // invalidate again
616
    return unit;
616
    return unit;
617
617
618
}
618
}
619
619
620
620
621
/// multithreaded execution of the carbon cycle routine
621
/// multithreaded execution of the carbon cycle routine
622
ResourceUnit *nc_carbonCycle(ResourceUnit *unit)
622
ResourceUnit *nc_carbonCycle(ResourceUnit *unit)
623
{
623
{
624
    try {
624
    try {
625
        // (1) do calculations on snag dynamics for the resource unit
625
        // (1) do calculations on snag dynamics for the resource unit
626
        unit->calculateCarbonCycle();
626
        unit->calculateCarbonCycle();
627
        // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N)
627
        // (2) do the soil carbon and nitrogen dynamics calculations (ICBM/2N)
628
    } catch (const IException& e) {
628
    } catch (const IException& e) {
629
        GlobalSettings::instance()->controller()->throwError(e.message());
629
        GlobalSettings::instance()->controller()->throwError(e.message());
630
    }
630
    }
631
631
632
    return unit;
632
    return unit;
633
}
633
}
634
634
635
/// beforeRun performs several steps before the models starts running.
635
/// beforeRun performs several steps before the models starts running.
636
/// inter alia: * setup of the stands
636
/// inter alia: * setup of the stands
637
///             * setup of the climates
637
///             * setup of the climates
638
void Model::beforeRun()
638
void Model::beforeRun()
639
{
639
{
640
    // setup outputs
640
    // setup outputs
641
    // setup output database
641
    // setup output database
642
    if (GlobalSettings::instance()->dbout().isOpen())
642
    if (GlobalSettings::instance()->dbout().isOpen())
643
        GlobalSettings::instance()->dbout().close();
643
        GlobalSettings::instance()->dbout().close();
644
    initOutputDatabase();
644
    initOutputDatabase();
645
    GlobalSettings::instance()->outputManager()->setup();
645
    GlobalSettings::instance()->outputManager()->setup();
646
    GlobalSettings::instance()->clearDebugLists();
646
    GlobalSettings::instance()->clearDebugLists();
647
647
648
    // initialize stands
648
    // initialize stands
649
    StandLoader loader(this);
649
    StandLoader loader(this);
650
    {
650
    {
651
    DebugTimer loadtrees("load trees");
651
    DebugTimer loadtrees("load trees");
652
    loader.processInit();
652
    loader.processInit();
653
    }
653
    }
654
    // initalization of ABE
654
    // initalization of ABE
655
    if (mABEManagement) {
655
    if (mABEManagement) {
656
        mABEManagement->setup();
656
        mABEManagement->setup();
657
        mABEManagement->runOnInit(true);
657
        mABEManagement->runOnInit(true);
658
    }
658
    }
659
659
660
    // load climate
660
    // load climate
661
    {
661
    {
662
        if (logLevelDebug()) qDebug() << "attempting to load climate..." ;
662
        if (logLevelDebug()) qDebug() << "attempting to load climate..." ;
663
        DebugTimer loadclim("load climate");
663
        DebugTimer loadclim("load climate");
664
        foreach(Climate *c, mClimates) {
664
        foreach(Climate *c, mClimates) {
665
            if (!c->isSetup())
665
            if (!c->isSetup())
666
                c->setup();
666
                c->setup();
667
        }
667
        }
668
        // load the first year of the climate database
668
        // load the first year of the climate database
669
        foreach(Climate *c, mClimates)
669
        foreach(Climate *c, mClimates)
670
            c->nextYear();
670
            c->nextYear();
671
671
672
    }
672
    }
673
673
674
674
675
    { DebugTimer loadinit("load standstatistics");
675
    { DebugTimer loadinit("load standstatistics");
676
    if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ;
676
    if (logLevelDebug()) qDebug() << "attempting to calculate initial stand statistics (incl. apply and read pattern)..." ;
677
    Tree::setGrid(mGrid, mHeightGrid);
677
    Tree::setGrid(mGrid, mHeightGrid);
678
    // debugCheckAllTrees(); // introduced for debugging session (2012-04-06)
678
    // debugCheckAllTrees(); // introduced for debugging session (2012-04-06)
679
    applyPattern();
679
    applyPattern();
680
    readPattern();
680
    readPattern();
681
    loader.processAfterInit(); // e.g. initialization of saplings
681
    loader.processAfterInit(); // e.g. initialization of saplings
682
682
683
    // force the compilation of initial stand statistics
683
    // force the compilation of initial stand statistics
684
    createStandStatistics();
684
    createStandStatistics();
685
    }
685
    }
686
686
687
    // initalization of ABE (now all stands are properly set up)
687
    // initalization of ABE (now all stands are properly set up)
688
    if (mABEManagement) {
688
    if (mABEManagement) {
689
        mABEManagement->initialize();
689
        mABEManagement->initialize();
690
        mABEManagement->runOnInit(false);
690
        mABEManagement->runOnInit(false);
691
    }
691
    }
692
692
693
693
694
    // outputs to create with inital state (without any growth) are called here:
694
    // outputs to create with inital state (without any growth) are called here:
695
    GlobalSettings::instance()->setCurrentYear(0); // set clock to "0" (for outputs with initial state)
695
    GlobalSettings::instance()->setCurrentYear(0); // set clock to "0" (for outputs with initial state)
696
696
697
    GlobalSettings::instance()->outputManager()->execute("stand"); // year=0
697
    GlobalSettings::instance()->outputManager()->execute("stand"); // year=0
698
    GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0
698
    GlobalSettings::instance()->outputManager()->execute("landscape"); // year=0
699
    GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0
699
    GlobalSettings::instance()->outputManager()->execute("sapling"); // year=0
700
    GlobalSettings::instance()->outputManager()->execute("tree"); // year=0
700
    GlobalSettings::instance()->outputManager()->execute("tree"); // year=0
701
    GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0
701
    GlobalSettings::instance()->outputManager()->execute("dynamicstand"); // year=0
702
702
703
    GlobalSettings::instance()->setCurrentYear(1); // set to first year
703
    GlobalSettings::instance()->setCurrentYear(1); // set to first year
704
704
705
}
705
}
706
706
707
/** Main model runner.
707
/** Main model runner.
708
  The sequence of actions is as follows:
708
  The sequence of actions is as follows:
709
  (1) Load the climate of the new year
709
  (1) Load the climate of the new year
710
  (2) Reset statistics for resource unit as well as for dead/managed trees
710
  (2) Reset statistics for resource unit as well as for dead/managed trees
711
  (3) Invoke Management.
711
  (3) Invoke Management.
712
  (4) *after* that, calculate Light patterns
712
  (4) *after* that, calculate Light patterns
713
  (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)
713
  (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)
714
  (6) execute Regeneration
714
  (6) execute Regeneration
715
  (7) invoke disturbance modules
715
  (7) invoke disturbance modules
716
  (8) calculate statistics for the year
716
  (8) calculate statistics for the year
717
  (9) write database outputs
717
  (9) write database outputs
718
  */
718
  */
719
void Model::runYear()
719
void Model::runYear()
720
{
720
{
721
    DebugTimer t("Model::runYear()");
721
    DebugTimer t("Model::runYear()");
722
    GlobalSettings::instance()->systemStatistics()->reset();
722
    GlobalSettings::instance()->systemStatistics()->reset();
723
    RandomGenerator::checkGenerator(); // see if we need to generate new numbers...
723
    RandomGenerator::checkGenerator(); // see if we need to generate new numbers...
724
    // initalization at start of year for external modules
724
    // initalization at start of year for external modules
725
    mModules->yearBegin();
725
    mModules->yearBegin();
726
726
727
    // execute scheduled events for the current year
727
    // execute scheduled events for the current year
728
    if (mTimeEvents)
728
    if (mTimeEvents)
729
        mTimeEvents->run();
729
        mTimeEvents->run();
730
730
731
    // load the next year of the climate database (except for the first year - the first climate year is loaded immediately
731
    // load the next year of the climate database (except for the first year - the first climate year is loaded immediately
732
    if (GlobalSettings::instance()->currentYear()>1) {
732
    if (GlobalSettings::instance()->currentYear()>1) {
733
        foreach(Climate *c, mClimates)
733
        foreach(Climate *c, mClimates)
734
            c->nextYear();
734
            c->nextYear();
735
    }
735
    }
736
736
737
    // reset statistics
737
    // reset statistics
738
    foreach(ResourceUnit *ru, mRU)
738
    foreach(ResourceUnit *ru, mRU)
739
        ru->newYear();
739
        ru->newYear();
740
740
741
    foreach(SpeciesSet *set, mSpeciesSets)
741
    foreach(SpeciesSet *set, mSpeciesSets)
742
        set->newYear();
742
        set->newYear();
743
    // management classic
743
    // management classic
744
    if (mManagement) {
744
    if (mManagement) {
745
        DebugTimer t("management");
745
        DebugTimer t("management");
746
        mManagement->run();
746
        mManagement->run();
747
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
747
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
748
    }
748
    }
749
    // ... or ABE (the agent based variant)
749
    // ... or ABE (the agent based variant)
750
    if (mABEManagement) {
750
    if (mABEManagement) {
751
        DebugTimer t("ABE:run");
751
        DebugTimer t("ABE:run");
752
        mABEManagement->run();
752
        mABEManagement->run();
753
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
753
        GlobalSettings::instance()->systemStatistics()->tManagement+=t.elapsed();
754
    }
754
    }
755
755
756
    // if trees are dead/removed because of management, the tree lists
756
    // if trees are dead/removed because of management, the tree lists
757
    // need to be cleaned (and the statistics need to be recreated)
757
    // need to be cleaned (and the statistics need to be recreated)
758
    cleanTreeLists();
758
    cleanTreeLists();
759
759
760
    // process a cycle of individual growth
760
    // process a cycle of individual growth
761
    applyPattern(); // create Light Influence Patterns
761
    applyPattern(); // create Light Influence Patterns
762
    readPattern(); // readout light state of individual trees
762
    readPattern(); // readout light state of individual trees
763
    grow(); // let the trees grow (growth on stand-level, tree-level, mortality)
763
    grow(); // let the trees grow (growth on stand-level, tree-level, mortality)
764
    mGrassCover->execute(); // evaluate the grass / herb cover (and its effect on regeneration)
764
    mGrassCover->execute(); // evaluate the grass / herb cover (and its effect on regeneration)
765
765
766
    // regeneration
766
    // regeneration
767
    if (settings().regenerationEnabled) {
767
    if (settings().regenerationEnabled) {
768
        // seed dispersal
768
        // seed dispersal
769
        DebugTimer tseed("Regeneration and Establishment");
769
        DebugTimer tseed("Regeneration and Establishment");
770
        foreach(SpeciesSet *set, mSpeciesSets)
770
        foreach(SpeciesSet *set, mSpeciesSets)
771
            set->regeneration(); // parallel execution for each species set
771
            set->regeneration(); // parallel execution for each species set
772
772
773
        GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed();
773
        GlobalSettings::instance()->systemStatistics()->tSeedDistribution+=tseed.elapsed();
774
        // establishment
774
        // establishment
775
        Sapling::updateBrowsingPressure();
775
        Sapling::updateBrowsingPressure();
776
776
777
        { DebugTimer t("saplingGrowthEstablishment");
777
        { DebugTimer t("saplingGrowthEstablishment");
778
        executePerResourceUnit( nc_sapling_growth_establishment, false /* true: force single thraeded operation */);
778
        executePerResourceUnit( nc_sapling_growth_establishment, false /* true: force single thraeded operation */);
779
        GlobalSettings::instance()->systemStatistics()->tSaplingAndEstablishment+=t.elapsed();
779
        GlobalSettings::instance()->systemStatistics()->tSaplingAndEstablishment+=t.elapsed();
780
        }
780
        }
781
781
782
        Establishment::debugInfo(); // debug test
782
        Establishment::debugInfo(); // debug test
783
783
784
    }
784
    }
785
785
786
    // calculate soil / snag dynamics
786
    // calculate soil / snag dynamics
787
    if (settings().carbonCycleEnabled) {
787
    if (settings().carbonCycleEnabled) {
788
        DebugTimer ccycle("carbon cylce");
788
        DebugTimer ccycle("carbon cylce");
789
        executePerResourceUnit( nc_carbonCycle, false /* true: force single threaded operation */);
789
        executePerResourceUnit( nc_carbonCycle, false /* true: force single threaded operation */);
790
        GlobalSettings::instance()->systemStatistics()->tCarbonCycle+=ccycle.elapsed();
790
        GlobalSettings::instance()->systemStatistics()->tCarbonCycle+=ccycle.elapsed();
791
791
792
    }
792
    }
793
793
794
    // external modules/disturbances
794
    // external modules/disturbances
795
    mModules->run();
795
    mModules->run();
796
796
797
    // cleanup of tree lists if external modules removed trees.
797
    // cleanup of tree lists if external modules removed trees.
798
    cleanTreeLists();
798
    cleanTreeLists();
799
799
800
800
801
    DebugTimer toutput("outputs");
801
    DebugTimer toutput("outputs");
802
    // calculate statistics
802
    // calculate statistics
803
    foreach(ResourceUnit *ru, mRU)
803
    foreach(ResourceUnit *ru, mRU)
804
        ru->yearEnd();
804
        ru->yearEnd();
805
805
806
    // create outputs
806
    // create outputs
807
    OutputManager *om = GlobalSettings::instance()->outputManager();
807
    OutputManager *om = GlobalSettings::instance()->outputManager();
808
    om->execute("tree"); // single tree output
808
    om->execute("tree"); // single tree output
809
    om->execute("treeremoval"); // single removed tree output
809
    om->execute("treeremoval"); // single removed tree output
810
    om->execute("stand"); //resource unit level x species
810
    om->execute("stand"); //resource unit level x species
811
    om->execute("landscape"); //landscape x species
811
    om->execute("landscape"); //landscape x species
-
 
812
    om->execute("landscape_removed"); //removed trees on landscape x species
812
    om->execute("sapling"); // sapling layer per RU x species
813
    om->execute("sapling"); // sapling layer per RU x species
813
    om->execute("production_month"); // 3pg responses growth per species x RU x month
814
    om->execute("production_month"); // 3pg responses growth per species x RU x month
814
    om->execute("dynamicstand"); // output with user-defined columns (based on species x RU)
815
    om->execute("dynamicstand"); // output with user-defined columns (based on species x RU)
815
    om->execute("standdead"); // resource unit level x species
816
    om->execute("standdead"); // resource unit level x species
816
    om->execute("management"); // resource unit level x species
817
    om->execute("management"); // resource unit level x species
817
    om->execute("carbon"); // resource unit level, carbon pools above and belowground
818
    om->execute("carbon"); // resource unit level, carbon pools above and belowground
818
    om->execute("carbonflow"); // resource unit level, GPP, NPP and total carbon flows (atmosphere, harvest, ...)
819
    om->execute("carbonflow"); // resource unit level, GPP, NPP and total carbon flows (atmosphere, harvest, ...)
819
820
820
    GlobalSettings::instance()->systemStatistics()->tWriteOutput+=toutput.elapsed();
821
    GlobalSettings::instance()->systemStatistics()->tWriteOutput+=toutput.elapsed();
821
    GlobalSettings::instance()->systemStatistics()->tTotalYear+=t.elapsed();
822
    GlobalSettings::instance()->systemStatistics()->tTotalYear+=t.elapsed();
822
    GlobalSettings::instance()->systemStatistics()->writeOutput();
823
    GlobalSettings::instance()->systemStatistics()->writeOutput();
823
824
824
    GlobalSettings::instance()->setCurrentYear(GlobalSettings::instance()->currentYear()+1);
825
    GlobalSettings::instance()->setCurrentYear(GlobalSettings::instance()->currentYear()+1);
825
826
826
    // try to clean up a bit of memory (useful if many large JS objects (e.g., grids) are used)
827
    // try to clean up a bit of memory (useful if many large JS objects (e.g., grids) are used)
827
    GlobalSettings::instance()->scriptEngine()->collectGarbage();
828
    GlobalSettings::instance()->scriptEngine()->collectGarbage();
828
}
829
}
829
830
830
831
831
832
832
void Model::afterStop()
833
void Model::afterStop()
833
{
834
{
834
    // do some cleanup
835
    // do some cleanup
835
}
836
}
836
837
837
/// multithreaded running function for LIP printing
838
/// multithreaded running function for LIP printing
838
ResourceUnit* nc_applyPattern(ResourceUnit *unit)
839
ResourceUnit* nc_applyPattern(ResourceUnit *unit)
839
{
840
{
840
841
841
    QVector<Tree>::iterator tit;
842
    QVector<Tree>::iterator tit;
842
    QVector<Tree>::iterator tend = unit->trees().end();
843
    QVector<Tree>::iterator tend = unit->trees().end();
843
844
844
    try {
845
    try {
845
846
846
        // light concurrence influence
847
        // light concurrence influence
847
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
848
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
848
            // height dominance grid
849
            // height dominance grid
849
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
850
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
850
                (*tit).heightGrid(); // just do it ;)
851
                (*tit).heightGrid(); // just do it ;)
851
852
852
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
853
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
853
                (*tit).applyLIP(); // just do it ;)
854
                (*tit).applyLIP(); // just do it ;)
854
855
855
        } else {
856
        } else {
856
            // height dominance grid
857
            // height dominance grid
857
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
858
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
858
                (*tit).heightGrid_torus(); // just do it ;)
859
                (*tit).heightGrid_torus(); // just do it ;)
859
860
860
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
861
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
861
                (*tit).applyLIP_torus(); // do it the wraparound way
862
                (*tit).applyLIP_torus(); // do it the wraparound way
862
        }
863
        }
863
        return unit;
864
        return unit;
864
    } catch (const IException &e) {
865
    } catch (const IException &e) {
865
        GlobalSettings::instance()->controller()->throwError(e.message());
866
        GlobalSettings::instance()->controller()->throwError(e.message());
866
    }
867
    }
867
    return unit;
868
    return unit;
868
}
869
}
869
870
870
/// multithreaded running function for LIP value extraction
871
/// multithreaded running function for LIP value extraction
871
ResourceUnit *nc_readPattern(ResourceUnit *unit)
872
ResourceUnit *nc_readPattern(ResourceUnit *unit)
872
{
873
{
873
    QVector<Tree>::iterator tit;
874
    QVector<Tree>::iterator tit;
874
    QVector<Tree>::iterator  tend = unit->trees().end();
875
    QVector<Tree>::iterator  tend = unit->trees().end();
875
    try {
876
    try {
876
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
877
        if (!GlobalSettings::instance()->settings().paramValueBool("torus")) {
877
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
878
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
878
                (*tit).readLIF(); // multipliactive approach
879
                (*tit).readLIF(); // multipliactive approach
879
        } else {
880
        } else {
880
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
881
            for (tit=unit->trees().begin(); tit!=tend; ++tit)
881
                (*tit).readLIF_torus(); // do it the wraparound way
882
                (*tit).readLIF_torus(); // do it the wraparound way
882
        }
883
        }
883
    } catch (const IException &e) {
884
    } catch (const IException &e) {
884
        GlobalSettings::instance()->controller()->throwError(e.message());
885
        GlobalSettings::instance()->controller()->throwError(e.message());
885
    }
886
    }
886
    return unit;
887
    return unit;
887
}
888
}
888
889
889
/// multithreaded running function for the growth of individual trees
890
/// multithreaded running function for the growth of individual trees
890
ResourceUnit *nc_grow(ResourceUnit *unit)
891
ResourceUnit *nc_grow(ResourceUnit *unit)
891
{
892
{
892
    QVector<Tree>::iterator tit;
893
    QVector<Tree>::iterator tit;
893
    QVector<Tree>::iterator  tend = unit->trees().end();
894
    QVector<Tree>::iterator  tend = unit->trees().end();
894
    try {
895
    try {
895
        unit->beforeGrow(); // reset statistics
896
        unit->beforeGrow(); // reset statistics
896
        // calculate light responses
897
        // calculate light responses
897
        // responses are based on *modified* values for LightResourceIndex
898
        // responses are based on *modified* values for LightResourceIndex
898
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
899
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
899
            (*tit).calcLightResponse();
900
            (*tit).calcLightResponse();
900
        }
901
        }
901
902
902
        unit->calculateInterceptedArea();
903
        unit->calculateInterceptedArea();
903
904
904
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
905
        for (tit=unit->trees().begin(); tit!=tend; ++tit) {
905
            (*tit).grow(); // actual growth of individual trees
906
            (*tit).grow(); // actual growth of individual trees
906
        }
907
        }
907
    } catch (const IException &e) {
908
    } catch (const IException &e) {
908
        GlobalSettings::instance()->controller()->throwError(e.message());
909
        GlobalSettings::instance()->controller()->throwError(e.message());
909
    }
910
    }
910
911
911
    GlobalSettings::instance()->systemStatistics()->treeCount+=unit->trees().count();
912
    GlobalSettings::instance()->systemStatistics()->treeCount+=unit->trees().count();
912
    return unit;
913
    return unit;
913
}
914
}
914
915
915
/// multithreaded running function for the resource level production
916
/// multithreaded running function for the resource level production
916
ResourceUnit *nc_production(ResourceUnit *unit)
917
ResourceUnit *nc_production(ResourceUnit *unit)
917
{
918
{
918
    try {
919
    try {
919
        unit->production();
920
        unit->production();
920
    } catch (const IException &e) {
921
    } catch (const IException &e) {
921
        GlobalSettings::instance()->controller()->throwError(e.message());
922
        GlobalSettings::instance()->controller()->throwError(e.message());
922
    }
923
    }
923
    return unit;
924
    return unit;
924
}
925
}
925
926
926
927
927
void Model::test()
928
void Model::test()
928
{
929
{
929
    // Test-funktion: braucht 1/3 time von readGrid()
930
    // Test-funktion: braucht 1/3 time von readGrid()
930
    DebugTimer t("test");
931
    DebugTimer t("test");
931
    FloatGrid averaged = mGrid->averaged(10);
932
    FloatGrid averaged = mGrid->averaged(10);
932
    int count = 0;
933
    int count = 0;
933
    float *end = averaged.end();
934
    float *end = averaged.end();
934
    for (float *p=averaged.begin(); p!=end; ++p)
935
    for (float *p=averaged.begin(); p!=end; ++p)
935
        if (*p > 0.9)
936
        if (*p > 0.9)
936
            count++;
937
            count++;
937
    qDebug() << count << "LIF>0.9 of " << averaged.count();
938
    qDebug() << count << "LIF>0.9 of " << averaged.count();
938
}
939
}
939
940
940
void Model::debugCheckAllTrees()
941
void Model::debugCheckAllTrees()
941
{
942
{
942
    AllTreeIterator at(this);
943
    AllTreeIterator at(this);
943
    bool has_errors = false; double dummy=0.;
944
    bool has_errors = false; double dummy=0.;
944
    while (Tree *t = at.next()) {
945
    while (Tree *t = at.next()) {
945
        // plausibility
946
        // plausibility
946
        if (t->dbh()<0 || t->dbh()>10000. || t->biomassFoliage()<0. || t->height()>1000. || t->height() < 0.
947
        if (t->dbh()<0 || t->dbh()>10000. || t->biomassFoliage()<0. || t->height()>1000. || t->height() < 0.
947
                || t->biomassFoliage() <0.)
948
                || t->biomassFoliage() <0.)
948
            has_errors = true;
949
            has_errors = true;
949
        // check for objects....
950
        // check for objects....
950
        dummy = t->stamp()->offset() + t->ru()->ruSpecies()[1]->statistics().count();
951
        dummy = t->stamp()->offset() + t->ru()->ruSpecies()[1]->statistics().count();
951
    }
952
    }
952
    if (has_errors)
953
    if (has_errors)
953
        qDebug() << "model: debugCheckAllTrees found problems" << dummy;
954
        qDebug() << "model: debugCheckAllTrees found problems" << dummy;
954
}
955
}
955
956
956
void Model::applyPattern()
957
void Model::applyPattern()
957
{
958
{
958
959
959
    DebugTimer t("applyPattern()");
960
    DebugTimer t("applyPattern()");
960
    // intialize grids...
961
    // intialize grids...
961
    initializeGrid();
962
    initializeGrid();
962
963
963
    // initialize height grid with a value of 4m. This is the height of the regeneration layer
964
    // initialize height grid with a value of 4m. This is the height of the regeneration layer
964
    for (HeightGridValue *h=mHeightGrid->begin();h!=mHeightGrid->end();++h) {
965
    for (HeightGridValue *h=mHeightGrid->begin();h!=mHeightGrid->end();++h) {
965
        h->resetCount(); // set count = 0, but do not touch the flags
966
        h->resetCount(); // set count = 0, but do not touch the flags
966
        h->height = 4.f;
967
        h->height = 4.f;
967
    }
968
    }
968
969
969
    threadRunner.run(nc_applyPattern);
970
    threadRunner.run(nc_applyPattern);
970
    GlobalSettings::instance()->systemStatistics()->tApplyPattern+=t.elapsed();
971
    GlobalSettings::instance()->systemStatistics()->tApplyPattern+=t.elapsed();
971
}
972
}
972
973
973
void Model::readPattern()
974
void Model::readPattern()
974
{
975
{
975
    DebugTimer t("readPattern()");
976
    DebugTimer t("readPattern()");
976
    threadRunner.run(nc_readPattern);
977
    threadRunner.run(nc_readPattern);
977
    GlobalSettings::instance()->systemStatistics()->tReadPattern+=t.elapsed();
978
    GlobalSettings::instance()->systemStatistics()->tReadPattern+=t.elapsed();
978
979
979
}
980
}
980
981
981
/** Main function for the growth of stands and trees.
982
/** Main function for the growth of stands and trees.
982
   This includes several steps.
983
   This includes several steps.
983
   (1) calculate the stocked area (i.e. count pixels in height grid)
984
   (1) calculate the stocked area (i.e. count pixels in height grid)
984
   (2) 3PG production (including response calculation, water cycle)
985
   (2) 3PG production (including response calculation, water cycle)
985
   (3) single tree growth (including mortality)
986
   (3) single tree growth (including mortality)
986
   (4) cleanup of tree lists (remove dead trees)
987
   (4) cleanup of tree lists (remove dead trees)
987
  */
988
  */
988
void Model::grow()
989
void Model::grow()
989
{
990
{
990
991
991
    if (!settings().growthEnabled)
992
    if (!settings().growthEnabled)
992
        return;
993
        return;
993
    { DebugTimer t("growRU()");
994
    { DebugTimer t("growRU()");
994
    calculateStockedArea();
995
    calculateStockedArea();
995
996
996
    // multithreaded: mutex for the message handler in mainwindow solved the crashes.
997
    // multithreaded: mutex for the message handler in mainwindow solved the crashes.
997
    threadRunner.run(nc_production);
998
    threadRunner.run(nc_production);
998
    }
999
    }
999
1000
1000
    DebugTimer t("growTrees()");
1001
    DebugTimer t("growTrees()");
1001
    threadRunner.run(nc_grow); // actual growth of individual trees
1002
    threadRunner.run(nc_grow); // actual growth of individual trees
1002
1003
1003
    foreach(ResourceUnit *ru, mRU) {
1004
    foreach(ResourceUnit *ru, mRU) {
1004
       ru->cleanTreeList();
1005
       ru->cleanTreeList();
1005
       ru->afterGrow();
1006
       ru->afterGrow();
1006
       //qDebug() << (b-n) << "trees died (of" << b << ").";
1007
       //qDebug() << (b-n) << "trees died (of" << b << ").";
1007
   }
1008
   }
1008
   GlobalSettings::instance()->systemStatistics()->tTreeGrowth+=t.elapsed();
1009
   GlobalSettings::instance()->systemStatistics()->tTreeGrowth+=t.elapsed();
1009
}
1010
}
1010
1011
1011
/** calculate for each resource unit the fraction of area which is stocked.
1012
/** calculate for each resource unit the fraction of area which is stocked.
1012
  This is done by checking the pixels of the global height grid.
1013
  This is done by checking the pixels of the global height grid.
1013
  */
1014
  */
1014
void Model::calculateStockedArea()
1015
void Model::calculateStockedArea()
1015
{
1016
{
1016
    // iterate over the whole heightgrid and count pixels for each ressource unit
1017
    // iterate over the whole heightgrid and count pixels for each ressource unit
1017
    HeightGridValue *end = mHeightGrid->end();
1018
    HeightGridValue *end = mHeightGrid->end();
1018
    QPointF cp;
1019
    QPointF cp;
1019
    ResourceUnit *ru;
1020
    ResourceUnit *ru;
1020
    for (HeightGridValue *i=mHeightGrid->begin(); i!=end; ++i) {
1021
    for (HeightGridValue *i=mHeightGrid->begin(); i!=end; ++i) {
1021
        cp = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
1022
        cp = mHeightGrid->cellCenterPoint(mHeightGrid->indexOf(i));
1022
        if (mRUmap.coordValid(cp)) {
1023
        if (mRUmap.coordValid(cp)) {
1023
            ru = mRUmap.valueAt(cp);
1024
            ru = mRUmap.valueAt(cp);
1024
            if (ru) {
1025
            if (ru) {
1025
                ru->countStockedPixel( (*i).count()>0 );
1026
                ru->countStockedPixel( (*i).count()>0 );
1026
            }
1027
            }
1027
        }
1028
        }
1028
1029
1029
    }
1030
    }
1030
}
1031
}
1031
1032
1032
/** calculate for each resource unit the stockable area.
1033
/** calculate for each resource unit the stockable area.
1033
  "stockability" is determined by the isValid flag of resource units which in turn
1034
  "stockability" is determined by the isValid flag of resource units which in turn
1034
  is derived from stand grid values.
1035
  is derived from stand grid values.
1035
  */
1036
  */
1036
void Model::calculateStockableArea()
1037
void Model::calculateStockableArea()
1037
{
1038
{
1038
1039
-
 
1040
    mTotalStockableArea = 0.;
1039
    foreach(ResourceUnit *ru, mRU) {
1041
    foreach(ResourceUnit *ru, mRU) {
1040
        // //
1042
        // //
1041
        //        if (ru->id()==-1) {
1043
        //        if (ru->id()==-1) {
1042
        //            ru->setStockableArea(0.);
1044
        //            ru->setStockableArea(0.);
1043
        //            continue;
1045
        //            continue;
1044
        //        }
1046
        //        }
1045
        GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
1047
        GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
1046
        int valid=0, total=0;
1048
        int valid=0, total=0;
1047
        while (runner.next()) {
1049
        while (runner.next()) {
1048
            if ( runner.current()->isValid() )
1050
            if ( runner.current()->isValid() )
1049
                valid++;
1051
                valid++;
1050
            total++;
1052
            total++;
1051
        }
1053
        }
1052
        if (total) {
1054
        if (total) {
1053
            ru->setStockableArea( cHeightPixelArea * valid);
-
 
-
 
1055
            ru->setStockableArea( cHeightPixelArea * valid); // in m2
-
 
1056
            mTotalStockableArea += cHeightPixelArea * valid / 10000.; // in ha
1054
            if (valid==0 && ru->id()>-1) {
1057
            if (valid==0 && ru->id()>-1) {
1055
                // invalidate this resource unit
1058
                // invalidate this resource unit
1056
                ru->setID(-1);
1059
                ru->setID(-1);
1057
            }
1060
            }
1058
            if (valid>0 && ru->id()==-1) {
1061
            if (valid>0 && ru->id()==-1) {
1059
                qDebug() << "Warning: a resource unit has id=-1 but stockable area (id was set to 0)!!! ru: " << ru->boundingBox() << "with index" << ru->index();
1062
                qDebug() << "Warning: a resource unit has id=-1 but stockable area (id was set to 0)!!! ru: " << ru->boundingBox() << "with index" << ru->index();
1060
                ru->setID(0);
1063
                ru->setID(0);
1061
                // test-code
1064
                // test-code
1062
                //GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
1065
                //GridRunner<HeightGridValue> runner(*mHeightGrid, ru->boundingBox());
1063
                //while (runner.next()) {
1066
                //while (runner.next()) {
1064
                //    qDebug() << mHeightGrid->cellCenterPoint(mHeightGrid->indexOf( runner.current() )) << ": " << runner.current()->isValid();
1067
                //    qDebug() << mHeightGrid->cellCenterPoint(mHeightGrid->indexOf( runner.current() )) << ": " << runner.current()->isValid();
1065
                //}
1068
                //}
1066
1069
1067
            }
1070
            }
1068
        } else
1071
        } else
1069
            throw IException("calculateStockableArea: resource unit without pixels!");
1072
            throw IException("calculateStockableArea: resource unit without pixels!");
1070
1073
1071
    }
1074
    }
1072
    // mark those pixels that are at the edge of a "forest-out-of-area"
1075
    // mark those pixels that are at the edge of a "forest-out-of-area"
1073
    GridRunner<HeightGridValue> runner(*mHeightGrid, mHeightGrid->metricRect());
1076
    GridRunner<HeightGridValue> runner(*mHeightGrid, mHeightGrid->metricRect());
1074
    HeightGridValue* neighbors[8];
1077
    HeightGridValue* neighbors[8];
1075
    while (runner.next()) {
1078
    while (runner.next()) {
1076
        if (runner.current()->isForestOutside()) {
1079
        if (runner.current()->isForestOutside()) {
1077
            // if the current pixel is a "radiating" border pixel,
1080
            // if the current pixel is a "radiating" border pixel,
1078
            // then check the neighbors and set a flag if the pixel is a neighbor of a in-project-area pixel.
1081
            // then check the neighbors and set a flag if the pixel is a neighbor of a in-project-area pixel.
1079
            runner.neighbors8(neighbors);
1082
            runner.neighbors8(neighbors);
1080
            for (int i=0;i<8;++i)
1083
            for (int i=0;i<8;++i)
1081
                if (neighbors[i] &&  neighbors[i]->isValid())
1084
                if (neighbors[i] &&  neighbors[i]->isValid())
1082
                    runner.current()->setIsRadiating();
1085
                    runner.current()->setIsRadiating();
1083
1086
1084
        }
1087
        }
1085
    }
1088
    }
-
 
1089
-
 
1090
    qDebug() << "Total stockable area of the landscape is" << mTotalStockableArea << "ha.";
1086
1091
1087
}
1092
}
1088
1093
1089
void Model::initializeGrid()
1094
void Model::initializeGrid()
1090
{
1095
{
1091
    // fill the whole grid with a value of "1."
1096
    // fill the whole grid with a value of "1."
1092
    mGrid->initialize(1.f);
1097
    mGrid->initialize(1.f);
1093
1098
1094
    // apply special values for grid cells border regions where out-of-area cells
1099
    // apply special values for grid cells border regions where out-of-area cells
1095
    // radiate into the main LIF grid.
1100
    // radiate into the main LIF grid.
1096
    QPoint p;
1101
    QPoint p;
1097
    int ix_min, ix_max, iy_min, iy_max, ix_center, iy_center;
1102
    int ix_min, ix_max, iy_min, iy_max, ix_center, iy_center;
1098
    const int px_offset = cPxPerHeight / 2; // for 5 px per height grid cell, the offset is 2
1103
    const int px_offset = cPxPerHeight / 2; // for 5 px per height grid cell, the offset is 2
1099
    const int max_radiate_distance = 7;
1104
    const int max_radiate_distance = 7;
1100
    const float step_width = 1.f / (float)max_radiate_distance;
1105
    const float step_width = 1.f / (float)max_radiate_distance;
1101
    int c_rad = 0;
1106
    int c_rad = 0;
1102
    for (HeightGridValue *hgv=mHeightGrid->begin(); hgv!=mHeightGrid->end(); ++hgv) {
1107
    for (HeightGridValue *hgv=mHeightGrid->begin(); hgv!=mHeightGrid->end(); ++hgv) {
1103
        if (hgv->isRadiating()) {
1108
        if (hgv->isRadiating()) {
1104
            p=mHeightGrid->indexOf(hgv);
1109
            p=mHeightGrid->indexOf(hgv);
1105
            ix_min = p.x() * cPxPerHeight - max_radiate_distance + px_offset;
1110
            ix_min = p.x() * cPxPerHeight - max_radiate_distance + px_offset;
1106
            ix_max = ix_min + 2*max_radiate_distance + 1;
1111
            ix_max = ix_min + 2*max_radiate_distance + 1;
1107
            ix_center = ix_min + max_radiate_distance;
1112
            ix_center = ix_min + max_radiate_distance;
1108
            iy_min = p.y() * cPxPerHeight - max_radiate_distance + px_offset;
1113
            iy_min = p.y() * cPxPerHeight - max_radiate_distance + px_offset;
1109
            iy_max = iy_min + 2*max_radiate_distance + 1;
1114
            iy_max = iy_min + 2*max_radiate_distance + 1;
1110
            iy_center = iy_min + max_radiate_distance;
1115
            iy_center = iy_min + max_radiate_distance;
1111
            for (int y=iy_min; y<=iy_max; ++y) {
1116
            for (int y=iy_min; y<=iy_max; ++y) {
1112
                for (int x=ix_min; x<=ix_max; ++x) {
1117
                for (int x=ix_min; x<=ix_max; ++x) {
1113
                    if (!mGrid->isIndexValid(x,y) ||  !(*mHeightGrid)(x/cPxPerHeight, y/cPxPerHeight).isValid())
1118
                    if (!mGrid->isIndexValid(x,y) ||  !(*mHeightGrid)(x/cPxPerHeight, y/cPxPerHeight).isValid())
1114
                        continue;
1119
                        continue;
1115
                    float value = qMax(qAbs(x-ix_center), qAbs(y-iy_center)) * step_width;
1120
                    float value = qMax(qAbs(x-ix_center), qAbs(y-iy_center)) * step_width;
1116
                    float &v = mGrid->valueAtIndex(x, y);
1121
                    float &v = mGrid->valueAtIndex(x, y);
1117
                    if (value>=0.f && v>value)
1122
                    if (value>=0.f && v>value)
1118
                        v = value;
1123
                        v = value;
1119
                }
1124
                }
1120
            }
1125
            }
1121
            c_rad++;
1126
            c_rad++;
1122
        }
1127
        }
1123
    }
1128
    }
1124
    if (logLevelDebug())
1129
    if (logLevelDebug())
1125
        qDebug() << "initialize grid:" << c_rad << "radiating pixels...";
1130
        qDebug() << "initialize grid:" << c_rad << "radiating pixels...";
1126
1131
1127
}
1132
}
1128
1133
1129
1134
1130
/// Force the creation of stand statistics.
1135
/// Force the creation of stand statistics.
1131
/// - stocked area (for resourceunit-areas)
1136
/// - stocked area (for resourceunit-areas)
1132
/// - ru - statistics
1137
/// - ru - statistics
1133
void Model::createStandStatistics()
1138
void Model::createStandStatistics()
1134
{
1139
{
1135
    calculateStockedArea();
1140
    calculateStockedArea();
1136
    foreach(ResourceUnit *ru, mRU) {
1141
    foreach(ResourceUnit *ru, mRU) {
1137
        ru->addTreeAgingForAllTrees();
1142
        ru->addTreeAgingForAllTrees();
1138
        ru->createStandStatistics();
1143
        ru->createStandStatistics();
1139
    }
1144
    }
1140
}
1145
}
1141
1146
1142
void Model::cleanTreeLists()
1147
void Model::cleanTreeLists()
1143
{
1148
{
1144
    foreach(ResourceUnit *ru, GlobalSettings::instance()->model()->ruList()) {
1149
    foreach(ResourceUnit *ru, GlobalSettings::instance()->model()->ruList()) {
1145
        if (ru->hasDiedTrees()) {
1150
        if (ru->hasDiedTrees()) {
1146
            ru->cleanTreeList();
1151
            ru->cleanTreeList();
1147
            ru->recreateStandStatistics();
1152
            ru->recreateStandStatistics();
1148
        }
1153
        }
1149
    }
1154
    }
1150
}
1155
}
1151
1156
1152
1157
1153
 
1158