Subversion Repositories public iLand

Rev

Rev 901 | Rev 911 | 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
 
210
void StandLoader::evaluateDebugTrees()
211
{
137 Werner 212
    // evaluate debugging
213
    QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
299 werner 214
    int counter=0;
137 Werner 215
    if (!dbg_str.isEmpty()) {
733 werner 216
        if (dbg_str == "debugstamp") {
734 werner 217
            qDebug() << "debug_tree = debugstamp: try touching all trees...";
733 werner 218
            // try to force an error if a stamp is invalid
219
            AllTreeIterator at(GlobalSettings::instance()->model());
753 werner 220
            double total_offset=0.;
733 werner 221
            while (Tree *t=at.next()) {
222
                total_offset += t->stamp()->offset();
734 werner 223
                if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
224
                    qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
733 werner 225
            }
734 werner 226
            qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
733 werner 227
            return;
228
        }
229
        TreeWrapper tw;
230
        Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
137 Werner 231
        AllTreeIterator at(GlobalSettings::instance()->model());
232
        double result;
233
        while (Tree *t = at.next()) {
140 Werner 234
            tw.setTree(t);
137 Werner 235
            result = dexp.execute();
299 werner 236
            if (result) {
137 Werner 237
                t->enableDebugging();
299 werner 238
                counter++;
239
            }
137 Werner 240
        }
732 werner 241
        qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
137 Werner 242
    }
106 Werner 243
}
244
 
904 werner 245
 
732 werner 246
/// load a single init file. Calls loadPicusFile() or loadiLandFile()
247
/// @param fileName file to load
248
/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
549 werner 249
int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
281 werner 250
{
284 werner 251
    QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
252
    if (!QFile::exists(pathFileName))
253
        throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
281 werner 254
 
549 werner 255
    if (type=="picus" || type=="single") {
256
        if (stand_id>0)
257
            throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
294 werner 258
        return loadPicusFile(pathFileName, ru);
549 werner 259
    }
384 werner 260
    if (type=="iland" || type=="distribution")
549 werner 261
        return loadiLandFile(pathFileName, ru, stand_id);
281 werner 262
 
384 werner 263
    throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
281 werner 264
}
265
 
393 werner 266
int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
106 Werner 267
{
294 werner 268
    QString content = Helper::loadTextFile(fileName);
269
    if (content.isEmpty()) {
270
        qDebug() << "file not found: " + fileName;
393 werner 271
        return 0;
294 werner 272
    }
393 werner 273
    return loadSingleTreeList(content, ru, fileName);
294 werner 274
}
275
 
389 werner 276
/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.
277
  returns the number of loaded trees.
278
  */
279
int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
294 werner 280
{
106 Werner 281
    if (!ru)
282
        ru = mModel->ru();
284 werner 283
    Q_ASSERT(ru!=0);
106 Werner 284
 
284 werner 285
    QPointF offset = ru->boundingBox().topLeft();
286
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
287
 
384 werner 288
    QString my_content(content);
289
    // cut out the <trees> </trees> part if present
290
    if (content.contains("<trees>")) {
291
        QRegExp rx(".*<trees>(.*)</trees>.*");
292
        rx.indexIn(content, 0);
293
        if (rx.capturedTexts().count()<1)
389 werner 294
            return 0;
384 werner 295
        my_content = rx.cap(1).trimmed();
296
    }
297
 
294 werner 298
    QStringList lines=my_content.split('\n');
106 Werner 299
    if (lines.count()<2)
389 werner 300
        return 0;
384 werner 301
    // drop comments
302
    while (!lines.isEmpty() && lines.front().startsWith('#') )
303
        lines.pop_front();
304
    while (!lines.isEmpty() && lines.last().isEmpty())
305
        lines.removeLast();
306
 
106 Werner 307
    char sep='\t';
308
    if (!lines[0].contains(sep))
309
        sep=';';
155 werner 310
    QStringList headers = lines[0].trimmed().split(sep);
138 Werner 311
 
247 werner 312
    int iID = headers.indexOf("id");
106 Werner 313
    int iX = headers.indexOf("x");
314
    int iY = headers.indexOf("y");
315
    int iBhd = headers.indexOf("bhdfrom");
384 werner 316
    if (iBhd<0)
317
        iBhd = headers.indexOf("dbh");
318
    double height_conversion = 100.;
106 Werner 319
    int iHeight = headers.indexOf("treeheight");
384 werner 320
    if (iHeight<0) {
321
        iHeight = headers.indexOf("height");
322
        height_conversion = 1.; // in meter
323
    }
106 Werner 324
    int iSpecies = headers.indexOf("species");
170 werner 325
    int iAge = headers.indexOf("age");
204 werner 326
    if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1)
384 werner 327
        throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
106 Werner 328
 
138 Werner 329
    double dbh;
384 werner 330
    bool ok;
389 werner 331
    int cnt=0;
384 werner 332
    QString speciesid;
106 Werner 333
    for (int i=1;i<lines.count();i++) {
334
        QString &line = lines[i];
138 Werner 335
        dbh = line.section(sep, iBhd, iBhd).toDouble();
336
        if (dbh<5.)
337
            continue;
285 werner 338
 
106 Werner 339
        QPointF f;
340
        if (iX>=0 && iY>=0) {
341
           f.setX( line.section(sep, iX, iX).toDouble() );
342
           f.setY( line.section(sep, iY, iY).toDouble() );
343
           f+=offset;
285 werner 344
 
106 Werner 345
        }
285 werner 346
        // position valid?
347
        if (!mModel->heightGrid()->valueAt(f).isValid())
348
            continue;
349
        Tree &tree = ru->newTree();
350
        tree.setPosition(f);
247 werner 351
        if (iID>=0)
352
            tree.setId(line.section(sep, iID, iID).toInt() );
106 Werner 353
 
138 Werner 354
        tree.setDbh(dbh);
384 werner 355
        tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
356
 
389 werner 357
        speciesid = line.section(sep, iSpecies, iSpecies).trimmed();
270 werner 358
        int picusid = speciesid.toInt(&ok);
359
        if (ok) {
384 werner 360
            int idx = picusSpeciesIds.indexOf(picusid);
270 werner 361
            if (idx==-1)
362
                throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
138 Werner 363
            speciesid = iLandSpeciesIds[idx];
270 werner 364
        }
138 Werner 365
        Species *s = speciesSet->species(speciesid);
366
        if (!ru || !s)
389 werner 367
            throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
388 werner 368
        tree.setSpecies(s);
107 Werner 369
 
388 werner 370
        ok = true;
371
        if (iAge>=0)
372
           tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a *real* age
373
        if (iAge<0 || !ok || tree.age()==0)
374
           tree.setAge(0, tree.height()); // no real tree age available
375
 
138 Werner 376
        tree.setRU(ru);
106 Werner 377
        tree.setup();
389 werner 378
        cnt++;
106 Werner 379
    }
389 werner 380
    return cnt;
106 Werner 381
    //qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
382
    //qDebug() << "lines: " << lines;
383
}
287 werner 384
 
393 werner 385
/** initialize trees on a resource unit based on dbh distributions.
386
  use a fairly clever algorithm to determine tree positions.
387
  see http://iland.boku.ac.at/initialize+trees
388
  @param content tree init file (including headers) in a string
389
  @param ru resource unit
390
  @param fileName source file name (for error reporting)
391
  @return number of trees added
392
  */
549 werner 393
int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
287 werner 394
{
904 werner 395
    int total_count = parseInitFile(content, fileName, ru);
396
    if (total_count==0)
397
        return 0;
398
 
399
 
400
    // setup the random distribution
401
    QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
402
    if (logLevelInfo())  qDebug() << "density function:" << density_func;
403
    if (!mRandom || (mRandom->densityFunction()!= density_func)) {
404
        if (mRandom)
405
            delete mRandom;
406
        mRandom=new RandomCustomPDF(density_func);
407
        if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
408
    }
409
    if (stand_id>0) {
410
        // execute stand based initialization
411
        executeiLandInitStand(stand_id);
412
    } else {
413
        // exeucte the initialization based on single resource units
414
        executeiLandInit(ru);
415
        ru->cleanTreeList();
416
    }
417
    return total_count;
418
 
419
}
420
 
421
int StandLoader::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru)
422
{
287 werner 423
    if (!ru)
424
        ru = mModel->ru();
425
    Q_ASSERT(ru!=0);
426
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
427
    Q_ASSERT(speciesSet!=0);
428
 
431 werner 429
    //DebugTimer t("StandLoader::loadiLandFile");
294 werner 430
    CSVFile infile;
431
    infile.loadFromString(content);
432
 
287 werner 433
    int icount = infile.columnIndex("count");
434
    int ispecies = infile.columnIndex("species");
435
    int idbh_from = infile.columnIndex("dbh_from");
436
    int idbh_to = infile.columnIndex("dbh_to");
437
    int ihd = infile.columnIndex("hd");
438
    int iage = infile.columnIndex("age");
311 werner 439
    int idensity = infile.columnIndex("density");
287 werner 440
    if (icount<0 || ispecies<0 || idbh_from<0 || idbh_to<0 || ihd<0 || iage<0)
441
        throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));
442
 
904 werner 443
    int istandid = infile.columnIndex("stand_id");
294 werner 444
    mInitItems.clear();
904 werner 445
    mStandInitItems.clear();
446
 
287 werner 447
    InitFileItem item;
384 werner 448
    bool ok;
393 werner 449
    int total_count = 0;
287 werner 450
    for (int row=0;row<infile.rowCount();row++) {
549 werner 451
        item.count = infile.value(row, icount).toInt();
452
        total_count += item.count;
453
        item.dbh_from = infile.value(row, idbh_from).toDouble();
454
        item.dbh_to = infile.value(row, idbh_to).toDouble();
455
        item.hd = infile.value(row, ihd).toDouble();
775 werner 456
        if (item.hd==0. || item.dbh_from / 100. * item.hd < 4.)
782 werner 457
            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) ;
458
            //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 459
        ok = true;
460
        if (iage>=0)
461
            item.age = infile.value(row, iage).toInt(&ok);
462
        if (iage<0 || !ok)
463
            item.age = 0;
384 werner 464
 
549 werner 465
        item.species = speciesSet->species(infile.value(row, ispecies).toString());
466
        if (idensity>=0)
467
            item.density = infile.value(row, idensity).toDouble();
468
        else
469
            item.density = 0.;
470
        if (item.density<-1 || item.density>1)
471
            throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
472
                             .arg(item.density)
473
                             .arg(fileName)
474
                             .arg(row));
475
        if (!item.species) {
476
            throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
477
                             .arg(infile.value(row, ispecies).toString())
478
                             .arg(fileName)
479
                             .arg(row));
480
        }
904 werner 481
        if (istandid>=0) {
482
            int standid = infile.value(row,istandid).toInt();
483
            mStandInitItems[standid].push_back(item);
484
        } else {
485
            mInitItems.push_back(item);
486
        }
287 werner 487
    }
393 werner 488
    return total_count;
287 werner 489
 
490
}
491
 
904 werner 492
 
549 werner 493
int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)
294 werner 494
{
495
    if (!QFile::exists(fileName))
496
        throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
497
    QString content = Helper::loadTextFile(fileName);
549 werner 498
    return loadDistributionList(content, ru, stand_id, fileName);
294 werner 499
}
500
 
287 werner 501
// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
502
// e.g. 12 = centerpixel, 0: upper left corner, ...
503
int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
504
                     1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
290 werner 505
int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
506
                       6,18,2,10,4,24,12,0,8,14,20,22};
507
 
508
 
549 werner 509
 
287 werner 510
// sort function
511
bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
311 werner 512
{
513
    return s1.second < s2.second;
514
}
515
 
549 werner 516
struct SInitPixel {
732 werner 517
    double basal_area; // accumulated basal area
518
    QPoint pixelOffset; // location of the pixel
519
    ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to
520
    double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)
521
    SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}
549 werner 522
};
523
 
524
bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)
525
{
526
    return s1.basal_area < s2.basal_area;
527
}
528
 
529
 
311 werner 530
/**
531
*/
549 werner 532
 
290 werner 533
void StandLoader::executeiLandInit(ResourceUnit *ru)
287 werner 534
{
290 werner 535
 
287 werner 536
    QPointF offset = ru->boundingBox().topLeft();
537
    QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);
538
 
539
    // a multimap holds a list for all trees.
540
    // key is the index of a 10x10m pixel within the resource unit
541
    QMultiMap<int, int> tree_map;
549 werner 542
    //QHash<int,SInitPixel> tcount;
543
 
287 werner 544
    QVector<QPair<int, double> > tcount; // counts
545
    for (int i=0;i<100;i++)
546
        tcount.push_back(QPair<int,double>(i,0.));
547
 
548
    int key;
311 werner 549
    double rand_val, rand_fraction;
312 werner 550
    int total_count = 0;
290 werner 551
    foreach(const InitFileItem &item, mInitItems) {
311 werner 552
        rand_fraction = fabs(double(item.density));
287 werner 553
        for (int i=0;i<item.count;i++) {
554
            // create trees
555
            int tree_idx = ru->newTreeIndex();
556
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
707 werner 557
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
287 werner 558
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
559
            tree.setSpecies(item.species);
384 werner 560
            if (item.age<=0)
388 werner 561
                tree.setAge(0,tree.height());
381 werner 562
            else
388 werner 563
                tree.setAge(item.age, tree.height());
287 werner 564
            tree.setRU(ru);
565
            tree.setup();
312 werner 566
            total_count++;
311 werner 567
 
568
            // calculate random value. "density" is from 1..-1.
569
            rand_val = mRandom->get();
570
            if (item.density<0)
571
                rand_val = 1. - rand_val;
705 werner 572
            rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
311 werner 573
 
287 werner 574
            // key: rank of target pixel
575
            // first: index of target pixel
576
            // second: sum of target pixel
311 werner 577
            key = limit(int(100*rand_val), 0, 99); // get from random number generator
287 werner 578
            tree_map.insert(tcount[key].first, tree_idx); // store tree in map
579
            tcount[key].second+=tree.basalArea(); // aggregate the basal area for each 10m pixel
312 werner 580
            if ( (total_count < 20 && i%2==0)
581
                || (total_count<100 && i%10==0 )
582
                || (i%30==0) ) {
287 werner 583
                qSort(tcount.begin(), tcount.end(), sortPairLessThan);
312 werner 584
            }
287 werner 585
        }
586
        qSort(tcount.begin(), tcount.end(), sortPairLessThan);
587
    }
588
 
589
    int bits, index, pos;
590
    int c;
591
    QList<int> trees;
592
    QPoint tree_pos;
289 werner 593
 
287 werner 594
    for (int i=0;i<100;i++) {
595
        trees = tree_map.values(i);
596
        c = trees.count();
290 werner 597
        QPointF pixel_center = ru->boundingBox().topLeft() + QPointF((i/10)*10. + 5., (i%10)*10. + 5.);
598
        if (!mModel->heightGrid()->valueAt(pixel_center).isValid()) {
599
            // no trees on that pixel: let trees die
600
            foreach(int tree_idx, trees) {
601
                ru->trees()[tree_idx].die();
602
            }
603
            continue;
604
        }
605
 
287 werner 606
        bits = 0;
607
        index = -1;
290 werner 608
        double r;
287 werner 609
        foreach(int tree_idx, trees) {
610
            if (c>18) {
611
                index = (index + 1)%25;
612
            } else {
613
                int stop=1000;
290 werner 614
                index = 0;
287 werner 615
                do {
290 werner 616
                    //r = drandom();
617
                    //if (r<0.5)  // skip position with a prob. of 50% -> adds a little "noise"
618
                    //    index++;
619
                    //index = (index + 1)%25; // increase and roll over
620
 
621
                    // search a random position
622
                    r = drandom();
623
                    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 624
                } while (isBitSet(bits, index)==true && stop--);
625
                if (!stop)
626
                    qDebug() << "executeiLandInit: found no free bit.";
627
                setBit(bits, index, true); // mark position as used
628
            }
290 werner 629
            // get position from fixed lists (one for even, one for uneven resource units)
630
            pos = ru->index()%2?evenlist[index]:unevenlist[index];
287 werner 631
            tree_pos = offsetIdx  // position of resource unit
632
                       + QPoint(5*(i/10), 5*(i%10)) // relative position of 10x10m pixel
633
                       + QPoint(pos/5, pos%5); // relative position within 10x10m pixel
634
            //qDebug() << tree_no++ << "to" << index;
635
            ru->trees()[tree_idx].setPosition(tree_pos);
636
        }
637
    }
638
}
549 werner 639
 
640
// provide a hashing function for the QPoint type (needed from stand init function below)
641
inline uint qHash(const QPoint &key)
642
 {
643
     return qHash(key.x()) ^ qHash(key.y());
644
 }
645
 
550 werner 646
 
647
// Initialization routine based on a stand map.
648
// Basically a list of 10m pixels for a given stand is retrieved
649
// and the filled with the same procedure as the resource unit based init
650
// see http://iland.boku.ac.at/initialize+trees
549 werner 651
void StandLoader::executeiLandInitStand(int stand_id)
652
{
653
 
654
    const MapGrid *grid = GlobalSettings::instance()->model()->standGrid();
901 werner 655
    if (mCurrentMap)
656
        grid = mCurrentMap;
549 werner 657
 
732 werner 658
    // get a list of positions of all pixels that belong to our stand
549 werner 659
    QList<int> indices = grid->gridIndices(stand_id);
550 werner 660
    if (indices.isEmpty()) {
661
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
662
        return;
549 werner 663
    }
664
    // a multiHash holds a list for all trees.
665
    // key is the location of the 10x10m pixel
666
    QMultiHash<QPoint, int> tree_map;
667
    QList<SInitPixel> pixel_list; // working list of all 10m pixels
732 werner 668
    pixel_list.reserve(indices.size());
549 werner 669
 
670
    foreach (int i, indices) {
671
       SInitPixel p;
672
       p.pixelOffset = grid->grid().indexOf(i); // index in the 10m grid
732 werner 673
       p.resource_unit = GlobalSettings::instance()->model()->ru( grid->grid().cellCenterPoint(p.pixelOffset));
674
       if (mInitHeightGrid)
675
           p.h_max = mInitHeightGrid->grid().constValueAtIndex(p.pixelOffset);
549 werner 676
       pixel_list.append(p);
677
    }
678
    double area_factor = grid->area(stand_id) / 10000.;
679
 
732 werner 680
    int key=0;
549 werner 681
    double rand_val, rand_fraction;
682
    int total_count = 0;
732 werner 683
    int total_tries = 0;
684
    int total_misses = 0;
685
    if (mInitHeightGrid && !mHeightGridResponse)
686
        throw IException("executeiLandInitStand: trying to initialize with height grid but without response function.");
549 werner 687
    foreach(const InitFileItem &item, mInitItems) {
688
        rand_fraction = fabs(double(item.density));
689
        int count = item.count * area_factor + 0.5; // round
732 werner 690
        double init_max_height = item.dbh_to/100. * item.hd;
549 werner 691
        for (int i=0;i<count;i++) {
692
 
732 werner 693
            bool found = false;
694
            int tries = mHeightGridTries;
695
            while (!found &&--tries) {
696
                // calculate random value. "density" is from 1..-1.
697
                rand_val = mRandom->get();
698
                if (item.density<0)
699
                    rand_val = 1. - rand_val;
700
                rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
701
                ++total_tries;
549 werner 702
 
732 werner 703
                // key: rank of target pixel
704
                key = limit(int(pixel_list.count()*rand_val), 0, pixel_list.count()-1); // get from random number generator
549 werner 705
 
732 werner 706
                if (mInitHeightGrid) {
707
                    // calculate how good the selected pixel fits w.r.t. the predefined height
708
                    double p_value = pixel_list[key].h_max>0.?mHeightGridResponse->calculate(init_max_height/pixel_list[key].h_max):0.;
709
                    if (drandom() < p_value)
710
                        found = true;
711
                } else {
712
                    found = true;
713
                }
714
            }
715
            if (tries<0) ++total_misses;
549 werner 716
 
717
            // create a tree
718
            ResourceUnit *ru = pixel_list[key].resource_unit;
719
            int tree_idx = ru->newTreeIndex();
720
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
721
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
722
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
723
            tree.setSpecies(item.species);
724
            if (item.age<=0)
725
                tree.setAge(0,tree.height());
726
            else
727
                tree.setAge(item.age, tree.height());
728
            tree.setRU(ru);
729
            tree.setup();
730
            total_count++;
731
 
732
            // store in the multiHash the position of the pixel and the tree_idx in the resepctive resource unit
733
            tree_map.insert(pixel_list[key].pixelOffset, tree_idx);
734
            pixel_list[key].basal_area+=tree.basalArea(); // aggregate the basal area for each 10m pixel
735
 
736
            // resort list
737
            if ( (total_count < 20 && i%2==0)
738
                || (total_count<100 && i%10==0 )
739
                || (i%30==0) ) {
740
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
741
            }
742
        }
743
        qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
744
    }
732 werner 745
    if (total_misses>0 || total_tries > total_count) {
744 werner 746
        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 747
    }
549 werner 748
 
749
    int bits, index, pos;
750
    int c;
751
    QList<int> trees;
752
    QPoint tree_pos;
753
 
754
    foreach(const SInitPixel &p, pixel_list) {
755
        trees = tree_map.values(p.pixelOffset);
756
        c = trees.count();
757
        bits = 0;
758
        index = -1;
759
        double r;
760
        foreach(int tree_idx, trees) {
761
            if (c>18) {
762
                index = (index + 1)%25;
763
            } else {
764
                int stop=1000;
765
                index = 0;
766
                do {
767
                    // search a random position
768
                    r = drandom();
769
                    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)
770
                } while (isBitSet(bits, index)==true && stop--);
771
                if (!stop)
772
                    qDebug() << "executeiLandInit: found no free bit.";
773
                setBit(bits, index, true); // mark position as used
774
            }
775
            // get position from fixed lists (one for even, one for uneven resource units)
776
            pos = p.resource_unit->index()%2?evenlist[index]:unevenlist[index];
777
            tree_pos = p.pixelOffset * cPxPerHeight; // convert to LIF index
600 werner 778
            tree_pos += QPoint(pos/cPxPerHeight, pos%cPxPerHeight);
550 werner 779
 
549 werner 780
            p.resource_unit->trees()[tree_idx].setPosition(tree_pos);
734 werner 781
            // test if tree position is valid..
782
            if (!GlobalSettings::instance()->model()->grid()->isIndexValid(tree_pos))
783
                qDebug() << "Standloader: invalid position!";
549 werner 784
        }
785
    }
550 werner 786
    if (logLevelInfo()) qDebug() << "init for stand" << stand_id << "with area" << "area (m2)" << grid->area(stand_id) << "count of 10m pixels:"  << indices.count() << "initialized trees:" << total_count;
787
 
549 werner 788
}
600 werner 789
 
790
/// a (hacky) way of adding saplings of a certain age to a stand defined by 'stand_id'.
791
int StandLoader::loadSaplings(const QString &content, int stand_id, const QString &fileName)
792
{
777 werner 793
    Q_UNUSED(fileName);
603 werner 794
    const MapGrid *stand_grid;
795
    if (mCurrentMap)
796
        stand_grid = mCurrentMap; // if set
797
    else
798
        stand_grid = GlobalSettings::instance()->model()->standGrid(); // default
799
 
600 werner 800
    QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels
801
    if (indices.isEmpty()) {
802
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
803
        return -1;
804
    }
805
    double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare)
806
 
807
    // parse the content of the init-file
808
    // species
809
    CSVFile init;
810
    init.loadFromString(content);
811
    int ispecies = init.columnIndex("species");
812
    int icount = init.columnIndex("count");
813
    int iheight = init.columnIndex("height");
814
    int iage = init.columnIndex("age");
815
    if (ispecies==-1 || icount==-1)
816
        throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!");
817
 
818
    const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet();
819
    double height, age;
820
    int total = 0;
821
    for (int row=0;row<init.rowCount();++row) {
822
        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)
823
        const Species *species = set->species(init.value(row, ispecies).toString());
824
        if (!species)
825
            throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString()));
826
        height = iheight==-1?0.05: init.value(row, iheight).toDouble();
827
        age = iage==-1?1:init.value(row,iage).toDouble();
828
 
829
        int misses = 0;
830
        int hits = 0;
831
        while (hits < pxcount) {
832
           int rnd_index = irandom(0, indices.count()-1);
833
           QPoint offset=stand_grid->grid().indexOf(indices[rnd_index]);
834
           ResourceUnit *ru = GlobalSettings::instance()->model()->ru(stand_grid->grid().cellCenterPoint(offset));
835
           //
836
           offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
837
           int in_p = irandom(0, cPxPerHeight*cPxPerHeight-1); // index of lif-pixel
838
           offset += QPoint(in_p / cPxPerHeight, in_p % cPxPerHeight);
839
           if (ru->saplingHeightForInit(offset) > height) {
840
               misses++;
841
           } else {
842
               // ok
843
               hits++;
844
               ru->resourceUnitSpecies(species).changeSapling().addSapling(offset);
845
           }
846
           if (misses > 3*pxcount) {
847
               qDebug() << "tried to add" << pxcount << "saplings at stand" << stand_id << "but failed in finding enough free positions. Added" << hits << "and stopped.";
848
               break;
849
           }
850
        }
851
        total += hits;
852
 
853
    }
854
    return total;
855
}