Subversion Repositories public iLand

Rev

Rev 384 | Rev 389 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
1
 
106 Werner 2
#include "global.h"
3
#include "standloader.h"
4
 
290 werner 5
 
106 Werner 6
#include "grid.h"
7
#include "model.h"
189 iland 8
#include "resourceunit.h"
106 Werner 9
#include "speciesset.h"
10
 
11
#include "helper.h"
290 werner 12
#include "random.h"
137 Werner 13
#include "expression.h"
140 Werner 14
#include "expressionwrapper.h"
284 werner 15
#include "environment.h"
287 werner 16
#include "csvfile.h"
106 Werner 17
 
381 werner 18
// provide a mapping between "Picus"-style and "iLand"-style species Ids
384 werner 19
QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
155 werner 20
QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";
290 werner 21
StandLoader::~StandLoader()
22
{
23
    if (mRandom)
24
        delete mRandom;
25
}
106 Werner 26
 
281 werner 27
void StandLoader::loadForUnit()
28
{
106 Werner 29
 
281 werner 30
}
106 Werner 31
 
281 werner 32
void StandLoader::copyTrees()
33
{
34
    // we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards
35
    const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();
36
    ResourceUnit **p = ruGrid.begin();
37
    if (!p)
38
        throw IException("Standloader: invalid resource unit pointer!");
39
    ++p; // skip the first...
40
    const QVector<Tree> &tocopy = mModel->ru()->trees();
41
    for (; p!=ruGrid.end(); ++p) {
42
        QRectF rect = (*p)->boundingBox();
43
        foreach(const Tree& tree, tocopy) {
44
            Tree &newtree = (*p)->newTree();
45
            newtree = tree; // copy tree data...
46
            newtree.setPosition(tree.position()+(*p)->boundingBox().topLeft());
47
            newtree.setRU(*p);
48
            newtree.setNewId();
49
        }
50
    }
51
    qDebug() << Tree::statCreated() << "trees loaded / copied.";
52
}
53
 
284 werner 54
/** main routine of the stand setup.
55
*/
106 Werner 56
void StandLoader::processInit()
57
{
58
    GlobalSettings *g = GlobalSettings::instance();
192 werner 59
    XmlHelper xml(g->settings().node("model.initialization"));
191 werner 60
 
281 werner 61
    QString copy_mode = xml.value("mode", "copy");
62
    QString type = xml.value("type", "");
192 werner 63
    QString  fileName = xml.value("file", "");
106 Werner 64
 
65
    Tree::resetStatistics();
281 werner 66
 
67
    // one global init-file for the whole area:
284 werner 68
    if (copy_mode=="single") {
281 werner 69
        loadInitFile(fileName, type);
299 werner 70
        evaluateDebugTrees();
281 werner 71
        return;
284 werner 72
    }
281 werner 73
 
74
    // copy trees from first unit to all other units:
75
    if (copy_mode=="copy") {
76
        loadInitFile(fileName, type);
77
        copyTrees();
299 werner 78
        evaluateDebugTrees();
281 werner 79
        return;
106 Werner 80
    }
137 Werner 81
 
281 werner 82
    // call a single tree init for each resource unit
83
    if (copy_mode=="unit") {
284 werner 84
        foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
85
            ResourceUnit *ru = const_cast<ResourceUnit*>(const_ru);
86
            // set environment
87
            g->model()->environment()->setPosition(ru->boundingBox().center());
88
            type = xml.value("type", "");
89
            fileName = xml.value("file", "");
319 werner 90
            if (fileName.isEmpty())
91
                continue;
284 werner 92
            loadInitFile(fileName, type, ru);
93
            qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
94
        }
299 werner 95
        evaluateDebugTrees();
284 werner 96
        return;
281 werner 97
    }
299 werner 98
 
281 werner 99
    throw IException("StandLoader::processInit: invalid initalization.mode!");
100
}
101
 
102
void StandLoader::evaluateDebugTrees()
103
{
137 Werner 104
    // evaluate debugging
105
    QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
299 werner 106
    int counter=0;
137 Werner 107
    if (!dbg_str.isEmpty()) {
140 Werner 108
       TreeWrapper tw;
109
       Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
137 Werner 110
        AllTreeIterator at(GlobalSettings::instance()->model());
111
        double result;
112
        while (Tree *t = at.next()) {
140 Werner 113
            tw.setTree(t);
137 Werner 114
            result = dexp.execute();
299 werner 115
            if (result) {
137 Werner 116
                t->enableDebugging();
299 werner 117
                counter++;
118
            }
137 Werner 119
        }
120
    }
299 werner 121
    qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
106 Werner 122
}
123
 
284 werner 124
void StandLoader::loadInitFile(const QString &fileName, const QString &type, ResourceUnit *ru)
281 werner 125
{
284 werner 126
    QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
127
    if (!QFile::exists(pathFileName))
128
        throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
281 werner 129
 
384 werner 130
    if (type=="picus" || type=="single")
294 werner 131
        return loadPicusFile(pathFileName, ru);
384 werner 132
    if (type=="iland" || type=="distribution")
287 werner 133
        return loadiLandFile(pathFileName, ru);
281 werner 134
 
384 werner 135
    throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
281 werner 136
}
137
 
294 werner 138
void StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
106 Werner 139
{
294 werner 140
    QString content = Helper::loadTextFile(fileName);
141
    if (content.isEmpty()) {
142
        qDebug() << "file not found: " + fileName;
143
        return;
144
    }
145
    loadSingleTreeList(content, ru, fileName);
146
}
147
 
148
void StandLoader::loadSingleTreeList(const QString &content, ResourceUnit*ru, const QString &fileName)
149
{
106 Werner 150
    if (!ru)
151
        ru = mModel->ru();
284 werner 152
    Q_ASSERT(ru!=0);
106 Werner 153
 
284 werner 154
    QPointF offset = ru->boundingBox().topLeft();
155
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
156
 
384 werner 157
    QString my_content(content);
158
    // cut out the <trees> </trees> part if present
159
    if (content.contains("<trees>")) {
160
        QRegExp rx(".*<trees>(.*)</trees>.*");
161
        rx.indexIn(content, 0);
162
        if (rx.capturedTexts().count()<1)
163
            return;
164
        my_content = rx.cap(1).trimmed();
165
    }
166
 
294 werner 167
    QStringList lines=my_content.split('\n');
106 Werner 168
    if (lines.count()<2)
169
        return;
384 werner 170
    // drop comments
171
    while (!lines.isEmpty() && lines.front().startsWith('#') )
172
        lines.pop_front();
173
    while (!lines.isEmpty() && lines.last().isEmpty())
174
        lines.removeLast();
175
 
106 Werner 176
    char sep='\t';
177
    if (!lines[0].contains(sep))
178
        sep=';';
155 werner 179
    QStringList headers = lines[0].trimmed().split(sep);
138 Werner 180
 
247 werner 181
    int iID = headers.indexOf("id");
106 Werner 182
    int iX = headers.indexOf("x");
183
    int iY = headers.indexOf("y");
184
    int iBhd = headers.indexOf("bhdfrom");
384 werner 185
    if (iBhd<0)
186
        iBhd = headers.indexOf("dbh");
187
    double height_conversion = 100.;
106 Werner 188
    int iHeight = headers.indexOf("treeheight");
384 werner 189
    if (iHeight<0) {
190
        iHeight = headers.indexOf("height");
191
        height_conversion = 1.; // in meter
192
    }
106 Werner 193
    int iSpecies = headers.indexOf("species");
170 werner 194
    int iAge = headers.indexOf("age");
204 werner 195
    if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1)
384 werner 196
        throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
106 Werner 197
 
138 Werner 198
    double dbh;
384 werner 199
    bool ok;
200
    QString speciesid;
106 Werner 201
    for (int i=1;i<lines.count();i++) {
202
        QString &line = lines[i];
138 Werner 203
        dbh = line.section(sep, iBhd, iBhd).toDouble();
204
        if (dbh<5.)
205
            continue;
285 werner 206
 
106 Werner 207
        QPointF f;
208
        if (iX>=0 && iY>=0) {
209
           f.setX( line.section(sep, iX, iX).toDouble() );
210
           f.setY( line.section(sep, iY, iY).toDouble() );
211
           f+=offset;
285 werner 212
 
106 Werner 213
        }
285 werner 214
        // position valid?
215
        if (!mModel->heightGrid()->valueAt(f).isValid())
216
            continue;
217
        Tree &tree = ru->newTree();
218
        tree.setPosition(f);
247 werner 219
        if (iID>=0)
220
            tree.setId(line.section(sep, iID, iID).toInt() );
106 Werner 221
 
138 Werner 222
        tree.setDbh(dbh);
384 werner 223
        tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
224
 
225
        speciesid = line.section(sep, iSpecies, iSpecies);
270 werner 226
        int picusid = speciesid.toInt(&ok);
227
        if (ok) {
384 werner 228
            int idx = picusSpeciesIds.indexOf(picusid);
270 werner 229
            if (idx==-1)
230
                throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
138 Werner 231
            speciesid = iLandSpeciesIds[idx];
270 werner 232
        }
138 Werner 233
        Species *s = speciesSet->species(speciesid);
234
        if (!ru || !s)
270 werner 235
            throw IException(QString("Loading init-file: either ressource unit or species invalid. Species: %1").arg(speciesid));
388 werner 236
        tree.setSpecies(s);
107 Werner 237
 
388 werner 238
        ok = true;
239
        if (iAge>=0)
240
           tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a *real* age
241
        if (iAge<0 || !ok || tree.age()==0)
242
           tree.setAge(0, tree.height()); // no real tree age available
243
 
138 Werner 244
        tree.setRU(ru);
106 Werner 245
        tree.setup();
246
    }
247
    //qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
248
    //qDebug() << "lines: " << lines;
249
}
287 werner 250
 
294 werner 251
void StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, const QString &fileName)
287 werner 252
{
253
    if (!ru)
254
        ru = mModel->ru();
255
    Q_ASSERT(ru!=0);
256
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
257
    Q_ASSERT(speciesSet!=0);
258
 
259
    DebugTimer t("StandLoader::loadiLandFile");
294 werner 260
    CSVFile infile;
261
    infile.loadFromString(content);
262
 
287 werner 263
    int icount = infile.columnIndex("count");
264
    int ispecies = infile.columnIndex("species");
265
    int idbh_from = infile.columnIndex("dbh_from");
266
    int idbh_to = infile.columnIndex("dbh_to");
267
    int ihd = infile.columnIndex("hd");
268
    int iage = infile.columnIndex("age");
311 werner 269
    int idensity = infile.columnIndex("density");
287 werner 270
    if (icount<0 || ispecies<0 || idbh_from<0 || idbh_to<0 || ihd<0 || iage<0)
271
        throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));
272
 
294 werner 273
    mInitItems.clear();
287 werner 274
    InitFileItem item;
384 werner 275
    bool ok;
287 werner 276
    for (int row=0;row<infile.rowCount();row++) {
277
         item.count = infile.value(row, icount).toInt();
278
         item.dbh_from = infile.value(row, idbh_from).toDouble();
279
         item.dbh_to = infile.value(row, idbh_to).toDouble();
280
         item.hd = infile.value(row, ihd).toDouble();
384 werner 281
         ok = true;
381 werner 282
         if (iage>=0)
384 werner 283
             item.age = infile.value(row, iage).toInt(&ok);
284
         if (iage<0 || !ok)
285
             item.age = 0;
286
 
287 werner 287
         item.species = speciesSet->species(infile.value(row, ispecies).toString());
311 werner 288
         if (idensity>=0)
289
             item.density = infile.value(row, idensity).toDouble();
290
         else
291
             item.density = 0.;
292
         if (item.density<-1 || item.density>1)
293
             throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
294
                              .arg(item.density)
295
                              .arg(fileName)
296
                              .arg(row));
287 werner 297
         if (!item.species) {
298
             throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
299
                              .arg(infile.value(row, ispecies).toString())
300
                              .arg(fileName)
301
                              .arg(row));
302
         }
290 werner 303
         mInitItems.push_back(item);
287 werner 304
    }
290 werner 305
    // setup the random distribution
306
    QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
307
    qDebug() << "density function:" << density_func;
308
    if (!mRandom || (mRandom->densityFunction()!= density_func)) {
309
        if (mRandom)
310
            delete mRandom;
311
        mRandom=new RandomCustomPDF(density_func);
312
        qDebug() << "new probabilty density function:" << density_func;
313
    }
314
 
287 werner 315
    // exeucte the
316
    executeiLandInit(ru);
290 werner 317
    ru->cleanTreeList();
287 werner 318
 
319
}
320
 
294 werner 321
void StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru)
322
{
323
    if (!QFile::exists(fileName))
324
        throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
325
    QString content = Helper::loadTextFile(fileName);
326
    loadDistributionList(content, ru, fileName);
327
}
328
 
287 werner 329
// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
330
// e.g. 12 = centerpixel, 0: upper left corner, ...
331
int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
332
                     1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
290 werner 333
int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
334
                       6,18,2,10,4,24,12,0,8,14,20,22};
335
 
336
 
287 werner 337
// sort function
338
bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
311 werner 339
{
340
    return s1.second < s2.second;
341
}
342
 
343
/**
344
*/
290 werner 345
void StandLoader::executeiLandInit(ResourceUnit *ru)
287 werner 346
{
290 werner 347
 
287 werner 348
    QPointF offset = ru->boundingBox().topLeft();
349
    QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);
350
 
351
    // a multimap holds a list for all trees.
352
    // key is the index of a 10x10m pixel within the resource unit
353
    QMultiMap<int, int> tree_map;
354
    QVector<QPair<int, double> > tcount; // counts
355
    for (int i=0;i<100;i++)
356
        tcount.push_back(QPair<int,double>(i,0.));
357
 
358
    int key;
311 werner 359
    double rand_val, rand_fraction;
312 werner 360
    int total_count = 0;
290 werner 361
    foreach(const InitFileItem &item, mInitItems) {
311 werner 362
        rand_fraction = fabs(double(item.density));
287 werner 363
        for (int i=0;i<item.count;i++) {
364
            // create trees
365
            int tree_idx = ru->newTreeIndex();
366
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
367
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
368
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
369
            tree.setSpecies(item.species);
384 werner 370
            if (item.age<=0)
388 werner 371
                tree.setAge(0,tree.height());
381 werner 372
            else
388 werner 373
                tree.setAge(item.age, tree.height());
287 werner 374
            tree.setRU(ru);
375
            tree.setup();
312 werner 376
            total_count++;
311 werner 377
 
378
            // calculate random value. "density" is from 1..-1.
379
            rand_val = mRandom->get();
380
            if (item.density<0)
381
                rand_val = 1. - rand_val;
312 werner 382
            rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
311 werner 383
 
287 werner 384
            // key: rank of target pixel
385
            // first: index of target pixel
386
            // second: sum of target pixel
311 werner 387
            key = limit(int(100*rand_val), 0, 99); // get from random number generator
287 werner 388
            tree_map.insert(tcount[key].first, tree_idx); // store tree in map
389
            tcount[key].second+=tree.basalArea(); // aggregate the basal area for each 10m pixel
312 werner 390
            if ( (total_count < 20 && i%2==0)
391
                || (total_count<100 && i%10==0 )
392
                || (i%30==0) ) {
287 werner 393
                qSort(tcount.begin(), tcount.end(), sortPairLessThan);
312 werner 394
            }
287 werner 395
        }
396
        qSort(tcount.begin(), tcount.end(), sortPairLessThan);
397
    }
398
 
399
    int bits, index, pos;
400
    int c;
401
    QList<int> trees;
402
    QPoint tree_pos;
289 werner 403
 
287 werner 404
    for (int i=0;i<100;i++) {
405
        trees = tree_map.values(i);
406
        c = trees.count();
290 werner 407
        QPointF pixel_center = ru->boundingBox().topLeft() + QPointF((i/10)*10. + 5., (i%10)*10. + 5.);
408
        if (!mModel->heightGrid()->valueAt(pixel_center).isValid()) {
409
            // no trees on that pixel: let trees die
410
            foreach(int tree_idx, trees) {
411
                ru->trees()[tree_idx].die();
412
            }
413
            continue;
414
        }
415
 
287 werner 416
        bits = 0;
417
        index = -1;
290 werner 418
        double r;
287 werner 419
        foreach(int tree_idx, trees) {
420
            if (c>18) {
421
                index = (index + 1)%25;
422
            } else {
423
                int stop=1000;
290 werner 424
                index = 0;
287 werner 425
                do {
290 werner 426
                    //r = drandom();
427
                    //if (r<0.5)  // skip position with a prob. of 50% -> adds a little "noise"
428
                    //    index++;
429
                    //index = (index + 1)%25; // increase and roll over
430
 
431
                    // search a random position
432
                    r = drandom();
433
                    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 434
                } while (isBitSet(bits, index)==true && stop--);
435
                if (!stop)
436
                    qDebug() << "executeiLandInit: found no free bit.";
437
                setBit(bits, index, true); // mark position as used
438
            }
290 werner 439
            // get position from fixed lists (one for even, one for uneven resource units)
440
            pos = ru->index()%2?evenlist[index]:unevenlist[index];
287 werner 441
            tree_pos = offsetIdx  // position of resource unit
442
                       + QPoint(5*(i/10), 5*(i%10)) // relative position of 10x10m pixel
443
                       + QPoint(pos/5, pos%5); // relative position within 10x10m pixel
444
            //qDebug() << tree_no++ << "to" << index;
445
            ru->trees()[tree_idx].setPosition(tree_pos);
446
        }
447
    }
448
}