WebSVN - public iLand - Diff - Rev 736 and 739


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));
        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;
            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();
        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) {
        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)
{
    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 736	Rev 739
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()+(*p)->boundingBox().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	height_grid.swap(p);	101	height_grid.swap(p);
102	mInitHeightGrid = height_grid.data();	102	mInitHeightGrid = height_grid.data();
103		103
104	QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");	104	QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");
105	mHeightGridResponse = new Expression(expr);	105	mHeightGridResponse = new Expression(expr);
106	mHeightGridResponse->setLinearizationEnabled(true);	106	mHeightGridResponse->setLinearizationEnabled(true);
107	}	107	}
108		108
109	Tree::resetStatistics();	109	Tree::resetStatistics();
110		110
111	// one global init-file for the whole area:	111	// one global init-file for the whole area:
112	if (copy_mode=="single") {	112	if (copy_mode=="single") {
113	loadInitFile(fileName, type);	113	loadInitFile(fileName, type);
114	evaluateDebugTrees();	114	evaluateDebugTrees();
115	return;	115	return;
116	}	116	}
117		117
118		118
119	// call a single tree init for each resource unit	119	// call a single tree init for each resource unit
120	if (copy_mode=="unit") {	120	if (copy_mode=="unit") {
121	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {	121	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
122	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);	122	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);
123	// set environment	123	// set environment
124	g->model()->environment()->setPosition(ru->boundingBox().center());	124	g->model()->environment()->setPosition(ru->boundingBox().center());
125	type = xml.value("type", "");	125	type = xml.value("type", "");
126	fileName = xml.value("file", "");	126	fileName = xml.value("file", "");
127	if (fileName.isEmpty())	127	if (fileName.isEmpty())
128	continue;	128	continue;
129	loadInitFile(fileName, type, 0, ru);	129	loadInitFile(fileName, type, 0, ru);
130	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";	130	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
131	}	131	}
132	evaluateDebugTrees();	132	evaluateDebugTrees();
133	return;	133	return;
134	}	134	}
135		135
136	// map-modus: load a init file for each "polygon" in the standgrid	136	// map-modus: load a init file for each "polygon" in the standgrid
137	if (copy_mode=="map") {	137	if (copy_mode=="map") {
138	if (!g->model()->standGrid() \|\| !g->model()->standGrid()->isValid())	138	if (!g->model()->standGrid() \|\| !g->model()->standGrid()->isValid())
139	throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));	139	throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));
140	QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");	140	QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");
141		141
142	CSVFile map_file(map_file_name);	142	CSVFile map_file(map_file_name);
143	if (map_file.rowCount()==0)	143	if (map_file.rowCount()==0)
144	throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));	144	throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));
145	int ikey = map_file.columnIndex("id");	145	int ikey = map_file.columnIndex("id");
146	int ivalue = map_file.columnIndex("filename");	146	int ivalue = map_file.columnIndex("filename");
147	if (ikey<0 \|\| ivalue<0)	147	if (ikey<0 \|\| ivalue<0)
148	throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));	148	throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));
149	QString file_name;	149	QString file_name;
150	for (int i=0;i<map_file.rowCount();i++) {	150	for (int i=0;i<map_file.rowCount();i++) {
151	int key = map_file.value(i, ikey).toInt();	151	int key = map_file.value(i, ikey).toInt();
152	if (key>0) {	152	if (key>0) {
153	file_name = map_file.value(i, ivalue).toString();	153	file_name = map_file.value(i, ivalue).toString();
154	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;	154	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;
155	if (!file_name.isEmpty())	155	if (!file_name.isEmpty())
156	loadInitFile(file_name, type, key, NULL);	156	loadInitFile(file_name, type, key, NULL);
157	}	157	}
158	}	158	}
159	mInitHeightGrid = 0;	159	mInitHeightGrid = 0;
160	evaluateDebugTrees();	160	evaluateDebugTrees();
161	return;	161	return;
162	}	162	}
163	if (copy_mode=="snapshot") {	163	if (copy_mode=="snapshot") {
164	// load a snapshot from a file	164	// load a snapshot from a file
165	Snapshot shot;	165	Snapshot shot;
166		166
167	QString input_db = GlobalSettings::instance()->path(fileName);	167	QString input_db = GlobalSettings::instance()->path(fileName);
168	shot.loadSnapshot(input_db);	168	shot.loadSnapshot(input_db);
169	return;	169	return;
170	}	170	}
171	throw IException("StandLoader::processInit: invalid initalization.mode!");	171	throw IException("StandLoader::processInit: invalid initalization.mode!");
172	}	172	}
173		173
174	void StandLoader::evaluateDebugTrees()	174	void StandLoader::evaluateDebugTrees()
175	{	175	{
176	// evaluate debugging	176	// evaluate debugging
177	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");	177	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
178	int counter=0;	178	int counter=0;
179	if (!dbg_str.isEmpty()) {	179	if (!dbg_str.isEmpty()) {
180	if (dbg_str == "debugstamp") {	180	if (dbg_str == "debugstamp") {
181	qDebug() << "debug_tree = debugstamp: try touching all trees...";	181	qDebug() << "debug_tree = debugstamp: try touching all trees...";
182	// try to force an error if a stamp is invalid	182	// try to force an error if a stamp is invalid
183	AllTreeIterator at(GlobalSettings::instance()->model());	183	AllTreeIterator at(GlobalSettings::instance()->model());
184	double total_offset;	184	double total_offset;
185	while (Tree *t=at.next()) {	185	while (Tree *t=at.next()) {
186	total_offset += t->stamp()->offset();	186	total_offset += t->stamp()->offset();
187	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))	187	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
188	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";	188	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
189	}	189	}
190	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";	190	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
191	return;	191	return;
192	}	192	}
193	TreeWrapper tw;	193	TreeWrapper tw;
194	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables	194	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
195	AllTreeIterator at(GlobalSettings::instance()->model());	195	AllTreeIterator at(GlobalSettings::instance()->model());
196	double result;	196	double result;
197	while (Tree *t = at.next()) {	197	while (Tree *t = at.next()) {
198	tw.setTree(t);	198	tw.setTree(t);
199	result = dexp.execute();	199	result = dexp.execute();
200	if (result) {	200	if (result) {
201	t->enableDebugging();	201	t->enableDebugging();
202	counter++;	202	counter++;
203	}	203	}
204	}	204	}
205	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";	205	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
206	}	206	}
207	}	207	}
208		208
209	/// load a single init file. Calls loadPicusFile() or loadiLandFile()	209	/// load a single init file. Calls loadPicusFile() or loadiLandFile()
210	/// @param fileName file to load	210	/// @param fileName file to load
211	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"	211	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
212	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)	212	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
213	{	213	{
214	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");	214	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
215	if (!QFile::exists(pathFileName))	215	if (!QFile::exists(pathFileName))
216	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));	216	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
217		217
218	if (type=="picus" \|\| type=="single") {	218	if (type=="picus" \|\| type=="single") {
219	if (stand_id>0)	219	if (stand_id>0)
220	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);	220	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
221	return loadPicusFile(pathFileName, ru);	221	return loadPicusFile(pathFileName, ru);
222	}	222	}
223	if (type=="iland" \|\| type=="distribution")	223	if (type=="iland" \|\| type=="distribution")
224	return loadiLandFile(pathFileName, ru, stand_id);	224	return loadiLandFile(pathFileName, ru, stand_id);
225		225
226	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);	226	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
227	}	227	}
228		228
229	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)	229	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
230	{	230	{
231	QString content = Helper::loadTextFile(fileName);	231	QString content = Helper::loadTextFile(fileName);
232	if (content.isEmpty()) {	232	if (content.isEmpty()) {
233	qDebug() << "file not found: " + fileName;	233	qDebug() << "file not found: " + fileName;
234	return 0;	234	return 0;
235	}	235	}
236	return loadSingleTreeList(content, ru, fileName);	236	return loadSingleTreeList(content, ru, fileName);
237	}	237	}
238		238
239	/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.	239	/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.
240	returns the number of loaded trees.	240	returns the number of loaded trees.
241	*/	241	*/
242	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)	242	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
243	{	243	{
244	if (!ru)	244	if (!ru)
245	ru = mModel->ru();	245	ru = mModel->ru();
246	Q_ASSERT(ru!=0);	246	Q_ASSERT(ru!=0);
247		247
248	QPointF offset = ru->boundingBox().topLeft();	248	QPointF offset = ru->boundingBox().topLeft();
249	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	249	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
250		250
251	QString my_content(content);	251	QString my_content(content);
252	// cut out the <trees> </trees> part if present	252	// cut out the <trees> </trees> part if present
253	if (content.contains("<trees>")) {	253	if (content.contains("<trees>")) {
254	QRegExp rx(".<trees>(.)</trees>.*");	254	QRegExp rx(".<trees>(.)</trees>.*");
255	rx.indexIn(content, 0);	255	rx.indexIn(content, 0);
256	if (rx.capturedTexts().count()<1)	256	if (rx.capturedTexts().count()<1)
257	return 0;	257	return 0;
258	my_content = rx.cap(1).trimmed();	258	my_content = rx.cap(1).trimmed();
259	}	259	}
260		260
261	QStringList lines=my_content.split('\n');	261	QStringList lines=my_content.split('\n');
262	if (lines.count()<2)	262	if (lines.count()<2)
263	return 0;	263	return 0;
264	// drop comments	264	// drop comments
265	while (!lines.isEmpty() && lines.front().startsWith('#') )	265	while (!lines.isEmpty() && lines.front().startsWith('#') )
266	lines.pop_front();	266	lines.pop_front();
267	while (!lines.isEmpty() && lines.last().isEmpty())	267	while (!lines.isEmpty() && lines.last().isEmpty())
268	lines.removeLast();	268	lines.removeLast();
269		269
270	char sep='\t';	270	char sep='\t';
271	if (!lines[0].contains(sep))	271	if (!lines[0].contains(sep))
272	sep=';';	272	sep=';';
273	QStringList headers = lines[0].trimmed().split(sep);	273	QStringList headers = lines[0].trimmed().split(sep);
274		274
275	int iID = headers.indexOf("id");	275	int iID = headers.indexOf("id");
276	int iX = headers.indexOf("x");	276	int iX = headers.indexOf("x");
277	int iY = headers.indexOf("y");	277	int iY = headers.indexOf("y");
278	int iBhd = headers.indexOf("bhdfrom");	278	int iBhd = headers.indexOf("bhdfrom");
279	if (iBhd<0)	279	if (iBhd<0)
280	iBhd = headers.indexOf("dbh");	280	iBhd = headers.indexOf("dbh");
281	double height_conversion = 100.;	281	double height_conversion = 100.;
282	int iHeight = headers.indexOf("treeheight");	282	int iHeight = headers.indexOf("treeheight");
283	if (iHeight<0) {	283	if (iHeight<0) {
284	iHeight = headers.indexOf("height");	284	iHeight = headers.indexOf("height");
285	height_conversion = 1.; // in meter	285	height_conversion = 1.; // in meter
286	}	286	}
287	int iSpecies = headers.indexOf("species");	287	int iSpecies = headers.indexOf("species");
288	int iAge = headers.indexOf("age");	288	int iAge = headers.indexOf("age");
289	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)	289	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)
290	throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));	290	throw IException(QString("Initfile %1 is not valid!\nObligatory columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
291		291
292	double dbh;	292	double dbh;
293	bool ok;	293	bool ok;
294	int cnt=0;	294	int cnt=0;
295	QString speciesid;	295	QString speciesid;
296	for (int i=1;i<lines.count();i++) {	296	for (int i=1;i<lines.count();i++) {
297	QString &line = lines[i];	297	QString &line = lines[i];
298	dbh = line.section(sep, iBhd, iBhd).toDouble();	298	dbh = line.section(sep, iBhd, iBhd).toDouble();
299	if (dbh<5.)	299	if (dbh<5.)
300	continue;	300	continue;
301		301
302	QPointF f;	302	QPointF f;
303	if (iX>=0 && iY>=0) {	303	if (iX>=0 && iY>=0) {
304	f.setX( line.section(sep, iX, iX).toDouble() );	304	f.setX( line.section(sep, iX, iX).toDouble() );
305	f.setY( line.section(sep, iY, iY).toDouble() );	305	f.setY( line.section(sep, iY, iY).toDouble() );
306	f+=offset;	306	f+=offset;
307		307
308	}	308	}
309	// position valid?	309	// position valid?
310	if (!mModel->heightGrid()->valueAt(f).isValid())	310	if (!mModel->heightGrid()->valueAt(f).isValid())
311	continue;	311	continue;
312	Tree &tree = ru->newTree();	312	Tree &tree = ru->newTree();
313	tree.setPosition(f);	313	tree.setPosition(f);
314	if (iID>=0)	314	if (iID>=0)
315	tree.setId(line.section(sep, iID, iID).toInt() );	315	tree.setId(line.section(sep, iID, iID).toInt() );
316		316
317	tree.setDbh(dbh);	317	tree.setDbh(dbh);
318	tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary	318	tree.setHeight(line.section(sep, iHeight, iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
319		319
320	speciesid = line.section(sep, iSpecies, iSpecies).trimmed();	320	speciesid = line.section(sep, iSpecies, iSpecies).trimmed();
321	int picusid = speciesid.toInt(&ok);	321	int picusid = speciesid.toInt(&ok);
322	if (ok) {	322	if (ok) {
323	int idx = picusSpeciesIds.indexOf(picusid);	323	int idx = picusSpeciesIds.indexOf(picusid);
324	if (idx==-1)	324	if (idx==-1)
325	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));	325	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
326	speciesid = iLandSpeciesIds[idx];	326	speciesid = iLandSpeciesIds[idx];
327	}	327	}
328	Species *s = speciesSet->species(speciesid);	328	Species *s = speciesSet->species(speciesid);
329	if (!ru \|\| !s)	329	if (!ru \|\| !s)
330	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));	330	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
331	tree.setSpecies(s);	331	tree.setSpecies(s);
332		332
333	ok = true;	333	ok = true;
334	if (iAge>=0)	334	if (iAge>=0)
335	tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a real age	335	tree.setAge(line.section(sep, iAge, iAge).toInt(&ok), tree.height()); // this is a real age
336	if (iAge<0 \|\| !ok \|\| tree.age()==0)	336	if (iAge<0 \|\| !ok \|\| tree.age()==0)
337	tree.setAge(0, tree.height()); // no real tree age available	337	tree.setAge(0, tree.height()); // no real tree age available
338		338
339	tree.setRU(ru);	339	tree.setRU(ru);
340	tree.setup();	340	tree.setup();
341	cnt++;	341	cnt++;
342	}	342	}
343	return cnt;	343	return cnt;
344	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";	344	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
345	//qDebug() << "lines: " << lines;	345	//qDebug() << "lines: " << lines;
346	}	346	}
347		347
348	/** initialize trees on a resource unit based on dbh distributions.	348	/** initialize trees on a resource unit based on dbh distributions.
349	use a fairly clever algorithm to determine tree positions.	349	use a fairly clever algorithm to determine tree positions.
350	see http://iland.boku.ac.at/initialize+trees	350	see http://iland.boku.ac.at/initialize+trees
351	@param content tree init file (including headers) in a string	351	@param content tree init file (including headers) in a string
352	@param ru resource unit	352	@param ru resource unit
353	@param fileName source file name (for error reporting)	353	@param fileName source file name (for error reporting)
354	@return number of trees added	354	@return number of trees added
355	*/	355	*/
356	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)	356	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
357	{	357	{
358	if (!ru)	358	if (!ru)
359	ru = mModel->ru();	359	ru = mModel->ru();
360	Q_ASSERT(ru!=0);	360	Q_ASSERT(ru!=0);
361	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	361	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
362	Q_ASSERT(speciesSet!=0);	362	Q_ASSERT(speciesSet!=0);
363		363
364	//DebugTimer t("StandLoader::loadiLandFile");	364	//DebugTimer t("StandLoader::loadiLandFile");
365	CSVFile infile;	365	CSVFile infile;
366	infile.loadFromString(content);	366	infile.loadFromString(content);
367		367
368	int icount = infile.columnIndex("count");	368	int icount = infile.columnIndex("count");
369	int ispecies = infile.columnIndex("species");	369	int ispecies = infile.columnIndex("species");
370	int idbh_from = infile.columnIndex("dbh_from");	370	int idbh_from = infile.columnIndex("dbh_from");
371	int idbh_to = infile.columnIndex("dbh_to");	371	int idbh_to = infile.columnIndex("dbh_to");
372	int ihd = infile.columnIndex("hd");	372	int ihd = infile.columnIndex("hd");
373	int iage = infile.columnIndex("age");	373	int iage = infile.columnIndex("age");
374	int idensity = infile.columnIndex("density");	374	int idensity = infile.columnIndex("density");
375	if (icount<0 \|\| ispecies<0 \|\| idbh_from<0 \|\| idbh_to<0 \|\| ihd<0 \|\| iage<0)	375	if (icount<0 \|\| ispecies<0 \|\| idbh_from<0 \|\| idbh_to<0 \|\| ihd<0 \|\| iage<0)
376	throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));	376	throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));
377		377
378	mInitItems.clear();	378	mInitItems.clear();
379	InitFileItem item;	379	InitFileItem item;
380	bool ok;	380	bool ok;
381	int total_count = 0;	381	int total_count = 0;
382	for (int row=0;row<infile.rowCount();row++) {	382	for (int row=0;row<infile.rowCount();row++) {
383	item.count = infile.value(row, icount).toInt();	383	item.count = infile.value(row, icount).toInt();
384	total_count += item.count;	384	total_count += item.count;
385	item.dbh_from = infile.value(row, idbh_from).toDouble();	385	item.dbh_from = infile.value(row, idbh_from).toDouble();
386	item.dbh_to = infile.value(row, idbh_to).toDouble();	386	item.dbh_to = infile.value(row, idbh_to).toDouble();
387	item.hd = infile.value(row, ihd).toDouble();	387	item.hd = infile.value(row, ihd).toDouble();
388	ok = true;	388	ok = true;
389	if (iage>=0)	389	if (iage>=0)
390	item.age = infile.value(row, iage).toInt(&ok);	390	item.age = infile.value(row, iage).toInt(&ok);
391	if (iage<0 \|\| !ok)	391	if (iage<0 \|\| !ok)
392	item.age = 0;	392	item.age = 0;
393		393
394	item.species = speciesSet->species(infile.value(row, ispecies).toString());	394	item.species = speciesSet->species(infile.value(row, ispecies).toString());
395	if (idensity>=0)	395	if (idensity>=0)
396	item.density = infile.value(row, idensity).toDouble();	396	item.density = infile.value(row, idensity).toDouble();
397	else	397	else
398	item.density = 0.;	398	item.density = 0.;
399	if (item.density<-1 \|\| item.density>1)	399	if (item.density<-1 \|\| item.density>1)
400	throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")	400	throw IException(QString("load-ini-file: invalid value for density. Allowed range is -1..1: '%1' in file '%2', line %3.")
401	.arg(item.density)	401	.arg(item.density)
402	.arg(fileName)	402	.arg(fileName)
403	.arg(row));	403	.arg(row));
404	if (!item.species) {	404	if (!item.species) {
405	throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")	405	throw IException(QString("load-ini-file: unknown speices '%1' in file '%2', line %3.")
406	.arg(infile.value(row, ispecies).toString())	406	.arg(infile.value(row, ispecies).toString())
407	.arg(fileName)	407	.arg(fileName)
408	.arg(row));	408	.arg(row));
409	}	409	}
410	mInitItems.push_back(item);	410	mInitItems.push_back(item);
411	}	411	}
412	// setup the random distribution	412	// setup the random distribution
413	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");	413	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
414	if (logLevelInfo()) qDebug() << "density function:" << density_func;	414	if (logLevelInfo()) qDebug() << "density function:" << density_func;
415	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {	415	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {
416	if (mRandom)	416	if (mRandom)
417	delete mRandom;	417	delete mRandom;
418	mRandom=new RandomCustomPDF(density_func);	418	mRandom=new RandomCustomPDF(density_func);
419	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;	419	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
420	}	420	}
421	if (stand_id>0) {	421	if (stand_id>0) {
422	// execute stand based initialization	422	// execute stand based initialization
423	executeiLandInitStand(stand_id);	423	executeiLandInitStand(stand_id);
424	} else {	424	} else {
425	// exeucte the initialization based on single resource units	425	// exeucte the initialization based on single resource units
426	executeiLandInit(ru);	426	executeiLandInit(ru);
427	ru->cleanTreeList();	427	ru->cleanTreeList();
428	}	428	}
429	return total_count;	429	return total_count;
430		430
431	}	431	}
432		432
433	int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)	433	int StandLoader::loadiLandFile(const QString &fileName, ResourceUnit *ru, int stand_id)
434	{	434	{
435	if (!QFile::exists(fileName))	435	if (!QFile::exists(fileName))
436	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));	436	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
437	QString content = Helper::loadTextFile(fileName);	437	QString content = Helper::loadTextFile(fileName);
438	return loadDistributionList(content, ru, stand_id, fileName);	438	return loadDistributionList(content, ru, stand_id, fileName);
439	}	439	}
440		440
441	// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.	441	// evenlist: tentative order of pixel-indices (within a 5x5 grid) used as tree positions.
442	// e.g. 12 = centerpixel, 0: upper left corner, ...	442	// e.g. 12 = centerpixel, 0: upper left corner, ...
443	int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,	443	int evenlist[25] = { 12, 6, 18, 16, 8, 22, 2, 10, 14, 0, 24, 20, 4,
444	1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};	444	1, 13, 15, 19, 21, 3, 7, 11, 17, 23, 5, 9};
445	int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,	445	int unevenlist[25] = { 11,13,7,17, 1,19,5,21, 9,23,3,15,
446	6,18,2,10,4,24,12,0,8,14,20,22};	446	6,18,2,10,4,24,12,0,8,14,20,22};
447		447
448		448
449		449
450	// sort function	450	// sort function
451	bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)	451	bool sortPairLessThan(const QPair<int, double> &s1, const QPair<int, double> &s2)
452	{	452	{
453	return s1.second < s2.second;	453	return s1.second < s2.second;
454	}	454	}
455		455
456	struct SInitPixel {	456	struct SInitPixel {
457	double basal_area; // accumulated basal area	457	double basal_area; // accumulated basal area
458	QPoint pixelOffset; // location of the pixel	458	QPoint pixelOffset; // location of the pixel
459	ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to	459	ResourceUnit *resource_unit; // pointer to the resource unit the pixel belongs to
460	double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)	460	double h_max; // predefined maximum height at given pixel (if available from LIDAR or so)
461	SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}	461	SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.) {}
462	};	462	};
463		463
464	bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)	464	bool sortInitPixelLessThan(const SInitPixel &s1, const SInitPixel &s2)
465	{	465	{
466	return s1.basal_area < s2.basal_area;	466	return s1.basal_area < s2.basal_area;
467	}	467	}
468		468
469		469
470	/**	470	/**
471	*/	471	*/
472		472
473	void StandLoader::executeiLandInit(ResourceUnit *ru)	473	void StandLoader::executeiLandInit(ResourceUnit *ru)
474	{	474	{
475		475
476	QPointF offset = ru->boundingBox().topLeft();	476	QPointF offset = ru->boundingBox().topLeft();
477	QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);	477	QPoint offsetIdx = GlobalSettings::instance()->model()->grid()->indexAt(offset);
478		478
479	// a multimap holds a list for all trees.	479	// a multimap holds a list for all trees.
480	// key is the index of a 10x10m pixel within the resource unit	480	// key is the index of a 10x10m pixel within the resource unit
481	QMultiMap<int, int> tree_map;	481	QMultiMap<int, int> tree_map;
482	//QHash<int,SInitPixel> tcount;	482	//QHash<int,SInitPixel> tcount;
483		483
484	QVector<QPair<int, double> > tcount; // counts	484	QVector<QPair<int, double> > tcount; // counts
485	for (int i=0;i<100;i++)	485	for (int i=0;i<100;i++)
486	tcount.push_back(QPair<int,double>(i,0.));	486	tcount.push_back(QPair<int,double>(i,0.));
487		487
488	int key;	488	int key;
489	double rand_val, rand_fraction;	489	double rand_val, rand_fraction;
490	int total_count = 0;	490	int total_count = 0;
491	foreach(const InitFileItem &item, mInitItems) {	491	foreach(const InitFileItem &item, mInitItems) {
492	rand_fraction = fabs(double(item.density));	492	rand_fraction = fabs(double(item.density));
493	for (int i=0;i<item.count;i++) {	493	for (int i=0;i<item.count;i++) {
494	// create trees	494	// create trees
495	int tree_idx = ru->newTreeIndex();	495	int tree_idx = ru->newTreeIndex();
496	Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree	496	Tree &tree = ru->trees()[tree_idx]; // get reference to modify tree
497	tree.setDbh(nrandom(item.dbh_from, item.dbh_to));	497	tree.setDbh(nrandom(item.dbh_from, item.dbh_to));
498	tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height	498	tree.setHeight(tree.dbh()/100. * item.hd); // dbh from cm->m, *hd-ratio -> meter height
499	tree.setSpecies(item.species);	499	tree.setSpecies(item.species);

Subversion Repositories public iLand

(root)/src/core/standloader.cpp @ 1222 - Rev 736 → 739