WebSVN - public iLand - Diff - Rev 739 and 743


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

Subversion Repositories public iLand

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