Subversion Repositories public iLand

Rev

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