WebSVN - public iLand - Diff - Rev 779 and 782


Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/standloader.cpp':
/********************************************************************************************
**    iLand - an individual based forest landscape and disturbance model
**    http://iland.boku.ac.at
**    Copyright (C) 2009-  Werner Rammer, Rupert Seidl
**
**    This program is free software: you can redistribute it and/or modify
**    it under the terms of the GNU General Public License as published by
**    the Free Software Foundation, either version 3 of the License, or
**    (at your option) any later version.
**
**    This program is distributed in the hope that it will be useful,
**    but WITHOUT ANY WARRANTY; without even the implied warranty of
**    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
**    GNU General Public License for more details.
**
**    You should have received a copy of the GNU General Public License
**    along with this program.  If not, see <http://www.gnu.org/licenses/>.
********************************************************************************************/

#include "global.h"
#include "standloader.h"


#include "grid.h"
#include "model.h"
#include "resourceunit.h"
#include "speciesset.h"

#include "helper.h"
#include "random.h"
#include "expression.h"
#include "expressionwrapper.h"
#include "environment.h"
#include "csvfile.h"
#include "mapgrid.h"
#include "snapshot.h"

/** @class StandLoader
    @ingroup tools
    loads (initializes) trees for a "stand" from various sources.
    StandLoader initializes trees on the landscape. It reads (usually) from text files, creates the
    trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid).

    See http://iland.boku.ac.at/initialize+trees
  */
// provide a mapping between "Picus"-style and "iLand"-style species Ids
QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";

StandLoader::~StandLoader()
{
    if (mRandom)
        delete mRandom;
    if (mHeightGridResponse)
        delete mHeightGridResponse;
}


void StandLoader::copyTrees()
{
    // we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards
    const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();
    ResourceUnit **p = ruGrid.begin();
    if (!p)
        throw IException("Standloader: invalid resource unit pointer!");
    ++p; // skip the first...
    const QVector<Tree> &tocopy = mModel->ru()->trees();
    for (; p!=ruGrid.end(); ++p) {
        QRectF rect = (*p)->boundingBox();
        foreach(const Tree& tree, tocopy) {
            Tree &newtree = (*p)->newTree();
            newtree = tree; // copy tree data...
            newtree.setPosition(tree.position()+rect.topLeft());
            newtree.setRU(*p);
            newtree.setNewId();
        }
    }
    if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";
}

/** main routine of the stand setup.
*/
void StandLoader::processInit()
{
    GlobalSettings *g = GlobalSettings::instance();
    XmlHelper xml(g->settings().node("model.initialization"));

    QString copy_mode = xml.value("mode", "copy");
    QString type = xml.value("type", "");
    QString  fileName = xml.value("file", "");

    bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);
    mHeightGridTries = xml.valueDouble("heightGrid.maxTries", 10.);
    QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function
    if (height_grid_enabled) {
        QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init");
        qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file;

        QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));
        if (!p->isValid()) {
            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));
        }
        height_grid.swap(p);
        mInitHeightGrid = height_grid.data();

        QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");
        mHeightGridResponse = new Expression(expr);
        mHeightGridResponse->setLinearizationEnabled(true);
    }

    Tree::resetStatistics();

    // one global init-file for the whole area:
    if (copy_mode=="single") {
        loadInitFile(fileName, type);
        evaluateDebugTrees();
        return;
    }


    // call a single tree init for each resource unit
    if (copy_mode=="unit") {
        foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
            ResourceUnit *ru = const_cast<ResourceUnit*>(const_ru);
            // set environment
            g->model()->environment()->setPosition(ru->boundingBox().center());
            type = xml.value("type", "");
            fileName = xml.value("file", "");
            if (fileName.isEmpty())
                continue;
            loadInitFile(fileName, type, 0, ru);
            if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
        }
        evaluateDebugTrees();
        return;
    }

    // map-modus: load a init file for each "polygon" in the standgrid
    if (copy_mode=="map") {
        if (!g->model()->standGrid() || !g->model()->standGrid()->isValid())
            throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));
        QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");

        CSVFile map_file(map_file_name);
        if (map_file.rowCount()==0)
            throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));
        int ikey = map_file.columnIndex("id");
        int ivalue = map_file.columnIndex("filename");
        if (ikey<0 || ivalue<0)
            throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));
        QString file_name;
        for (int i=0;i<map_file.rowCount();i++) {
            int key = map_file.value(i, ikey).toInt();
            if (key>0) {
                file_name = map_file.value(i, ivalue).toString();
                if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;
                if (!file_name.isEmpty())
                    loadInitFile(file_name, type, key, NULL);
            }
        }
        mInitHeightGrid = 0;
        evaluateDebugTrees();
        return;
    }
    if (copy_mode=="snapshot") {
        // load a snapshot from a file
        Snapshot shot;

        QString input_db = GlobalSettings::instance()->path(fileName);
        shot.loadSnapshot(input_db);
        return;
    }
    throw IException("StandLoader::processInit: invalid initalization.mode!");
}

void StandLoader::evaluateDebugTrees()
{
    // evaluate debugging
    QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
    int counter=0;
    if (!dbg_str.isEmpty()) {
        if (dbg_str == "debugstamp") {
            qDebug() << "debug_tree = debugstamp: try touching all trees...";
            // try to force an error if a stamp is invalid
            AllTreeIterator at(GlobalSettings::instance()->model());
            double total_offset=0.;
            while (Tree *t=at.next()) {
                total_offset += t->stamp()->offset();
                if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
                    qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
            }
            qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
            return;
        }
        TreeWrapper tw;
        Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
        AllTreeIterator at(GlobalSettings::instance()->model());
        double result;
        while (Tree *t = at.next()) {
            tw.setTree(t);
            result = dexp.execute();
            if (result) {
                t->enableDebugging();
                counter++;
            }
        }
        qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
    }
}

/// load a single init file. Calls loadPicusFile() or loadiLandFile()
/// @param fileName file to load
/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
{
    QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
    if (!QFile::exists(pathFileName))
        throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));

    if (type=="picus" || type=="single") {
        if (stand_id>0)
            throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
        return loadPicusFile(pathFileName, ru);
    }
    if (type=="iland" || type=="distribution")
        return loadiLandFile(pathFileName, ru, stand_id);

    throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
}

int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
{
    QString content = Helper::loadTextFile(fileName);
    if (content.isEmpty()) {
        qDebug() << "file not found: " + fileName;
        return 0;
    }
    return loadSingleTreeList(content, ru, fileName);
}

/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.
  returns the number of loaded trees.
  */
int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
{
    if (!ru)
        ru = mModel->ru();
    Q_ASSERT(ru!=0);

    QPointF offset = ru->boundingBox().topLeft();
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU

    QString my_content(content);
    // cut out the <trees> </trees> part if present
    if (content.contains("<trees>")) {
        QRegExp rx(".*<trees>(.*)</trees>.*");
        rx.indexIn(content, 0);
        if (rx.capturedTexts().count()<1)
            return 0;
        my_content = rx.cap(1).trimmed();
    }

    QStringList lines=my_content.split('\n');
    if (lines.count()<2)
        return 0;
    // drop comments
    while (!lines.isEmpty() && lines.front().startsWith('#') )
        lines.pop_front();
    while (!lines.isEmpty() && lines.last().isEmpty())
        lines.removeLast();

    char sep='\t';
    if (!lines[0].contains(sep))
        sep=';';
    QStringList headers = lines[0].trimmed().split(sep);

    int iID = headers.indexOf("id");
    int iX = headers.indexOf("x");
    int iY = headers.indexOf("y");
    int iBhd = headers.indexOf("bhdfrom");
    if (iBhd<0)
        iBhd = headers.indexOf("dbh");
    double height_conversion = 100.;
    int iHeight = headers.indexOf("treeheight");
    if (iHeight<0) {
        iHeight = headers.indexOf("height");
        height_conversion = 1.; // in meter
    }
    int iSpecies = headers.indexOf("species");
    int iAge = headers.indexOf("age");
    if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1)
        throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));

    double dbh;
    bool ok;
    int cnt=0;
    QString speciesid;
    for (int i=1;i<lines.count();i++) {
        QString &line = lines[i];
        dbh = line.section(sep, iBhd, iBhd).toDouble();
        if (dbh<5.)
            continue;

        QPointF f;
        if (iX>=0 && iY>=0) {
           f.setX( line.section(sep, iX, iX).toDouble() );
           f.setY( line.section(sep, iY, iY).toDouble() );
           f+=offset;

        }
        // position valid?
        if (!mModel->heightGrid()->valueAt(f).isValid())
            continue;
        Tree &tree = ru->newTree();
        tree.setPosition(f);
        if (iID>=0)
            tree.setId(line.section(sep, iID, iID).toInt() );

        tree.setDbh(dbh);
        tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary

        speciesid = line.section(sep, iSpecies, iSpecies).trimmed();
        int picusid = speciesid.toInt(&ok);
        if (ok) {
            int idx = picusSpeciesIds.indexOf(picusid);
            if (idx==-1)
                throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
            speciesid = iLandSpeciesIds[idx];
        }
        Species *s = speciesSet->species(speciesid);
        if (!ru || !s)
            throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
        tree.setSpecies(s);

        ok = true;
        if (iAge>=0)
           tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a *real* age
        if (iAge<0 || !ok || tree.age()==0)
           tree.setAge(0, tree.height()); // no real tree age available

        tree.setRU(ru);
        tree.setup();
        cnt++;
    }
    return cnt;
    //qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
    //qDebug() << "lines: " << lines;
}

/** initialize trees on a resource unit based on dbh distributions.
  use a fairly clever algorithm to determine tree positions.
  see http://iland.boku.ac.at/initialize+trees
  @param content tree init file (including headers) in a string
  @param ru resource unit
  @param fileName source file name (for error reporting)
  @return number of trees added
  */
int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
{
    if (!ru)
        ru = mModel->ru();
    Q_ASSERT(ru!=0);
    SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
    Q_ASSERT(speciesSet!=0);

    //DebugTimer t("StandLoader::loadiLandFile");
    CSVFile infile;
    infile.loadFromString(content);

    int icount = infile.columnIndex("count");
    int ispecies = infile.columnIndex("species");
    int idbh_from = infile.columnIndex("dbh_from");
    int idbh_to = infile.columnIndex("dbh_to");
    int ihd = infile.columnIndex("hd");
    int iage = infile.columnIndex("age");
    int idensity = infile.columnIndex("density");
    if (icount<0 || ispecies<0 || idbh_from<0 || idbh_to<0 || ihd<0 || iage<0)
        throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));

    mInitItems.clear();
    InitFileItem item;
    bool ok;
    int total_count = 0;
    for (int row=0;row<infile.rowCount();row++) {
        item.count = infile.value(row, icount).toInt();
        total_count += item.count;
        item.dbh_from = infile.value(row, idbh_from).toDouble();
        item.dbh_to = infile.value(row, idbh_to).toDouble();
        item.hd = infile.value(row, ihd).toDouble();
        if (item.hd==0. || item.dbh_from / 100. * item.hd < 4.)
 
            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) ;
            //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) );
        ok = true;
        if (iage>=0)
            item.age = infile.value(row, iage).toInt(&ok);
        if (iage<0 || !ok)
            item.age = 0;

        item.species = speciesSet->species(infile.value(row, ispecies).toString());
        if (idensity>=0)
            item.density = infile.value(row, idensity).toDouble();
        else
            item.density = 0.;
        if (item.density<-1 || item.density>1)
            throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
                             .arg(item.density)
                             .arg(fileName)
                             .arg(row));
        if (!item.species) {
            throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
                             .arg(infile.value(row, ispecies).toString())
                             .arg(fileName)
                             .arg(row));
        }
        mInitItems.push_back(item);
    }
    // setup the random distribution
    QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
    if (logLevelInfo())  qDebug() << "density function:" << density_func;
    if (!mRandom || (mRandom->densityFunction()!= density_func)) {
        if (mRandom)
            delete mRandom;
        mRandom=new RandomCustomPDF(density_func);
        if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
    }
    if (stand_id>0) {
        // execute stand based initialization
        executeiLandInitStand(stand_id);
    } else {
        // exeucte the initialization based on single resource units
        executeiLandInit(ru);
        ru->cleanTreeList();
    }
    return total_count;

}

int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)
{
    if (!QFile::exists(fileName))
        throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
    QString content = Helper::loadTextFile(fileName);
    return loadDistributionList(content, ru, stand_id, fileName);
}

// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
// e.g. 12 = centerpixel, 0: upper left corner, ...
int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
                     1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
                       6,18,2,10,4,24,12,0,8,14,20,22};



// sort function
bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
{
    return s1.second < s2.second;
}

struct SInitPixel {
    double basal_area; // accumulated basal area
    QPoint pixelOffset; // location of the pixel
    ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to
    double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)
    SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}
};

bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)
{
    return s1.basal_area < s2.basal_area;
}


/**
*/

void StandLoader::executeiLandInit(ResourceUnit *ru)
{

    QPointF offset = ru->boundingBox().topLeft();
    QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);

    // a multimap holds a list for all trees.
    // key is the index of a 10x10m pixel within the resource unit
    QMultiMap<int, int> tree_map;
    //QHash<int,SInitPixel> tcount;

    QVector<QPair<int, double> > tcount; // counts
    for (int i=0;i<100;i++)
        tcount.push_back(QPair<int,double>(i,0.));

    int key;
    double rand_val, rand_fraction;
    int total_count = 0;
    foreach(const InitFileItem &item, mInitItems) {
        rand_fraction = fabs(double(item.density));
        for (int i=0;i<item.count;i++) {
            // create trees
            int tree_idx = ru->newTreeIndex();
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
            tree.setSpecies(item.species);
            if (item.age<=0)
                tree.setAge(0,tree.height());
            else
                tree.setAge(item.age, tree.height());
            tree.setRU(ru);
            tree.setup();
            total_count++;

            // calculate random value. "density" is from 1..-1.
            rand_val = mRandom->get();
            if (item.density<0)
                rand_val = 1. - rand_val;
            rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);

            // key: rank of target pixel
            // first: index of target pixel
            // second: sum of target pixel
            key = limit(int(100*rand_val), 0, 99); // get from random number generator
            tree_map.insert(tcount[key].first, tree_idx); // store tree in map
            tcount[key].second+=tree.basalArea(); // aggregate the basal area for each 10m pixel
            if ( (total_count < 20 && i%2==0)
                || (total_count<100 && i%10==0 )
                || (i%30==0) ) {
                qSort(tcount.begin(), tcount.end(), sortPairLessThan);
            }
        }
        qSort(tcount.begin(), tcount.end(), sortPairLessThan);
    }

    int bits, index, pos;
    int c;
    QList<int> trees;
    QPoint tree_pos;

    for (int i=0;i<100;i++) {
        trees = tree_map.values(i);
        c = trees.count();
        QPointF pixel_center = ru->boundingBox().topLeft() + QPointF((i/10)*10. + 5., (i%10)*10. + 5.);
        if (!mModel->heightGrid()->valueAt(pixel_center).isValid()) {
            // no trees on that pixel: let trees die
            foreach(int tree_idx, trees) {
                ru->trees()[tree_idx].die();
            }
            continue;
        }

        bits = 0;
        index = -1;
        double r;
        foreach(int tree_idx, trees) {
            if (c>18) {
                index = (index + 1)%25;
            } else {
                int stop=1000;
                index = 0;
                do {
                    //r = drandom();
                    //if (r<0.5)  // skip position with a prob. of 50% -> adds a little "noise"
                    //    index++;
                    //index = (index + 1)%25; // increase and roll over

                    // search a random position
                    r = drandom();
                    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)
                } while (isBitSet(bits, index)==true && stop--);
                if (!stop)
                    qDebug() << "executeiLandInit: found no free bit.";
                setBit(bits, index, true); // mark position as used
            }
            // get position from fixed lists (one for even, one for uneven resource units)
            pos = ru->index()%2?evenlist[index]:unevenlist[index];
            tree_pos = offsetIdx  // position of resource unit
                       + QPoint(5*(i/10), 5*(i%10)) // relative position of 10x10m pixel
                       + QPoint(pos/5, pos%5); // relative position within 10x10m pixel
            //qDebug() << tree_no++ << "to" << index;
            ru->trees()[tree_idx].setPosition(tree_pos);
        }
    }
}

// provide a hashing function for the QPoint type (needed from stand init function below)
inline uint qHash(const QPoint &key)
 {
     return qHash(key.x()) ^ qHash(key.y());
 }


// Initialization routine based on a stand map.
// Basically a list of 10m pixels for a given stand is retrieved
// and the filled with the same procedure as the resource unit based init
// see http://iland.boku.ac.at/initialize+trees
void StandLoader::executeiLandInitStand(int stand_id)
{

    const MapGrid *grid = GlobalSettings::instance()->model()->standGrid();

    // get a list of positions of all pixels that belong to our stand
    QList<int> indices = grid->gridIndices(stand_id);
    if (indices.isEmpty()) {
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
        return;
    }
    // a multiHash holds a list for all trees.
    // key is the location of the 10x10m pixel
    QMultiHash<QPoint, int> tree_map;
    QList<SInitPixel> pixel_list; // working list of all 10m pixels
    pixel_list.reserve(indices.size());

    foreach (int i, indices) {
       SInitPixel p;
       p.pixelOffset = grid->grid().indexOf(i); // index in the 10m grid
       p.resource_unit = GlobalSettings::instance()->model()->ru( grid->grid().cellCenterPoint(p.pixelOffset));
       if (mInitHeightGrid)
           p.h_max = mInitHeightGrid->grid().constValueAtIndex(p.pixelOffset);
       pixel_list.append(p);
    }
    double area_factor = grid->area(stand_id) / 10000.;

    int key=0;
    double rand_val, rand_fraction;
    int total_count = 0;
    int total_tries = 0;
    int total_misses = 0;
    if (mInitHeightGrid && !mHeightGridResponse)
        throw IException("executeiLandInitStand: trying to initialize with height grid but without response function.");
    foreach(const InitFileItem &item, mInitItems) {
        rand_fraction = fabs(double(item.density));
        int count = item.count * area_factor + 0.5; // round
        double init_max_height = item.dbh_to/100. * item.hd;
        for (int i=0;i<count;i++) {

            bool found = false;
            int tries = mHeightGridTries;
            while (!found &&--tries) {
                // calculate random value. "density" is from 1..-1.
                rand_val = mRandom->get();
                if (item.density<0)
                    rand_val = 1. - rand_val;
                rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
                ++total_tries;

                // key: rank of target pixel
                key = limit(int(pixel_list.count()*rand_val), 0, pixel_list.count()-1); // get from random number generator

                if (mInitHeightGrid) {
                    // calculate how good the selected pixel fits w.r.t. the predefined height
                    double p_value = pixel_list[key].h_max>0.?mHeightGridResponse->calculate(init_max_height/pixel_list[key].h_max):0.;
                    if (drandom() < p_value)
                        found = true;
                } else {
                    found = true;
                }
            }
            if (tries<0) ++total_misses;

            // create a tree
            ResourceUnit *ru = pixel_list[key].resource_unit;
            int tree_idx = ru->newTreeIndex();
            Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
            tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
            tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
            tree.setSpecies(item.species);
            if (item.age<=0)
                tree.setAge(0,tree.height());
            else
                tree.setAge(item.age, tree.height());
            tree.setRU(ru);
            tree.setup();
            total_count++;

            // store in the multiHash the position of the pixel and the tree_idx in the resepctive resource unit
            tree_map.insert(pixel_list[key].pixelOffset, tree_idx);
            pixel_list[key].basal_area+=tree.basalArea(); // aggregate the basal area for each 10m pixel

            // resort list
            if ( (total_count < 20 && i%2==0)
                || (total_count<100 && i%10==0 )
                || (i%30==0) ) {
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
            }
        }
        qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
    }
    if (total_misses>0 || total_tries > total_count) {
        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);
    }

    int bits, index, pos;
    int c;
    QList<int> trees;
    QPoint tree_pos;

    foreach(const SInitPixel &p, pixel_list) {
        trees = tree_map.values(p.pixelOffset);
        c = trees.count();
        bits = 0;
        index = -1;
        double r;
        foreach(int tree_idx, trees) {
            if (c>18) {
                index = (index + 1)%25;
            } else {
                int stop=1000;
                index = 0;
                do {
                    // search a random position
                    r = drandom();
                    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)
                } while (isBitSet(bits, index)==true && stop--);
                if (!stop)
                    qDebug() << "executeiLandInit: found no free bit.";
                setBit(bits, index, true); // mark position as used
            }
            // get position from fixed lists (one for even, one for uneven resource units)
            pos = p.resource_unit->index()%2?evenlist[index]:unevenlist[index];
            tree_pos = p.pixelOffset * cPxPerHeight; // convert to LIF index
            tree_pos += QPoint(pos/cPxPerHeight, pos%cPxPerHeight);

            p.resource_unit->trees()[tree_idx].setPosition(tree_pos);
            // test if tree position is valid..
            if (!GlobalSettings::instance()->model()->grid()->isIndexValid(tree_pos))
                qDebug() << "Standloader: invalid position!";
        }
    }
    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;

}

/// a (hacky) way of adding saplings of a certain age to a stand defined by 'stand_id'.
int StandLoader::loadSaplings(const QString &content, int stand_id, const QString &fileName)
{
    Q_UNUSED(fileName);
    const MapGrid *stand_grid;
    if (mCurrentMap)
        stand_grid = mCurrentMap; // if set
    else
        stand_grid = GlobalSettings::instance()->model()->standGrid(); // default

    QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels
    if (indices.isEmpty()) {
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
        return -1;
    }
    double area_factor = stand_grid->area(stand_id) / 10000.; // multiplier for grid (e.g. 2 if stand has area of 2 hectare)

    // parse the content of the init-file
    // species
    CSVFile init;
    init.loadFromString(content);
    int ispecies = init.columnIndex("species");
    int icount = init.columnIndex("count");
    int iheight = init.columnIndex("height");
    int iage = init.columnIndex("age");
    if (ispecies==-1 || icount==-1)
        throw IException("Error while loading saplings: columns 'species' or 'count' are missing!!");

    const SpeciesSet *set = GlobalSettings::instance()->model()->ru()->speciesSet();
    double height, age;
    int total = 0;
    for (int row=0;row<init.rowCount();++row) {
        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)
        const Species *species = set->species(init.value(row, ispecies).toString());
        if (!species)
            throw IException(QString("Error while loading saplings: invalid species '%1'.").arg(init.value(row, ispecies).toString()));
        height = iheight==-1?0.05: init.value(row, iheight).toDouble();
        age = iage==-1?1:init.value(row,iage).toDouble();

        int misses = 0;
        int hits = 0;
        while (hits < pxcount) {
           int rnd_index = irandom(0, indices.count()-1);
           QPoint offset=stand_grid->grid().indexOf(indices[rnd_index]);
           ResourceUnit *ru = GlobalSettings::instance()->model()->ru(stand_grid->grid().cellCenterPoint(offset));
           //
           offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
           int in_p = irandom(0, cPxPerHeight*cPxPerHeight-1); // index of lif-pixel
           offset += QPoint(in_p / cPxPerHeight, in_p % cPxPerHeight);
           if (ru->saplingHeightForInit(offset) > height) {
               misses++;
           } else {
               // ok
               hits++;
               ru->resourceUnitSpecies(species).changeSapling().addSapling(offset);
           }
           if (misses > 3*pxcount) {
               qDebug() << "tried to add" << pxcount << "saplings at stand" << stand_id << "but failed in finding enough free positions. Added" << hits << "and stopped.";
               break;
           }
        }
        total += hits;

    }
    return total;
}
 

Rev 779	Rev 782
1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/standloader.cpp':	1	Redirecting to URL 'https://iland.boku.ac.at/svn/iland/tags/release_1.0/src/core/standloader.cpp':
2	/********************************************************************************************	2	/********************************************************************************************
3	** iLand - an individual based forest landscape and disturbance model	3	** iLand - an individual based forest landscape and disturbance model
4	** http://iland.boku.ac.at	4	** http://iland.boku.ac.at
5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl	5	** Copyright (C) 2009- Werner Rammer, Rupert Seidl
6	**	6	**
7	** This program is free software: you can redistribute it and/or modify	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	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	9	** the Free Software Foundation, either version 3 of the License, or
10	** (at your option) any later version.	10	** (at your option) any later version.
11	**	11	**
12	** This program is distributed in the hope that it will be useful,	12	** This program is distributed in the hope that it will be useful,
13	** but WITHOUT ANY WARRANTY; without even the implied warranty of	13	** but WITHOUT ANY WARRANTY; without even the implied warranty of
14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the	14	** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15	** GNU General Public License for more details.	15	** GNU General Public License for more details.
16	**	16	**
17	** You should have received a copy of the GNU General Public License	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/>.	18	** along with this program. If not, see <http://www.gnu.org/licenses/>.
19	********************************************************************************************/	19	********************************************************************************************/
20		20
21	#include "global.h"	21	#include "global.h"
22	#include "standloader.h"	22	#include "standloader.h"
23		23
24		24
25	#include "grid.h"	25	#include "grid.h"
26	#include "model.h"	26	#include "model.h"
27	#include "resourceunit.h"	27	#include "resourceunit.h"
28	#include "speciesset.h"	28	#include "speciesset.h"
29		29
30	#include "helper.h"	30	#include "helper.h"
31	#include "random.h"	31	#include "random.h"
32	#include "expression.h"	32	#include "expression.h"
33	#include "expressionwrapper.h"	33	#include "expressionwrapper.h"
34	#include "environment.h"	34	#include "environment.h"
35	#include "csvfile.h"	35	#include "csvfile.h"
36	#include "mapgrid.h"	36	#include "mapgrid.h"
37	#include "snapshot.h"	37	#include "snapshot.h"
38		38
39	/** @class StandLoader	39	/** @class StandLoader
40	@ingroup tools	40	@ingroup tools
41	loads (initializes) trees for a "stand" from various sources.	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	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).	43	trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid).
44		44
45	See http://iland.boku.ac.at/initialize+trees	45	See http://iland.boku.ac.at/initialize+trees
46	*/	46	*/
47	// provide a mapping between "Picus"-style and "iLand"-style species Ids	47	// provide a mapping between "Picus"-style and "iLand"-style species Ids
48	QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;	48	QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
49	QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";	49	QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";
50		50
51	StandLoader::~StandLoader()	51	StandLoader::~StandLoader()
52	{	52	{
53	if (mRandom)	53	if (mRandom)
54	delete mRandom;	54	delete mRandom;
55	if (mHeightGridResponse)	55	if (mHeightGridResponse)
56	delete mHeightGridResponse;	56	delete mHeightGridResponse;
57	}	57	}
58		58
59		59
60	void StandLoader::copyTrees()	60	void StandLoader::copyTrees()
61	{	61	{
62	// we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards	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();	63	const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();
64	ResourceUnit **p = ruGrid.begin();	64	ResourceUnit **p = ruGrid.begin();
65	if (!p)	65	if (!p)
66	throw IException("Standloader: invalid resource unit pointer!");	66	throw IException("Standloader: invalid resource unit pointer!");
67	++p; // skip the first...	67	++p; // skip the first...
68	const QVector<Tree> &tocopy = mModel->ru()->trees();	68	const QVector<Tree> &tocopy = mModel->ru()->trees();
69	for (; p!=ruGrid.end(); ++p) {	69	for (; p!=ruGrid.end(); ++p) {
70	QRectF rect = (*p)->boundingBox();	70	QRectF rect = (*p)->boundingBox();
71	foreach(const Tree& tree, tocopy) {	71	foreach(const Tree& tree, tocopy) {
72	Tree &newtree = (*p)->newTree();	72	Tree &newtree = (*p)->newTree();
73	newtree = tree; // copy tree data...	73	newtree = tree; // copy tree data...
74	newtree.setPosition(tree.position()+rect.topLeft());	74	newtree.setPosition(tree.position()+rect.topLeft());
75	newtree.setRU(*p);	75	newtree.setRU(*p);
76	newtree.setNewId();	76	newtree.setNewId();
77	}	77	}
78	}	78	}
79	if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";	79	if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";
80	}	80	}
81		81
82	/** main routine of the stand setup.	82	/** main routine of the stand setup.
83	*/	83	*/
84	void StandLoader::processInit()	84	void StandLoader::processInit()
85	{	85	{
86	GlobalSettings *g = GlobalSettings::instance();	86	GlobalSettings *g = GlobalSettings::instance();
87	XmlHelper xml(g->settings().node("model.initialization"));	87	XmlHelper xml(g->settings().node("model.initialization"));
88		88
89	QString copy_mode = xml.value("mode", "copy");	89	QString copy_mode = xml.value("mode", "copy");
90	QString type = xml.value("type", "");	90	QString type = xml.value("type", "");
91	QString fileName = xml.value("file", "");	91	QString fileName = xml.value("file", "");
92		92
93	bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);	93	bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);
94	mHeightGridTries = xml.valueDouble("heightGrid.maxTries", 10.);	94	mHeightGridTries = xml.valueDouble("heightGrid.maxTries", 10.);
95	QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function	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) {	96	if (height_grid_enabled) {
97	QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init");	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;	98	qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file;
99		99
100	QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));	100	QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));
101	if (!p->isValid()) {	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));	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	}	103	}
104	height_grid.swap(p);	104	height_grid.swap(p);
105	mInitHeightGrid = height_grid.data();	105	mInitHeightGrid = height_grid.data();
106		106
107	QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");	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);	108	mHeightGridResponse = new Expression(expr);
109	mHeightGridResponse->setLinearizationEnabled(true);	109	mHeightGridResponse->setLinearizationEnabled(true);
110	}	110	}
111		111
112	Tree::resetStatistics();	112	Tree::resetStatistics();
113		113
114	// one global init-file for the whole area:	114	// one global init-file for the whole area:
115	if (copy_mode=="single") {	115	if (copy_mode=="single") {
116	loadInitFile(fileName, type);	116	loadInitFile(fileName, type);
117	evaluateDebugTrees();	117	evaluateDebugTrees();
118	return;	118	return;
119	}	119	}
120		120
121		121
122	// call a single tree init for each resource unit	122	// call a single tree init for each resource unit
123	if (copy_mode=="unit") {	123	if (copy_mode=="unit") {
124	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {	124	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
125	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);	125	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);
126	// set environment	126	// set environment
127	g->model()->environment()->setPosition(ru->boundingBox().center());	127	g->model()->environment()->setPosition(ru->boundingBox().center());
128	type = xml.value("type", "");	128	type = xml.value("type", "");
129	fileName = xml.value("file", "");	129	fileName = xml.value("file", "");
130	if (fileName.isEmpty())	130	if (fileName.isEmpty())
131	continue;	131	continue;
132	loadInitFile(fileName, type, 0, ru);	132	loadInitFile(fileName, type, 0, ru);
133	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";	133	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
134	}	134	}
135	evaluateDebugTrees();	135	evaluateDebugTrees();
136	return;	136	return;
137	}	137	}
138		138
139	// map-modus: load a init file for each "polygon" in the standgrid	139	// map-modus: load a init file for each "polygon" in the standgrid
140	if (copy_mode=="map") {	140	if (copy_mode=="map") {
141	if (!g->model()->standGrid() \|\| !g->model()->standGrid()->isValid())	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)"));	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");	143	QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");
144		144
145	CSVFile map_file(map_file_name);	145	CSVFile map_file(map_file_name);
146	if (map_file.rowCount()==0)	146	if (map_file.rowCount()==0)
147	throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));	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");	148	int ikey = map_file.columnIndex("id");
149	int ivalue = map_file.columnIndex("filename");	149	int ivalue = map_file.columnIndex("filename");
150	if (ikey<0 \|\| ivalue<0)	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));	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;	152	QString file_name;
153	for (int i=0;i<map_file.rowCount();i++) {	153	for (int i=0;i<map_file.rowCount();i++) {
154	int key = map_file.value(i, ikey).toInt();	154	int key = map_file.value(i, ikey).toInt();
155	if (key>0) {	155	if (key>0) {
156	file_name = map_file.value(i, ivalue).toString();	156	file_name = map_file.value(i, ivalue).toString();
157	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;	157	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;
158	if (!file_name.isEmpty())	158	if (!file_name.isEmpty())
159	loadInitFile(file_name, type, key, NULL);	159	loadInitFile(file_name, type, key, NULL);
160	}	160	}
161	}	161	}
162	mInitHeightGrid = 0;	162	mInitHeightGrid = 0;
163	evaluateDebugTrees();	163	evaluateDebugTrees();
164	return;	164	return;
165	}	165	}
166	if (copy_mode=="snapshot") {	166	if (copy_mode=="snapshot") {
167	// load a snapshot from a file	167	// load a snapshot from a file
168	Snapshot shot;	168	Snapshot shot;
169		169
170	QString input_db = GlobalSettings::instance()->path(fileName);	170	QString input_db = GlobalSettings::instance()->path(fileName);
171	shot.loadSnapshot(input_db);	171	shot.loadSnapshot(input_db);
172	return;	172	return;
173	}	173	}
174	throw IException("StandLoader::processInit: invalid initalization.mode!");	174	throw IException("StandLoader::processInit: invalid initalization.mode!");
175	}	175	}
176		176
177	void StandLoader::evaluateDebugTrees()	177	void StandLoader::evaluateDebugTrees()
178	{	178	{
179	// evaluate debugging	179	// evaluate debugging
180	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");	180	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
181	int counter=0;	181	int counter=0;
182	if (!dbg_str.isEmpty()) {	182	if (!dbg_str.isEmpty()) {
183	if (dbg_str == "debugstamp") {	183	if (dbg_str == "debugstamp") {
184	qDebug() << "debug_tree = debugstamp: try touching all trees...";	184	qDebug() << "debug_tree = debugstamp: try touching all trees...";
185	// try to force an error if a stamp is invalid	185	// try to force an error if a stamp is invalid
186	AllTreeIterator at(GlobalSettings::instance()->model());	186	AllTreeIterator at(GlobalSettings::instance()->model());
187	double total_offset=0.;	187	double total_offset=0.;
188	while (Tree *t=at.next()) {	188	while (Tree *t=at.next()) {
189	total_offset += t->stamp()->offset();	189	total_offset += t->stamp()->offset();
190	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))	190	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
191	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";	191	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
192	}	192	}
193	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";	193	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
194	return;	194	return;
195	}	195	}
196	TreeWrapper tw;	196	TreeWrapper tw;
197	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables	197	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
198	AllTreeIterator at(GlobalSettings::instance()->model());	198	AllTreeIterator at(GlobalSettings::instance()->model());
199	double result;	199	double result;
200	while (Tree *t = at.next()) {	200	while (Tree *t = at.next()) {
201	tw.setTree(t);	201	tw.setTree(t);
202	result = dexp.execute();	202	result = dexp.execute();
203	if (result) {	203	if (result) {
204	t->enableDebugging();	204	t->enableDebugging();
205	counter++;	205	counter++;
206	}	206	}
207	}	207	}
208	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";	208	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
209	}	209	}
210	}	210	}
211		211
212	/// load a single init file. Calls loadPicusFile() or loadiLandFile()	212	/// load a single init file. Calls loadPicusFile() or loadiLandFile()
213	/// @param fileName file to load	213	/// @param fileName file to load
214	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"	214	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
215	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)	215	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
216	{	216	{
217	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");	217	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
218	if (!QFile::exists(pathFileName))	218	if (!QFile::exists(pathFileName))
219	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));	219	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
220		220
221	if (type=="picus" \|\| type=="single") {	221	if (type=="picus" \|\| type=="single") {
222	if (stand_id>0)	222	if (stand_id>0)
223	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);	223	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
224	return loadPicusFile(pathFileName, ru);	224	return loadPicusFile(pathFileName, ru);
225	}	225	}
226	if (type=="iland" \|\| type=="distribution")	226	if (type=="iland" \|\| type=="distribution")
227	return loadiLandFile(pathFileName, ru, stand_id);	227	return loadiLandFile(pathFileName, ru, stand_id);
228		228
229	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);	229	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
230	}	230	}
231		231
232	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)	232	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
233	{	233	{
234	QString content = Helper::loadTextFile(fileName);	234	QString content = Helper::loadTextFile(fileName);
235	if (content.isEmpty()) {	235	if (content.isEmpty()) {
236	qDebug() << "file not found: " + fileName;	236	qDebug() << "file not found: " + fileName;
237	return 0;	237	return 0;
238	}	238	}
239	return loadSingleTreeList(content, ru, fileName);	239	return loadSingleTreeList(content, ru, fileName);
240	}	240	}
241		241
242	/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.	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.	243	returns the number of loaded trees.
244	*/	244	*/
245	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)	245	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
246	{	246	{
247	if (!ru)	247	if (!ru)
248	ru = mModel->ru();	248	ru = mModel->ru();
249	Q_ASSERT(ru!=0);	249	Q_ASSERT(ru!=0);
250		250
251	QPointF offset = ru->boundingBox().topLeft();	251	QPointF offset = ru->boundingBox().topLeft();
252	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	252	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
253		253
254	QString my_content(content);	254	QString my_content(content);
255	// cut out the <trees> </trees> part if present	255	// cut out the <trees> </trees> part if present
256	if (content.contains("<trees>")) {	256	if (content.contains("<trees>")) {
257	QRegExp rx(".<trees>(.)</trees>.*");	257	QRegExp rx(".<trees>(.)</trees>.*");
258	rx.indexIn(content, 0);	258	rx.indexIn(content, 0);
259	if (rx.capturedTexts().count()<1)	259	if (rx.capturedTexts().count()<1)
260	return 0;	260	return 0;
261	my_content = rx.cap(1).trimmed();	261	my_content = rx.cap(1).trimmed();
262	}	262	}
263		263
264	QStringList lines=my_content.split('\n');	264	QStringList lines=my_content.split('\n');
265	if (lines.count()<2)	265	if (lines.count()<2)
266	return 0;	266	return 0;
267	// drop comments	267	// drop comments
268	while (!lines.isEmpty() && lines.front().startsWith('#') )	268	while (!lines.isEmpty() && lines.front().startsWith('#') )
269	lines.pop_front();	269	lines.pop_front();
270	while (!lines.isEmpty() && lines.last().isEmpty())	270	while (!lines.isEmpty() && lines.last().isEmpty())
271	lines.removeLast();	271	lines.removeLast();
272		272
273	char sep='\t';	273	char sep='\t';
274	if (!lines[0].contains(sep))	274	if (!lines[0].contains(sep))
275	sep=';';	275	sep=';';
276	QStringList headers = lines[0].trimmed().split(sep);	276	QStringList headers = lines[0].trimmed().split(sep);
277		277
278	int iID = headers.indexOf("id");	278	int iID = headers.indexOf("id");
279	int iX = headers.indexOf("x");	279	int iX = headers.indexOf("x");
280	int iY = headers.indexOf("y");	280	int iY = headers.indexOf("y");
281	int iBhd = headers.indexOf("bhdfrom");	281	int iBhd = headers.indexOf("bhdfrom");
282	if (iBhd<0)	282	if (iBhd<0)
283	iBhd = headers.indexOf("dbh");	283	iBhd = headers.indexOf("dbh");
284	double height_conversion = 100.;	284	double height_conversion = 100.;
285	int iHeight = headers.indexOf("treeheight");	285	int iHeight = headers.indexOf("treeheight");
286	if (iHeight<0) {	286	if (iHeight<0) {
287	iHeight = headers.indexOf("height");	287	iHeight = headers.indexOf("height");
288	height_conversion = 1.; // in meter	288	height_conversion = 1.; // in meter
289	}	289	}
290	int iSpecies = headers.indexOf("species");	290	int iSpecies = headers.indexOf("species");
291	int iAge = headers.indexOf("age");	291	int iAge = headers.indexOf("age");
292	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)	292	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)
293	throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));	293	throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
294		294
295	double dbh;	295	double dbh;
296	bool ok;	296	bool ok;
297	int cnt=0;	297	int cnt=0;
298	QString speciesid;	298	QString speciesid;
299	for (int i=1;i<lines.count();i++) {	299	for (int i=1;i<lines.count();i++) {
300	QString &line = lines[i];	300	QString &line = lines[i];
301	dbh = line.section(sep, iBhd, iBhd).toDouble();	301	dbh = line.section(sep, iBhd, iBhd).toDouble();
302	if (dbh<5.)	302	if (dbh<5.)
303	continue;	303	continue;
304		304
305	QPointF f;	305	QPointF f;
306	if (iX>=0 && iY>=0) {	306	if (iX>=0 && iY>=0) {
307	f.setX( line.section(sep, iX, iX).toDouble() );	307	f.setX( line.section(sep, iX, iX).toDouble() );
308	f.setY( line.section(sep, iY, iY).toDouble() );	308	f.setY( line.section(sep, iY, iY).toDouble() );
309	f+=offset;	309	f+=offset;
310		310
311	}	311	}
312	// position valid?	312	// position valid?
313	if (!mModel->heightGrid()->valueAt(f).isValid())	313	if (!mModel->heightGrid()->valueAt(f).isValid())
314	continue;	314	continue;
315	Tree &tree = ru->newTree();	315	Tree &tree = ru->newTree();
316	tree.setPosition(f);	316	tree.setPosition(f);
317	if (iID>=0)	317	if (iID>=0)
318	tree.setId(line.section(sep, iID, iID).toInt() );	318	tree.setId(line.section(sep, iID, iID).toInt() );
319		319
320	tree.setDbh(dbh);	320	tree.setDbh(dbh);
321	tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary	321	tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
322		322
323	speciesid = line.section(sep, iSpecies, iSpecies).trimmed();	323	speciesid = line.section(sep, iSpecies, iSpecies).trimmed();
324	int picusid = speciesid.toInt(&ok);	324	int picusid = speciesid.toInt(&ok);
325	if (ok) {	325	if (ok) {
326	int idx = picusSpeciesIds.indexOf(picusid);	326	int idx = picusSpeciesIds.indexOf(picusid);
327	if (idx==-1)	327	if (idx==-1)
328	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));	328	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
329	speciesid = iLandSpeciesIds[idx];	329	speciesid = iLandSpeciesIds[idx];
330	}	330	}
331	Species *s = speciesSet->species(speciesid);	331	Species *s = speciesSet->species(speciesid);
332	if (!ru \|\| !s)	332	if (!ru \|\| !s)
333	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));	333	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
334	tree.setSpecies(s);	334	tree.setSpecies(s);
335		335
336	ok = true;	336	ok = true;
337	if (iAge>=0)	337	if (iAge>=0)
338	tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a real age	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)	339	if (iAge<0 \|\| !ok \|\| tree.age()==0)
340	tree.setAge(0, tree.height()); // no real tree age available	340	tree.setAge(0, tree.height()); // no real tree age available
341		341
342	tree.setRU(ru);	342	tree.setRU(ru);
343	tree.setup();	343	tree.setup();
344	cnt++;	344	cnt++;
345	}	345	}
346	return cnt;	346	return cnt;
347	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";	347	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
348	//qDebug() << "lines: " << lines;	348	//qDebug() << "lines: " << lines;
349	}	349	}
350		350
351	/** initialize trees on a resource unit based on dbh distributions.	351	/** initialize trees on a resource unit based on dbh distributions.
352	use a fairly clever algorithm to determine tree positions.	352	use a fairly clever algorithm to determine tree positions.
353	see http://iland.boku.ac.at/initialize+trees	353	see http://iland.boku.ac.at/initialize+trees
354	@param content tree init file (including headers) in a string	354	@param content tree init file (including headers) in a string
355	@param ru resource unit	355	@param ru resource unit
356	@param fileName source file name (for error reporting)	356	@param fileName source file name (for error reporting)
357	@return number of trees added	357	@return number of trees added
358	*/	358	*/
359	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)	359	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
360	{	360	{
361	if (!ru)	361	if (!ru)
362	ru = mModel->ru();	362	ru = mModel->ru();
363	Q_ASSERT(ru!=0);	363	Q_ASSERT(ru!=0);
364	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	364	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
365	Q_ASSERT(speciesSet!=0);	365	Q_ASSERT(speciesSet!=0);
366		366
367	//DebugTimer t("StandLoader::loadiLandFile");	367	//DebugTimer t("StandLoader::loadiLandFile");
368	CSVFile infile;	368	CSVFile infile;
369	infile.loadFromString(content);	369	infile.loadFromString(content);
370		370
371	int icount = infile.columnIndex("count");	371	int icount = infile.columnIndex("count");
372	int ispecies = infile.columnIndex("species");	372	int ispecies = infile.columnIndex("species");
373	int idbh_from = infile.columnIndex("dbh_from");	373	int idbh_from = infile.columnIndex("dbh_from");
374	int idbh_to = infile.columnIndex("dbh_to");	374	int idbh_to = infile.columnIndex("dbh_to");
375	int ihd = infile.columnIndex("hd");	375	int ihd = infile.columnIndex("hd");
376	int iage = infile.columnIndex("age");	376	int iage = infile.columnIndex("age");
377	int idensity = infile.columnIndex("density");	377	int idensity = infile.columnIndex("density");
378	if (icount<0 \|\| ispecies<0 \|\| idbh_from<0 \|\| idbh_to<0 \|\| ihd<0 \|\| iage<0)	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));	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		380
381	mInitItems.clear();	381	mInitItems.clear();
382	InitFileItem item;	382	InitFileItem item;
383	bool ok;	383	bool ok;
384	int total_count = 0;	384	int total_count = 0;
385	for (int row=0;row<infile.rowCount();row++) {	385	for (int row=0;row<infile.rowCount();row++) {
386	item.count = infile.value(row, icount).toInt();	386	item.count = infile.value(row, icount).toInt();
387	total_count += item.count;	387	total_count += item.count;
388	item.dbh_from = infile.value(row, idbh_from).toDouble();	388	item.dbh_from = infile.value(row, idbh_from).toDouble();
389	item.dbh_to = infile.value(row, idbh_to).toDouble();	389	item.dbh_to = infile.value(row, idbh_to).toDouble();
390	item.hd = infile.value(row, ihd).toDouble();	390	item.hd = infile.value(row, ihd).toDouble();
391	if (item.hd==0. \|\| item.dbh_from / 100. * item.hd < 4.)	391	if (item.hd==0. \|\| item.dbh_from / 100. * item.hd < 4.)
392	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) );	-
-		392	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) ;
-		393	//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) );
393	ok = true;	394	ok = true;
394	if (iage>=0)	395	if (iage>=0)
395	item.age = infile.value(row, iage).toInt(&ok);	396	item.age = infile.value(row, iage).toInt(&ok);
396	if (iage<0 \|\| !ok)	397	if (iage<0 \|\| !ok)
397	item.age = 0;	398	item.age = 0;
398		399
399	item.species = speciesSet->species(infile.value(row, ispecies).toString());	400	item.species = speciesSet->species(infile.value(row, ispecies).toString());
400	if (idensity>=0)	401	if (idensity>=0)
401	item.density = infile.value(row, idensity).toDouble();	402	item.density = infile.value(row, idensity).toDouble();
402	else	403	else
403	item.density = 0.;	404	item.density = 0.;
404	if (item.density<-1 \|\| item.density>1)	405	if (item.density<-1 \|\| item.density>1)
405	throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")	406	throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
406	.arg(item.density)	407	.arg(item.density)
407	.arg(fileName)	408	.arg(fileName)
408	.arg(row));	409	.arg(row));
409	if (!item.species) {	410	if (!item.species) {
410	throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")	411	throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
411	.arg(infile.value(row, ispecies).toString())	412	.arg(infile.value(row, ispecies).toString())
412	.arg(fileName)	413	.arg(fileName)
413	.arg(row));	414	.arg(row));
414	}	415	}
415	mInitItems.push_back(item);	416	mInitItems.push_back(item);
416	}	417	}
417	// setup the random distribution	418	// setup the random distribution
418	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");	419	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
419	if (logLevelInfo()) qDebug() << "density function:" << density_func;	420	if (logLevelInfo()) qDebug() << "density function:" << density_func;
420	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {	421	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {
421	if (mRandom)	422	if (mRandom)
422	delete mRandom;	423	delete mRandom;
423	mRandom=new RandomCustomPDF(density_func);	424	mRandom=new RandomCustomPDF(density_func);
424	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;	425	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
425	}	426	}
426	if (stand_id>0) {	427	if (stand_id>0) {
427	// execute stand based initialization	428	// execute stand based initialization
428	executeiLandInitStand(stand_id);	429	executeiLandInitStand(stand_id);
429	} else {	430	} else {
430	// exeucte the initialization based on single resource units	431	// exeucte the initialization based on single resource units
431	executeiLandInit(ru);	432	executeiLandInit(ru);
432	ru->cleanTreeList();	433	ru->cleanTreeList();
433	}	434	}
434	return total_count;	435	return total_count;
435		436
436	}	437	}
437		438
438	int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)	439	int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)
439	{	440	{
440	if (!QFile::exists(fileName))	441	if (!QFile::exists(fileName))
441	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));	442	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
442	QString content = Helper::loadTextFile(fileName);	443	QString content = Helper::loadTextFile(fileName);
443	return loadDistributionList(content, ru, stand_id, fileName);	444	return loadDistributionList(content, ru, stand_id, fileName);
444	}	445	}
445		446
446	// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.	447	// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
447	// e.g. 12 = centerpixel, 0: upper left corner, ...	448	// e.g. 12 = centerpixel, 0: upper left corner, ...
448	int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,	449	int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
449	1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};	450	1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
450	int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,	451	int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
451	6,18,2,10,4,24,12,0,8,14,20,22};	452	6,18,2,10,4,24,12,0,8,14,20,22};
452		453
453		454
454		455
455	// sort function	456	// sort function
456	bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)	457	bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
457	{	458	{
458	return s1.second < s2.second;	459	return s1.second < s2.second;
459	}	460	}
460		461
461	struct SInitPixel {	462	struct SInitPixel {
462	double basal_area; // accumulated basal area	463	double basal_area; // accumulated basal area
463	QPoint pixelOffset; // location of the pixel	464	QPoint pixelOffset; // location of the pixel
464	ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to	465	ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to
465	double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)	466	double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)
466	SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}	467	SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}
467	};	468	};
468		469
469	bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)	470	bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)
470	{	471	{
471	return s1.basal_area < s2.basal_area;	472	return s1.basal_area < s2.basal_area;
472	}	473	}
473		474
474		475
475	/**	476	/**
476	*/	477	*/
477		478
478	void StandLoader::executeiLandInit(ResourceUnit *ru)	479	void StandLoader::executeiLandInit(ResourceUnit *ru)
479	{	480	{
480		481
481	QPointF offset = ru->boundingBox().topLeft();	482	QPointF offset = ru->boundingBox().topLeft();
482	QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);	483	QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);
483		484
484	// a multimap holds a list for all trees.	485	// a multimap holds a list for all trees.
485	// key is the index of a 10x10m pixel within the resource unit	486	// key is the index of a 10x10m pixel within the resource unit
486	QMultiMap<int, int> tree_map;	487	QMultiMap<int, int> tree_map;
487	//QHash<int,SInitPixel> tcount;	488	//QHash<int,SInitPixel> tcount;
488		489
489	QVector<QPair<int, double> > tcount; // counts	490	QVector<QPair<int, double> > tcount; // counts
490	for (int i=0;i<100;i++)	491	for (int i=0;i<100;i++)
491	tcount.push_back(QPair<int,double>(i,0.));	492	tcount.push_back(QPair<int,double>(i,0.));
492		493
493	int key;	494	int key;
494	double rand_val, rand_fraction;	495	double rand_val, rand_fraction;
495	int total_count = 0;	496	int total_count = 0;
496	foreach(const InitFileItem &item, mInitItems) {	497	foreach(const InitFileItem &item, mInitItems) {
497	rand_fraction = fabs(double(item.density));	498	rand_fraction = fabs(double(item.density));
498	for (int i=0;i<item.count;i++) {	499	for (int i=0;i<item.count;i++) {

Subversion Repositories public iLand

(root)/src/core/standloader.cpp @ 1222 - Rev 779 → 782