Subversion Repositories public iLand

Rev

Rev 971 | Rev 1057 | 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
 
106 Werner 21
#include "global.h"
22
#include "standloader.h"
23
 
290 werner 24
 
106 Werner 25
#include "grid.h"
26
#include "model.h"
189 iland 27
#include "resourceunit.h"
106 Werner 28
#include "speciesset.h"
29
 
30
#include "helper.h"
290 werner 31
#include "random.h"
137 Werner 32
#include "expression.h"
140 Werner 33
#include "expressionwrapper.h"
284 werner 34
#include "environment.h"
287 werner 35
#include "csvfile.h"
549 werner 36
#include "mapgrid.h"
699 werner 37
#include "snapshot.h"
106 Werner 38
 
697 werner 39
/** @class StandLoader
40
    @ingroup tools
41
    loads (initializes) trees for a "stand" from various sources.
42
    StandLoader initializes trees on the landscape. It reads (usually) from text files, creates the
43
    trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid).
44
 
45
    See http://iland.boku.ac.at/initialize+trees
46
  */
381 werner 47
// provide a mapping between "Picus"-style and "iLand"-style species Ids
384 werner 48
QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
155 werner 49
QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";
732 werner 50
 
290 werner 51
StandLoader::~StandLoader()
52
{
53
    if (mRandom)
54
        delete mRandom;
732 werner 55
    if (mHeightGridResponse)
56
        delete mHeightGridResponse;
290 werner 57
}
106 Werner 58
 
59
 
281 werner 60
void StandLoader::copyTrees()
61
{
62
    // we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards
63
    const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();
64
    ResourceUnit **p = ruGrid.begin();
65
    if (!p)
66
        throw IException("Standloader: invalid resource unit pointer!");
67
    ++p; // skip the first...
68
    const QVector<Tree> &tocopy = mModel->ru()->trees();
69
    for (; p!=ruGrid.end(); ++p) {
70
        QRectF rect = (*p)->boundingBox();
71
        foreach(const Tree& tree, tocopy) {
72
            Tree &newtree = (*p)->newTree();
73
            newtree = tree; // copy tree data...
739 werner 74
            newtree.setPosition(tree.position()+rect.topLeft());
281 werner 75
            newtree.setRU(*p);
76
            newtree.setNewId();
77
        }
78
    }
431 werner 79
    if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";
281 werner 80
}
81
 
284 werner 82
/** main routine of the stand setup.
83
*/
106 Werner 84
void StandLoader::processInit()
85
{
86
    GlobalSettings *g = GlobalSettings::instance();
192 werner 87
    XmlHelper xml(g->settings().node("model.initialization"));
191 werner 88
 
281 werner 89
    QString copy_mode = xml.value("mode", "copy");
90
    QString type = xml.value("type", "");
192 werner 91
    QString  fileName = xml.value("file", "");
106 Werner 92
 
732 werner 93
    bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);
736 werner 94
    mHeightGridTries = xml.valueDouble("heightGrid.maxTries", 10.);
732 werner 95
    QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function
96
    if (height_grid_enabled) {
97
        QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init");
98
        qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file;
99
 
100
        QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));
743 werner 101
        if (!p->isValid()) {
102
            throw IException(QString("Error when loading grid with tree heights for stand initalization: file %1 not found or not valid.").arg(init_height_grid_file));
103
        }
732 werner 104
        height_grid.swap(p);
105
        mInitHeightGrid = height_grid.data();
106
 
107
        QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");
108
        mHeightGridResponse = new Expression(expr);
109
        mHeightGridResponse->setLinearizationEnabled(true);
110
    }
111
 
106 Werner 112
    Tree::resetStatistics();
281 werner 113
 
114
    // one global init-file for the whole area:
284 werner 115
    if (copy_mode=="single") {
281 werner 116
        loadInitFile(fileName, type);
299 werner 117
        evaluateDebugTrees();
281 werner 118
        return;
284 werner 119
    }
281 werner 120
 
137 Werner 121
 
281 werner 122
    // call a single tree init for each resource unit
123
    if (copy_mode=="unit") {
284 werner 124
        foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
125
            ResourceUnit *ru = const_cast<ResourceUnit*>(const_ru);
126
            // set environment
127
            g->model()->environment()->setPosition(ru->boundingBox().center());
128
            type = xml.value("type", "");
129
            fileName = xml.value("file", "");
319 werner 130
            if (fileName.isEmpty())
131
                continue;
549 werner 132
            loadInitFile(fileName, type, 0, ru);
431 werner 133
            if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
284 werner 134
        }
299 werner 135
        evaluateDebugTrees();
284 werner 136
        return;
281 werner 137
    }
299 werner 138
 
732 werner 139
    // map-modus: load a init file for each "polygon" in the standgrid
549 werner 140
    if (copy_mode=="map") {
141
        if (!g->model()->standGrid() || !g->model()->standGrid()->isValid())
142
            throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));
143
        QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");
144
 
145
        CSVFile map_file(map_file_name);
146
        if (map_file.rowCount()==0)
147
            throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));
148
        int ikey = map_file.columnIndex("id");
149
        int ivalue = map_file.columnIndex("filename");
150
        if (ikey<0 || ivalue<0)
151
            throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));
152
        QString file_name;
153
        for (int i=0;i<map_file.rowCount();i++) {
154
            int key = map_file.value(i, ikey).toInt();
155
            if (key>0) {
156
                file_name = map_file.value(i, ivalue).toString();
550 werner 157
                if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;
734 werner 158
                if (!file_name.isEmpty())
159
                    loadInitFile(file_name, type, key, NULL);
549 werner 160
            }
161
        }
732 werner 162
        mInitHeightGrid = 0;
549 werner 163
        evaluateDebugTrees();
164
        return;
165
    }
904 werner 166
 
167
    // standgrid mode: load one large init file
168
    if (copy_mode=="standgrid") {
169
        fileName = GlobalSettings::instance()->path(fileName, "init");
170
        if (!QFile::exists(fileName))
171
            throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
172
        QString content = Helper::loadTextFile(fileName);
173
        // this processes the init file (also does the checking) and
174
        // stores in a QHash datastrucutre
175
        parseInitFile(content, fileName);
176
 
177
        // setup the random distribution
178
        QString density_func = xml.value("model.initialization.randomFunction", "1-x^2");
179
        if (logLevelInfo())  qDebug() << "density function:" << density_func;
180
        if (!mRandom || (mRandom->densityFunction()!= density_func)) {
181
            if (mRandom)
182
                delete mRandom;
183
            mRandom=new RandomCustomPDF(density_func);
184
            if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
185
        }
186
 
187
        if (mStandInitItems.isEmpty())
188
            throw IException("StandLoader::processInit: 'mode' is 'standgrid' but the init file is either empty or contains no 'stand_id'-column.");
189
        QHash<int, QVector<InitFileItem> >::const_iterator it = mStandInitItems.constBegin();
190
        while (it!=mStandInitItems.constEnd()) {
191
            mInitItems = it.value(); // copy the items...
192
            executeiLandInitStand(it.key());
193
            ++it;
194
        }
195
        qDebug() << "finished setup of trees.";
196
        return;
197
 
198
    }
699 werner 199
    if (copy_mode=="snapshot") {
200
        // load a snapshot from a file
201
        Snapshot shot;
549 werner 202
 
699 werner 203
        QString input_db = GlobalSettings::instance()->path(fileName);
204
        shot.loadSnapshot(input_db);
205
        return;
206
    }
281 werner 207
    throw IException("StandLoader::processInit: invalid initalization.mode!");
208
}
209
 
967 werner 210
void StandLoader::processAfterInit()
211
{
212
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.initialization"));
213
 
214
    QString mode = xml.value("mode", "copy");
215
    if (mode=="standgrid") {
216
        // load a file with saplings per polygon
217
        QString  filename = xml.value("saplingFile", "");
218
        if (filename.isEmpty())
219
            return;
220
        filename = GlobalSettings::instance()->path(filename, "init");
221
        if (!QFile::exists(filename))
222
            throw IException(QString("load-sapling-ini-file: file '%1' does not exist.").arg(filename));
223
        CSVFile init_file(filename);
224
        int istandid = init_file.columnIndex("stand_id");
225
        if (istandid==-1)
226
            throw IException("Sapling-Init: the init file contains no 'stand_id' column (required in 'standgrid' mode).");
227
 
228
        int stand_id = -99999;
229
        int ilow = -1, ihigh = 0;
230
        int total = 0;
231
        for (int i=0;i<init_file.rowCount();++i) {
232
            int row_stand = init_file.value(i, istandid).toInt();
233
            if (row_stand != stand_id) {
234
                if (stand_id>=0) {
235
                    // process stand
236
                    ihigh = i-1; // up to the last
237
                    total += loadSaplingsLIF(stand_id, init_file, ilow, ihigh);
238
                }
239
                ilow = i; // mark beginning of new stand
240
                stand_id = row_stand;
241
            }
242
        }
243
        if (stand_id>=0)
244
            total += loadSaplingsLIF(stand_id, init_file, ilow, init_file.rowCount()-1); // the last stand
245
        qDebug() << "initialization of sapling: total created:" << total;
246
 
247
    }
248
 
249
}
250
 
281 werner 251
void StandLoader::evaluateDebugTrees()
252
{
137 Werner 253
    // evaluate debugging
254
    QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
299 werner 255
    int counter=0;
137 Werner 256
    if (!dbg_str.isEmpty()) {
733 werner 257
        if (dbg_str == "debugstamp") {
734 werner 258
            qDebug() << "debug_tree = debugstamp: try touching all trees...";
733 werner 259
            // try to force an error if a stamp is invalid
260
            AllTreeIterator at(GlobalSettings::instance()->model());
753 werner 261
            double total_offset=0.;
733 werner 262
            while (Tree *t=at.next()) {
263
                total_offset += t->stamp()->offset();
734 werner 264
                if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
265
                    qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
733 werner 266
            }
734 werner 267
            qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
733 werner 268
            return;
269
        }
270
        TreeWrapper tw;
271
        Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
137 Werner 272
        AllTreeIterator at(GlobalSettings::instance()->model());
273
        double result;
274
        while (Tree *t = at.next()) {
140 Werner 275
            tw.setTree(t);
137 Werner 276
            result = dexp.execute();
299 werner 277
            if (result) {
137 Werner 278
                t->enableDebugging();
299 werner 279
                counter++;
280
            }
137 Werner 281
        }
732 werner 282
        qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
137 Werner 283
    }
106 Werner 284
}
285
 
904 werner 286
 
732 werner 287
/// load a single init file. Calls loadPicusFile() or loadiLandFile()
288
/// @param fileName file to load
289
/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
549 werner 290
int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
281 werner 291
{
284 werner 292
    QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
293
    if (!QFile::exists(pathFileName))
294
        throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
281 werner 295
 
549 werner 296
    if (type=="picus" || type=="single") {
297
        if (stand_id>0)
298
            throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
294 werner 299
        return loadPicusFile(pathFileName, ru);
549 werner 300
    }
384 werner 301
    if (type=="iland" || type=="distribution")
549 werner 302
        return loadiLandFile(pathFileName, ru, stand_id);
281 werner 303
 
384 werner 304
    throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
281 werner 305
}
306
 
393 werner 307
int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
106 Werner 308
{
294 werner 309
    QString content = Helper::loadTextFile(fileName);
310
    if (content.isEmpty()) {
311
        qDebug() << "file not found: " + fileName;
393 werner 312
        return 0;
294 werner 313
    }
393 werner 314
    return loadSingleTreeList(content, ru, fileName);
294 werner 315
}
316
 
389 werner 317
/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.
318
  returns the number of loaded trees.
319
  */
320
int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
294 werner 321
{
106 Werner 322
    if (!ru)
323
        ru = mModel->ru();
284 werner 324
    Q_ASSERT(ru!=0);
106 Werner 325
 
284 werner 326
    QPointF offset = ru->boundingBox().topLeft();
327
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
328
 
384 werner 329
    QString my_content(content);
330
    // cut out the <trees> </trees> part if present
331
    if (content.contains("<trees>")) {
332
        QRegExp rx(".*<trees>(.*)</trees>.*");
333
        rx.indexIn(content, 0);
334
        if (rx.capturedTexts().count()<1)
389 werner 335
            return 0;
384 werner 336
        my_content = rx.cap(1).trimmed();
337
    }
338
 
948 werner 339
    CSVFile infile;
340
    infile.loadFromString(my_content);
384 werner 341
 
138 Werner 342
 
948 werner 343
    int iID =  infile.columnIndex("id");
344
    int iX = infile.columnIndex("x");
345
    int iY = infile.columnIndex("y");
346
    int iBhd = infile.columnIndex("bhdfrom");
384 werner 347
    if (iBhd<0)
948 werner 348
        iBhd = infile.columnIndex("dbh");
384 werner 349
    double height_conversion = 100.;
948 werner 350
    int iHeight =infile.columnIndex("treeheight");
384 werner 351
    if (iHeight<0) {
948 werner 352
        iHeight = infile.columnIndex("height");
384 werner 353
        height_conversion = 1.; // in meter
354
    }
948 werner 355
    int iSpecies = infile.columnIndex("species");
356
    int iAge = infile.columnIndex("age");
204 werner 357
    if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1)
948 werner 358
        throw IException(QString("Initfile %1 is not valid!\nRequired columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
106 Werner 359
 
138 Werner 360
    double dbh;
384 werner 361
    bool ok;
389 werner 362
    int cnt=0;
384 werner 363
    QString speciesid;
948 werner 364
    for (int i=1;i<infile.rowCount();i++) {
365
        dbh = infile.value(i, iBhd).toDouble();
285 werner 366
 
948 werner 367
        //if (dbh<5.)
368
        //    continue;
369
 
106 Werner 370
        QPointF f;
371
        if (iX>=0 && iY>=0) {
948 werner 372
           f.setX( infile.value(i, iX).toDouble() );
373
           f.setY( infile.value(i, iY).toDouble() );
106 Werner 374
           f+=offset;
285 werner 375
 
106 Werner 376
        }
285 werner 377
        // position valid?
378
        if (!mModel->heightGrid()->valueAt(f).isValid())
379
            continue;
380
        Tree &tree = ru->newTree();
381
        tree.setPosition(f);
247 werner 382
        if (iID>=0)
948 werner 383
            tree.setId(infile.value(i, iID).toInt());
106 Werner 384
 
138 Werner 385
        tree.setDbh(dbh);
948 werner 386
        tree.setHeight(infile.value(i,iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
384 werner 387
 
948 werner 388
        speciesid = infile.value(i,iSpecies).toString();
270 werner 389
        int picusid = speciesid.toInt(&ok);
390
        if (ok) {
384 werner 391
            int idx = picusSpeciesIds.indexOf(picusid);
270 werner 392
            if (idx==-1)
393
                throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
138 Werner 394
            speciesid = iLandSpeciesIds[idx];
270 werner 395
        }
138 Werner 396
        Species *s = speciesSet->species(speciesid);
397
        if (!ru || !s)
389 werner 398
            throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
388 werner 399
        tree.setSpecies(s);
107 Werner 400
 
388 werner 401
        ok = true;
402
        if (iAge>=0)
948 werner 403
           tree.setAge(infile.value(i, iAge).toInt(&ok), tree.height()); // this is a *real* age
388 werner 404
        if (iAge<0 || !ok || tree.age()==0)
405
           tree.setAge(0, tree.height()); // no real tree age available
406
 
138 Werner 407
        tree.setRU(ru);
106 Werner 408
        tree.setup();
389 werner 409
        cnt++;
106 Werner 410
    }
389 werner 411
    return cnt;
106 Werner 412
    //qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
413
    //qDebug() << "lines: " << lines;
414
}
287 werner 415
 
393 werner 416
/** initialize trees on a resource unit based on dbh distributions.
417
  use a fairly clever algorithm to determine tree positions.
418
  see http://iland.boku.ac.at/initialize+trees
419
  @param content tree init file (including headers) in a string
420
  @param ru resource unit
421
  @param fileName source file name (for error reporting)
422
  @return number of trees added
423
  */
549 werner 424
int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
287 werner 425
{
904 werner 426
    int total_count = parseInitFile(content, fileName, ru);
427
    if (total_count==0)
428
        return 0;
429
 
430
 
431
    // setup the random distribution
432
    QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
433
    if (logLevelInfo())  qDebug() << "density function:" << density_func;
434
    if (!mRandom || (mRandom->densityFunction()!= density_func)) {
435
        if (mRandom)
436
            delete mRandom;
437
        mRandom=new RandomCustomPDF(density_func);
438
        if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
439
    }
440
    if (stand_id>0) {
441
        // execute stand based initialization
442
        executeiLandInitStand(stand_id);
443
    } else {
444
        // exeucte the initialization based on single resource units
445
        executeiLandInit(ru);
446
        ru->cleanTreeList();
447
    }
448
    return total_count;
449
 
450
}
451
 
452
int StandLoader::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru)
453
{
287 werner 454
    if (!ru)
455
        ru = mModel->ru();
456
    Q_ASSERT(ru!=0);
457
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
458
    Q_ASSERT(speciesSet!=0);
459
 
431 werner 460
    //DebugTimer t("StandLoader::loadiLandFile");
294 werner 461
    CSVFile infile;
462
    infile.loadFromString(content);
463
 
287 werner 464
    int icount = infile.columnIndex("count");
465
    int ispecies = infile.columnIndex("species");
466
    int idbh_from = infile.columnIndex("dbh_from");
467
    int idbh_to = infile.columnIndex("dbh_to");
468
    int ihd = infile.columnIndex("hd");
469
    int iage = infile.columnIndex("age");
311 werner 470
    int idensity = infile.columnIndex("density");
287 werner 471
    if (icount<0 || ispecies<0 || idbh_from<0 || idbh_to<0 || ihd<0 || iage<0)
472
        throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));
473
 
904 werner 474
    int istandid = infile.columnIndex("stand_id");
294 werner 475
    mInitItems.clear();
904 werner 476
    mStandInitItems.clear();
477
 
287 werner 478
    InitFileItem item;
384 werner 479
    bool ok;
393 werner 480
    int total_count = 0;
287 werner 481
    for (int row=0;row<infile.rowCount();row++) {
968 werner 482
        item.count = infile.value(row, icount).toDouble();
549 werner 483
        total_count += item.count;
484
        item.dbh_from = infile.value(row, idbh_from).toDouble();
485
        item.dbh_to = infile.value(row, idbh_to).toDouble();
486
        item.hd = infile.value(row, ihd).toDouble();
775 werner 487
        if (item.hd==0. || item.dbh_from / 100. * item.hd < 4.)
782 werner 488
            qWarning() << QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ;
489
            //throw IException(QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) );
549 werner 490
        ok = true;
491
        if (iage>=0)
492
            item.age = infile.value(row, iage).toInt(&ok);
493
        if (iage<0 || !ok)
494
            item.age = 0;
384 werner 495
 
549 werner 496
        item.species = speciesSet->species(infile.value(row, ispecies).toString());
497
        if (idensity>=0)
498
            item.density = infile.value(row, idensity).toDouble();
499
        else
500
            item.density = 0.;
971 werner 501
        if (item.density<-1)
549 werner 502
            throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
503
                             .arg(item.density)
504
                             .arg(fileName)
505
                             .arg(row));
506
        if (!item.species) {
507
            throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
508
                             .arg(infile.value(row, ispecies).toString())
509
                             .arg(fileName)
510
                             .arg(row));
511
        }
904 werner 512
        if (istandid>=0) {
513
            int standid = infile.value(row,istandid).toInt();
514
            mStandInitItems[standid].push_back(item);
515
        } else {
516
            mInitItems.push_back(item);
517
        }
287 werner 518
    }
393 werner 519
    return total_count;
287 werner 520
 
521
}
522
 
904 werner 523
 
549 werner 524
int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)
294 werner 525
{
526
    if (!QFile::exists(fileName))
527
        throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
528
    QString content = Helper::loadTextFile(fileName);
549 werner 529
    return loadDistributionList(content, ru, stand_id, fileName);
294 werner 530
}
531
 
287 werner 532
// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
533
// e.g. 12 = centerpixel, 0: upper left corner, ...
534
int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
535
                     1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
290 werner 536
int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
537
                       6,18,2,10,4,24,12,0,8,14,20,22};
538
 
539
 
549 werner 540
 
287 werner 541
// sort function
542
bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
311 werner 543
{
544
    return s1.second < s2.second;
545
}
546
 
549 werner 547
struct SInitPixel {
732 werner 548
    double basal_area; // accumulated basal area
549
    QPoint pixelOffset; // location of the pixel
550
    ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to
551
    double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)
971 werner 552
    bool locked; // pixel is dedicated to a single species
553
    SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.), locked(false) {}
549 werner 554
};
555
 
556
bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)
557
{
558
    return s1.basal_area < s2.basal_area;
559
}
560
 
971 werner 561
bool sortInitPixelUnlocked(const SInitPixel &s1, const SInitPixel &s2)
562
{
563
    return !s1.locked && s2.locked;
564
}
549 werner 565
 
311 werner 566
/**
567
*/
549 werner 568
 
290 werner 569
void StandLoader::executeiLandInit(ResourceUnit *ru)
287 werner 570
{
290 werner 571
 
287 werner 572
    QPointF offset = ru->boundingBox().topLeft();
573
    QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);
574
 
575
    // a multimap holds a list for all trees.
576
    // key is the index of a 10x10m pixel within the resource unit
577
    QMultiMap<int, int> tree_map;
549 werner 578
    //QHash<int,SInitPixel> tcount;
579
 
287 werner 580
    QVector<QPair<int, double> > tcount; // counts
581
    for (int i=0;i<100;i++)
582
        tcount.push_back(QPair<int,double>(i,0.));
583
 
584
    int key;
311 werner 585
    double rand_val, rand_fraction;
312 werner 586
    int total_count = 0;
290 werner 587
    foreach(const InitFileItem &item, mInitItems) {
311 werner 588
        rand_fraction = fabs(double(item.density));
287 werner 589
        for (int i=0;i<item.count;i++) {
590
            // create trees
591
            int tree_idx = ru->newTreeIndex();
592
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
707 werner 593
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
287 werner 594
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
595
            tree.setSpecies(item.species);
384 werner 596
            if (item.age<=0)
388 werner 597
                tree.setAge(0,tree.height());
381 werner 598
            else
388 werner 599
                tree.setAge(item.age, tree.height());
287 werner 600
            tree.setRU(ru);
601
            tree.setup();
312 werner 602
            total_count++;
311 werner 603
 
604
            // calculate random value. "density" is from 1..-1.
605
            rand_val = mRandom->get();
606
            if (item.density<0)
607
                rand_val = 1. - rand_val;
705 werner 608
            rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
311 werner 609
 
287 werner 610
            // key: rank of target pixel
611
            // first: index of target pixel
612
            // second: sum of target pixel
311 werner 613
            key = limit(int(100*rand_val), 0, 99); // get from random number generator
287 werner 614
            tree_map.insert(tcount[key].first, tree_idx); // store tree in map
615
            tcount[key].second+=tree.basalArea(); // aggregate the basal area for each 10m pixel
312 werner 616
            if ( (total_count < 20 && i%2==0)
617
                || (total_count<100 && i%10==0 )
618
                || (i%30==0) ) {
287 werner 619
                qSort(tcount.begin(), tcount.end(), sortPairLessThan);
312 werner 620
            }
287 werner 621
        }
622
        qSort(tcount.begin(), tcount.end(), sortPairLessThan);
623
    }
624
 
625
    int bits, index, pos;
626
    int c;
627
    QList<int> trees;
628
    QPoint tree_pos;
289 werner 629
 
287 werner 630
    for (int i=0;i<100;i++) {
631
        trees = tree_map.values(i);
632
        c = trees.count();
290 werner 633
        QPointF pixel_center = ru->boundingBox().topLeft() + QPointF((i/10)*10. + 5., (i%10)*10. + 5.);
634
        if (!mModel->heightGrid()->valueAt(pixel_center).isValid()) {
635
            // no trees on that pixel: let trees die
636
            foreach(int tree_idx, trees) {
637
                ru->trees()[tree_idx].die();
638
            }
639
            continue;
640
        }
641
 
287 werner 642
        bits = 0;
643
        index = -1;
290 werner 644
        double r;
287 werner 645
        foreach(int tree_idx, trees) {
646
            if (c>18) {
647
                index = (index + 1)%25;
648
            } else {
649
                int stop=1000;
290 werner 650
                index = 0;
287 werner 651
                do {
290 werner 652
                    //r = drandom();
653
                    //if (r<0.5)  // skip position with a prob. of 50% -> adds a little "noise"
654
                    //    index++;
655
                    //index = (index + 1)%25; // increase and roll over
656
 
657
                    // search a random position
658
                    r = drandom();
659
                    index = limit(int(25 *  r*r), 0, 24); // use rnd()^2 to search for locations -> higher number of low indices (i.e. 50% of lookups in first 25% of locations)
287 werner 660
                } while (isBitSet(bits, index)==true && stop--);
661
                if (!stop)
662
                    qDebug() << "executeiLandInit: found no free bit.";
663
                setBit(bits, index, true); // mark position as used
664
            }
290 werner 665
            // get position from fixed lists (one for even, one for uneven resource units)
666
            pos = ru->index()%2?evenlist[index]:unevenlist[index];
287 werner 667
            tree_pos = offsetIdx  // position of resource unit
668
                       + QPoint(5*(i/10), 5*(i%10)) // relative position of 10x10m pixel
669
                       + QPoint(pos/5, pos%5); // relative position within 10x10m pixel
670
            //qDebug() << tree_no++ << "to" << index;
671
            ru->trees()[tree_idx].setPosition(tree_pos);
672
        }
673
    }
674
}
549 werner 675
 
676
 
550 werner 677
 
678
// Initialization routine based on a stand map.
679
// Basically a list of 10m pixels for a given stand is retrieved
680
// and the filled with the same procedure as the resource unit based init
681
// see http://iland.boku.ac.at/initialize+trees
549 werner 682
void StandLoader::executeiLandInitStand(int stand_id)
683
{
684
 
685
    const MapGrid *grid = GlobalSettings::instance()->model()->standGrid();
901 werner 686
    if (mCurrentMap)
687
        grid = mCurrentMap;
549 werner 688
 
732 werner 689
    // get a list of positions of all pixels that belong to our stand
549 werner 690
    QList<int> indices = grid->gridIndices(stand_id);
550 werner 691
    if (indices.isEmpty()) {
692
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
693
        return;
549 werner 694
    }
695
    // a multiHash holds a list for all trees.
696
    // key is the location of the 10x10m pixel
697
    QMultiHash<QPoint, int> tree_map;
698
    QList<SInitPixel> pixel_list; // working list of all 10m pixels
732 werner 699
    pixel_list.reserve(indices.size());
549 werner 700
 
701
    foreach (int i, indices) {
702
       SInitPixel p;
703
       p.pixelOffset = grid->grid().indexOf(i); // index in the 10m grid
732 werner 704
       p.resource_unit = GlobalSettings::instance()->model()->ru( grid->grid().cellCenterPoint(p.pixelOffset));
705
       if (mInitHeightGrid)
706
           p.h_max = mInitHeightGrid->grid().constValueAtIndex(p.pixelOffset);
549 werner 707
       pixel_list.append(p);
708
    }
709
    double area_factor = grid->area(stand_id) / 10000.;
710
 
732 werner 711
    int key=0;
549 werner 712
    double rand_val, rand_fraction;
713
    int total_count = 0;
732 werner 714
    int total_tries = 0;
715
    int total_misses = 0;
716
    if (mInitHeightGrid && !mHeightGridResponse)
717
        throw IException("executeiLandInitStand: trying to initialize with height grid but without response function.");
971 werner 718
 
719
    Species *last_locked_species=0;
549 werner 720
    foreach(const InitFileItem &item, mInitItems) {
971 werner 721
        if (item.density>1.) {
722
            // special case with single-species-area
723
            if (total_count==0) {
724
                // randomize the pixels
725
                for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it)
726
                    it->basal_area = drandom();
727
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
728
                for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it)
729
                    it->basal_area = 0.;
730
            }
731
 
732
            if (item.species != last_locked_species) {
733
                last_locked_species=item.species;
734
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelUnlocked);
735
            }
736
        } else {
737
            qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
738
            last_locked_species=0;
739
        }
740
        rand_fraction = item.density;
549 werner 741
        int count = item.count * area_factor + 0.5; // round
732 werner 742
        double init_max_height = item.dbh_to/100. * item.hd;
549 werner 743
        for (int i=0;i<count;i++) {
744
 
732 werner 745
            bool found = false;
746
            int tries = mHeightGridTries;
747
            while (!found &&--tries) {
748
                // calculate random value. "density" is from 1..-1.
971 werner 749
                if (item.density <= 1.) {
750
                    rand_val = mRandom->get();
751
                    if (item.density<0)
752
                        rand_val = 1. - rand_val;
753
                    rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
754
                } else {
755
                    // limited area: limit potential area using the "density" input parameter
756
                    rand_val = drandom() * qMin(item.density/100., 1.);
757
                }
732 werner 758
                ++total_tries;
549 werner 759
 
732 werner 760
                // key: rank of target pixel
761
                key = limit(int(pixel_list.count()*rand_val), 0, pixel_list.count()-1); // get from random number generator
549 werner 762
 
732 werner 763
                if (mInitHeightGrid) {
764
                    // calculate how good the selected pixel fits w.r.t. the predefined height
765
                    double p_value = pixel_list[key].h_max>0.?mHeightGridResponse->calculate(init_max_height/pixel_list[key].h_max):0.;
766
                    if (drandom() < p_value)
767
                        found = true;
768
                } else {
769
                    found = true;
770
                }
971 werner 771
                if (!last_locked_species && pixel_list[key].locked)
772
                    found = false;
732 werner 773
            }
774
            if (tries<0) ++total_misses;
549 werner 775
 
776
            // create a tree
777
            ResourceUnit *ru = pixel_list[key].resource_unit;
778
            int tree_idx = ru->newTreeIndex();
779
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
780
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
781
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
782
            tree.setSpecies(item.species);
783
            if (item.age<=0)
784
                tree.setAge(0,tree.height());
785
            else
786
                tree.setAge(item.age, tree.height());
787
            tree.setRU(ru);
788
            tree.setup();
789
            total_count++;
790
 
791
            // store in the multiHash the position of the pixel and the tree_idx in the resepctive resource unit
792
            tree_map.insert(pixel_list[key].pixelOffset, tree_idx);
793
            pixel_list[key].basal_area+=tree.basalArea(); // aggregate the basal area for each 10m pixel
971 werner 794
            if (last_locked_species)
795
                pixel_list[key].locked = true;
549 werner 796
 
797
            // resort list
971 werner 798
            if (last_locked_species==0 && ((total_count < 20 && i%2==0)
549 werner 799
                || (total_count<100 && i%10==0 )
971 werner 800
                || (i%30==0)) ) {
549 werner 801
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
802
            }
803
        }
804
    }
732 werner 805
    if (total_misses>0 || total_tries > total_count) {
744 werner 806
        if (logLevelInfo()) qDebug() << "init for stand" << stand_id << "treecount:" << total_count << ", tries:" << total_tries << ", misses:" << total_misses << ", %miss:" << qRound(total_misses*100 / (double)total_count);
732 werner 807
    }
549 werner 808
 
809
    int bits, index, pos;
810
    int c;
811
    QList<int> trees;
812
    QPoint tree_pos;
813
 
814
    foreach(const SInitPixel &p, pixel_list) {
815
        trees = tree_map.values(p.pixelOffset);
816
        c = trees.count();
817
        bits = 0;
818
        index = -1;
819
        double r;
820
        foreach(int tree_idx, trees) {
821
            if (c>18) {
822
                index = (index + 1)%25;
823
            } else {
824
                int stop=1000;
825
                index = 0;
826
                do {
827
                    // search a random position
828
                    r = drandom();
829
                    index = limit(int(25 *  r*r), 0, 24); // use rnd()^2 to search for locations -> higher number of low indices (i.e. 50% of lookups in first 25% of locations)
830
                } while (isBitSet(bits, index)==true && stop--);
831
                if (!stop)
832
                    qDebug() << "executeiLandInit: found no free bit.";
833
                setBit(bits, index, true); // mark position as used
834
            }
835
            // get position from fixed lists (one for even, one for uneven resource units)
836
            pos = p.resource_unit->index()%2?evenlist[index]:unevenlist[index];
837
            tree_pos = p.pixelOffset * cPxPerHeight; // convert to LIF index
600 werner 838
            tree_pos += QPoint(pos/cPxPerHeight, pos%cPxPerHeight);
550 werner 839
 
549 werner 840
            p.resource_unit->trees()[tree_idx].setPosition(tree_pos);
734 werner 841
            // test if tree position is valid..
842
            if (!GlobalSettings::instance()->model()->grid()->isIndexValid(tree_pos))
843
                qDebug() << "Standloader: invalid position!";
549 werner 844
        }
845
    }
971 werner 846
    if (logLevelInfo())
847
        qDebug() << "init for stand" << stand_id << "with area" << "area (m2)" << grid->area(stand_id) << "count of 10m pixels:"  << indices.count() << "initialized trees:" << total_count;
550 werner 848
 
549 werner 849
}
600 werner 850
 
851
/// a (hacky) way of adding saplings of a certain age to a stand defined by 'stand_id'.
852
int StandLoader::loadSaplings(const QString &content, int stand_id, const QString &fileName)
853
{
777 werner 854
    Q_UNUSED(fileName);
603 werner 855
    const MapGrid *stand_grid;
856
    if (mCurrentMap)
857
        stand_grid = mCurrentMap; // if set
858
    else
859
        stand_grid = GlobalSettings::instance()->model()->standGrid(); // default
860
 
600 werner 861
    QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels
862
    if (indices.isEmpty()) {
863
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
864
        return -1;
865
    }
866
    double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare)
867
 
868
    // parse the content of the init-file
869
    // species
870
    CSVFile init;
871
    init.loadFromString(content);
872
    int ispecies = init.columnIndex("species");
873
    int icount = init.columnIndex("count");
874
    int iheight = init.columnIndex("height");
875
    int iage = init.columnIndex("age");
876
    if (ispecies==-1 || icount==-1)
877
        throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!");
878
 
879
    const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet();
880
    double height, age;
881
    int total = 0;
882
    for (int row=0;row<init.rowCount();++row) {
883
        int pxcount = init.value(row, icount).toDouble() * area_factor + 0.5; // no. of pixels that should be filled (sapling grid is the same resolution as the lif-grid)
884
        const Species *species = set->species(init.value(row, ispecies).toString());
885
        if (!species)
886
            throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString()));
887
        height = iheight==-1?0.05: init.value(row, iheight).toDouble();
888
        age = iage==-1?1:init.value(row,iage).toDouble();
889
 
890
        int misses = 0;
891
        int hits = 0;
892
        while (hits < pxcount) {
893
           int rnd_index = irandom(0, indices.count()-1);
894
           QPoint offset=stand_grid->grid().indexOf(indices[rnd_index]);
895
           ResourceUnit *ru = GlobalSettings::instance()->model()->ru(stand_grid->grid().cellCenterPoint(offset));
896
           //
897
           offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
898
           int in_p = irandom(0, cPxPerHeight*cPxPerHeight-1); // index of lif-pixel
899
           offset += QPoint(in_p / cPxPerHeight, in_p % cPxPerHeight);
1042 werner 900
           if (!ru || ru->saplingHeightForInit(offset) > height) {
600 werner 901
               misses++;
902
           } else {
903
               // ok
904
               hits++;
951 werner 905
               ru->resourceUnitSpecies(species).changeSapling().addSapling(offset, height, age);
600 werner 906
           }
907
           if (misses > 3*pxcount) {
908
               qDebug() << "tried to add" << pxcount << "saplings at stand" << stand_id << "but failed in finding enough free positions. Added" << hits << "and stopped.";
909
               break;
910
           }
911
        }
912
        total += hits;
913
 
914
    }
915
    return total;
916
}
966 werner 917
 
967 werner 918
bool LIFValueHigher(const float *a, const float *b)
966 werner 919
{
967 werner 920
    return *a > *b;
966 werner 921
}
922
 
967 werner 923
int StandLoader::loadSaplingsLIF(int stand_id, const CSVFile &init, int low_index, int high_index)
966 werner 924
{
925
    const MapGrid *stand_grid;
926
    if (mCurrentMap)
927
        stand_grid = mCurrentMap; // if set
928
    else
929
        stand_grid = GlobalSettings::instance()->model()->standGrid(); // default
930
 
967 werner 931
    if (!stand_grid->isValid(stand_id))
932
        return 0;
933
 
966 werner 934
    QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels
935
    if (indices.isEmpty()) {
936
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
967 werner 937
        return 0;
966 werner 938
    }
939
    // prepare space for LIF-pointers (2m Pixel)
940
    QVector<float*> lif_ptrs;
967 werner 941
    lif_ptrs.reserve(indices.size() * cPxPerHeight * cPxPerHeight);
966 werner 942
    for (int l=0;l<indices.size();++l){
943
        QPoint offset=stand_grid->grid().indexOf(indices[l]);
944
        offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
967 werner 945
        for (int y=0;y<cPxPerHeight;++y)
946
            for(int x=0;x<cPxPerHeight;++x)
947
                lif_ptrs.push_back( GlobalSettings::instance()->model()->grid()->ptr(offset.x()+x, offset.y()+y) );
966 werner 948
    }
949
    // sort based on LIF-Value
967 werner 950
    std::sort(lif_ptrs.begin(), lif_ptrs.end(), LIFValueHigher); // higher: highest values first
966 werner 951
 
952
 
953
    double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare)
954
 
955
    // parse the content of the init-file
956
    // species
957
    int ispecies = init.columnIndex("species");
958
    int icount = init.columnIndex("count");
959
    int iheight = init.columnIndex("height");
960
    int iheightfrom = init.columnIndex("height_from");
961
    int iheightto = init.columnIndex("height_to");
962
    int iage = init.columnIndex("age");
963
    int itopage = init.columnIndex("age4m");
964
    int iminlif = init.columnIndex("min_lif");
965
    if ((iheightfrom==-1) ^ (iheightto==-1))
966
        throw IException("Error while loading saplings: height not correctly provided. Use either 'height' or 'height_from' and 'height_to'.");
967
    if (ispecies==-1 || icount==-1)
968
        throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!");
969
 
970
    const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet();
971
    double height, age;
972
    int total = 0;
967 werner 973
    for (int row=low_index;row<=high_index;++row) {
974
        double pxcount = init.value(row, icount).toDouble() * area_factor; // no. of pixels that should be filled (sapling grid is the same resolution as the lif-grid)
966 werner 975
        const Species *species = set->species(init.value(row, ispecies).toString());
976
        if (!species)
977
            throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString()));
978
        height = iheight==-1?0.05: init.value(row, iheight).toDouble();
979
        age = iage==-1?1:init.value(row,iage).toDouble();
980
        double age4m = itopage==-1?10:init.value(row, itopage).toDouble();
981
        double height_from = iheightfrom==-1?-1.: init.value(row, iheightfrom).toDouble();
982
        double height_to = iheightto==-1?-1.: init.value(row, iheightto).toDouble();
983
        double min_lif = iminlif==-1?1.: init.value(row, iminlif).toDouble();
984
        // find LIF-level in the pixels
985
        int min_lif_index = 0;
986
        for (QVector<float*>::ConstIterator it=lif_ptrs.constBegin(); it!=lif_ptrs.constEnd(); ++it, ++min_lif_index)
967 werner 987
            if (**it <= min_lif)
966 werner 988
                break;
989
 
990
        if (min_lif_index < pxcount) {
991
            // not enough LIF pixels available
992
            min_lif_index = pxcount; // try the brightest pixels (ie with the largest value for the LIF)
993
        }
994
 
967 werner 995
        double hits = 0.;
966 werner 996
        while (hits < pxcount) {
997
           int rnd_index = irandom(0, min_lif_index-1);
998
           if (iheightfrom!=-1) {
999
               height = limit(nrandom(height_from, height_to), 0.05,4.);
1000
               age = qMax(qRound(height/4. * age4m),1); // assume a linear relationship between height and age
1001
           }
1002
           QPoint offset = GlobalSettings::instance()->model()->grid()->indexOf(lif_ptrs[rnd_index]);
1003
 
1004
           ResourceUnit *ru = GlobalSettings::instance()->model()->ru(GlobalSettings::instance()->model()->grid()->cellCenterPoint(offset));
1042 werner 1005
           if (ru) {
1006
               ru->resourceUnitSpecies(species).changeSapling().addSapling(offset, height, age);
1007
               hits += ru->resourceUnitSpecies(species).sapling().representedStemNumber(height);
1008
           } else {
1009
               hits++; // avoid an infinite loop
1010
           }
967 werner 1011
 
1012
 
966 werner 1013
        }
1014
        total += pxcount;
1015
 
1016
    }
1017
    return total;
1018
}