Subversion Repositories public iLand

Rev

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