WebSVN - public iLand - Diff - Rev 1190 and 1203


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 "species.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"
#include "grasscover.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
static QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
static 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.valueInt("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->linearize(0., 2.);
    }

    Tree::resetStatistics();

    // one global init-file for the whole area:
    if (copy_mode=="single") {
        // useful for 1ha simulations only...
        if (GlobalSettings::instance()->model()->ruList().size()>1)
            throw IException("Error initialization: 'mode' is 'single' but more than one resource unit is simulated (consider using another 'mode').");

        loadInitFile(fileName, type, 0, GlobalSettings::instance()->model()->ru()); // this is the first resource unit
        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;
    }

    // standgrid mode: load one large init file
    if (copy_mode=="standgrid") {
        fileName = GlobalSettings::instance()->path(fileName, "init");
        if (!QFile::exists(fileName))
            throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
        QString content = Helper::loadTextFile(fileName);
        // this processes the init file (also does the checking) and
        // stores in a QHash datastrucutre
        parseInitFile(content, fileName);

        // setup the random distribution
        QString density_func = xml.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 (mStandInitItems.isEmpty()) {
            qDebug() << "Initialize trees ('standgrid'-mode): no items to process (empty landscape).";
            return;
            //throw IException("StandLoader::processInit: 'mode' is 'standgrid' but the init file is either empty or contains no 'stand_id'-column.");
        }
        QHash<int, QVector<InitFileItem> >::const_iterator it = mStandInitItems.constBegin();
        while (it!=mStandInitItems.constEnd()) {
            mInitItems = it.value(); // copy the items...
            executeiLandInitStand(it.key());
            ++it;
        }
        qDebug() << "finished setup of trees.";
        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::processAfterInit()
{
    XmlHelper xml(GlobalSettings::instance()->settings().node("model.initialization"));

    QString mode = xml.value("mode", "copy");
    if (mode=="standgrid") {
        // load a file with saplings per polygon
        QString  filename = xml.value("saplingFile", "");
        if (filename.isEmpty())
            return;
        filename = GlobalSettings::instance()->path(filename, "init");
        if (!QFile::exists(filename))
            throw IException(QString("load-sapling-ini-file: file '%1' does not exist.").arg(filename));
        CSVFile init_file(filename);
        int istandid = init_file.columnIndex("stand_id");
        if (istandid==-1)
            throw IException("Sapling-Init: the init file contains no 'stand_id' column (required in 'standgrid' mode).");

        int stand_id = -99999;
        int ilow = -1, ihigh = 0;
        int total = 0;
        for (int i=0;i<init_file.rowCount();++i) {
            int row_stand = init_file.value(i, istandid).toInt();
            if (row_stand != stand_id) {
                if (stand_id>=0) {
                    // process stand
                    ihigh = i-1; // up to the last
                    total += loadSaplingsLIF(stand_id, init_file, ilow, ihigh);
                }
                ilow = i; // mark beginning of new stand
                stand_id = row_stand;
            }
        }
        if (stand_id>=0)
            total += loadSaplingsLIF(stand_id, init_file, ilow, init_file.rowCount()-1); // the last stand
        qDebug() << "initialization of sapling: total created:" << total;

    }

}

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


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

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

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

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

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

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

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

    CSVFile infile;
    infile.loadFromString(my_content);


    int iID =  infile.columnIndex("id");
    int iX = infile.columnIndex("x");
    int iY = infile.columnIndex("y");
    int iBhd = infile.columnIndex("bhdfrom");
    if (iBhd<0)
        iBhd = infile.columnIndex("dbh");
    double height_conversion = 100.;
    int iHeight =infile.columnIndex("treeheight");
    if (iHeight<0) {
        iHeight = infile.columnIndex("height");
        height_conversion = 1.; // in meter
    }
    int iSpecies = infile.columnIndex("species");
    int iAge = infile.columnIndex("age");
    if (iX==-1 || iY==-1 || iBhd==-1 || iSpecies==-1 || iHeight==-1)
        throw IException(QString("Initfile %1 is not valid!\nRequired 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<infile.rowCount();i++) {
        dbh = infile.value(i, iBhd).toDouble();

        //if (dbh<5.)
        //    continue;

        QPointF f;
        if (iX>=0 && iY>=0) {
           f.setX( infile.value(i, iX).toDouble() );
           f.setY( infile.value(i, 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(infile.value(i, iID).toInt());

        tree.setDbh(dbh);
        tree.setHeight(infile.value(i,iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary

        speciesid = infile.value(i,iSpecies).toString();
        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(infile.value(i, 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)
{
    int total_count = parseInitFile(content, fileName, ru);
    if (total_count==0)
        return 0;


    // 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::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru)
{
    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));

    int istandid = infile.columnIndex("stand_id");
    mInitItems.clear();
    mStandInitItems.clear();

    InitFileItem item;
    bool ok;
    int total_count = 0;
    for (int row=0;row<infile.rowCount();row++) {
        item.count = infile.value(row, icount).toDouble();
        total_count += item.count;
        item.dbh_from = infile.value(row, idbh_from).toDouble();
        item.dbh_to = infile.value(row, idbh_to).toDouble();
        item.hd = infile.value(row, ihd).toDouble();
        if (item.hd==0. || item.dbh_from / 100. * item.hd < 4.)
            qWarning() << QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ;
            //throw IException(QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) );
        ok = true;
        if (iage>=0)
            item.age = infile.value(row, iage).toInt(&ok);
        if (iage<0 || !ok)
            item.age = 0;

        item.species = speciesSet->species(infile.value(row, ispecies).toString());
        if (idensity>=0)
            item.density = infile.value(row, idensity).toDouble();
        else
            item.density = 0.;
        if (item.density<-1)
            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));
        }
        if (istandid>=0) {
            int standid = infile.value(row,istandid).toInt();
            mStandInitItems[standid].push_back(item);
        } else {
            mInitItems.push_back(item);
        }
    }
    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)
    bool locked; // pixel is dedicated to a single species
    SInitPixel(): basal_area(0.), resource_unit(0), h_max(-1.), locked(false) {}
};

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

bool sortInitPixelUnlocked(const SInitPixel &s1, const SInitPixel &s2)
{
    return !s1.locked && s2.locked;
}

/**
*/

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);
        }
    }
}



// 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();
    if (mCurrentMap)
        grid = mCurrentMap;

    // 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) / cRUArea;

    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.");

    Species *last_locked_species=0;
    foreach(const InitFileItem &item, mInitItems) {
        if (item.density>1.) {
            // special case with single-species-area
            if (total_count==0) {
                // randomize the pixels
                for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it)
                    it->basal_area = drandom();
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
                for (QList<SInitPixel>::iterator it=pixel_list.begin();it!=pixel_list.end();++it)
                    it->basal_area = 0.;
            }

            if (item.species != last_locked_species) {
                last_locked_species=item.species;
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelUnlocked);
            }
        } else {
            qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
            last_locked_species=0;
        }
        rand_fraction = 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.
                if (item.density <= 1.) {
                    rand_val = mRandom->get();
                    if (item.density<0)
                        rand_val = 1. - rand_val;
                    rand_val = rand_val * rand_fraction + drandom()*(1.-rand_fraction);
                } else {
                    // limited area: limit potential area using the "density" input parameter
                    rand_val = drandom() * qMin(item.density/100., 1.);
                }
                ++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 (!last_locked_species && pixel_list[key].locked)
                    found = false;
            }
            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
            if (last_locked_species)
                pixel_list[key].locked = true;

            // resort list
            if (last_locked_species==0 && ((total_count < 20 && i%2==0)
                || (total_count<100 && i%10==0 )
                || (i%30==0)) ) {
                qSort(pixel_list.begin(), pixel_list.end(), sortInitPixelLessThan);
            }
        }
    }
    if (total_misses>0 || total_tries > total_count) {
        if (logLevelInfo()) qDebug() << "init for stand" << stand_id << "treecount:" << total_count << ", tries:" << total_tries << ", misses:" << total_misses << ", %miss:" << qRound(total_misses*100 / (double)total_count);
    }

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

    foreach(const SInitPixel &p, pixel_list) {
        trees = tree_map.values(p.pixelOffset);
        c = trees.count();
        bits = 0;
        index = -1;
        double r;
        foreach(int tree_idx, trees) {
            if (c>18) {
                index = (index + 1)%25;
            } else {
                int stop=1000;
                index = 0;
                do {
                    // search a random position
                    r = drandom();
                    index = limit(int(25 *  r*r), 0, 24); // use rnd()^2 to search for locations -> higher number of low indices (i.e. 50% of lookups in first 25% of locations)
                } while (isBitSet(bits, index)==true && stop--);
                if (!stop)
                    qDebug() << "executeiLandInit: found no free bit.";
                setBit(bits, index, true); // mark position as used
            }
            // get position from fixed lists (one for even, one for uneven resource units)
            pos = p.resource_unit->index()%2?evenlist[index]:unevenlist[index];
            tree_pos = p.pixelOffset * cPxPerHeight; // convert to LIF index
            tree_pos += QPoint(pos/cPxPerHeight, pos%cPxPerHeight);

            p.resource_unit->trees()[tree_idx].setPosition(tree_pos);
            // test if tree position is valid..
            if (!GlobalSettings::instance()->model()->grid()->isIndexValid(tree_pos))
                qDebug() << "Standloader: invalid position!";
        }
    }
    if (logLevelInfo())
        qDebug() << "init for stand" << stand_id << "with area" << "area (m2)" << grid->area(stand_id) << "count of 10m pixels:"  << indices.count() << "initialized trees:" << total_count;

}

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

    QList<int> indices = stand_grid->gridIndices(stand_id); // list of 10x10m pixels
    if (indices.isEmpty()) {
        qDebug() << "stand" << stand_id << "not in project area. No init performed.";
        return -1;
    }
    double area_factor = stand_grid->area(stand_id) / cRUArea; // 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 = qRound(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());
           QPoint offset=stand_grid->grid().indexOf(indices[rnd_index]);
 
           //
           offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
           int in_p = irandom(0, cPxPerHeight*cPxPerHeight); // index of lif-pixel
           offset += QPoint(in_p / cPxPerHeight, in_p % cPxPerHeight);
           SaplingCell *sc = GlobalSettings::instance()->model()->saplings()->cell(offset);
           if (sc && sc->max_height()>height) {
           //if (!ru || ru->saplingHeightForInit(offset) > height) {
               misses++;
           } else {
               // ok
               hits++;
               if (sc)
                   sc->addSapling(height, age, species->index());
               //ru->resourceUnitSpecies(species).changeSapling().addSapling(offset, height, age);
           }
           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;
}

bool LIFValueHigher(const float *a, const float *b)
{
    return *a > *b;
}

int StandLoader::loadSaplingsLIF(int stand_id, const CSVFile &init, int low_index, int high_index)
{
    const MapGrid *stand_grid;
    if (mCurrentMap)
        stand_grid = mCurrentMap; // if set
    else
        stand_grid = GlobalSettings::instance()->model()->standGrid(); // default

    if (!stand_grid->isValid(stand_id))
        return 0;

    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 0;
    }
    // prepare space for LIF-pointers (2m Pixel)
    QVector<float*> lif_ptrs;
    lif_ptrs.reserve(indices.size() * cPxPerHeight * cPxPerHeight);
    for (int l=0;l<indices.size();++l){
        QPoint offset=stand_grid->grid().indexOf(indices[l]);
        offset = offset * cPxPerHeight; // index of 10m patch -> to lif pixel coordinates
        for (int y=0;y<cPxPerHeight;++y)
            for(int x=0;x<cPxPerHeight;++x)
                lif_ptrs.push_back( GlobalSettings::instance()->model()->grid()->ptr(offset.x()+x, offset.y()+y) );
    }
    // sort based on LIF-Value
    std::sort(lif_ptrs.begin(), lif_ptrs.end(), LIFValueHigher); // higher: highest values first


    double area_factor = stand_grid->area(stand_id) / cRUArea; // multiplier for grid (e.g. 2 if stand has area of 2 hectare)

    // parse the content of the init-file
    // species
    int ispecies = init.columnIndex("species");
    int icount = init.columnIndex("count");
    int iheight = init.columnIndex("height");
    int iheightfrom = init.columnIndex("height_from");
    int iheightto = init.columnIndex("height_to");
    int iage = init.columnIndex("age");
    int itopage = init.columnIndex("age4m");
    int iminlif = init.columnIndex("min_lif");
    if ((iheightfrom==-1) ^ (iheightto==-1))
        throw IException("Error while loading saplings: height not correctly provided. Use either 'height' or 'height_from' and 'height_to'.");
    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=low_index;row<=high_index;++row) {
        double pxcount = init.value(row, icount).toDouble() * area_factor; // 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();
        double age4m = itopage==-1?10:init.value(row, itopage).toDouble();
        double height_from = iheightfrom==-1?-1.: init.value(row, iheightfrom).toDouble();
        double height_to = iheightto==-1?-1.: init.value(row, iheightto).toDouble();
        double min_lif = iminlif==-1?1.: init.value(row, iminlif).toDouble();
        // find LIF-level in the pixels
        int min_lif_index = 0;
        if (min_lif < 1.) {
            for (QVector<float*>::ConstIterator it=lif_ptrs.constBegin(); it!=lif_ptrs.constEnd(); ++it, ++min_lif_index)
                if (**it <= min_lif)
                    break;
            if (pxcount < min_lif_index) {
                // not enough LIF pixels available
                min_lif_index = static_cast<int>(pxcount); // try the brightest pixels (ie with the largest value for the LIF)
            }
        } else {
            // No LIF threshold: the full range of pixels is valid
            min_lif_index = lif_ptrs.size();
        }



        double hits = 0.;
        while (hits < pxcount) {
            int rnd_index = irandom(0, min_lif_index);
            if (iheightfrom!=-1) {
                height = limit(nrandom(height_from, height_to), 0.05,4.);
                if (age<=1.)
                    age = qMax(qRound(height/4. * age4m),1); // assume a linear relationship between height and age
            }
            QPoint offset = GlobalSettings::instance()->model()->grid()->indexOf(lif_ptrs[rnd_index]);
            ResourceUnit *ru;
            SaplingCell *sc = GlobalSettings::instance()->model()->saplings()->cell(offset, true, &ru);
            if (sc) {
                if (SaplingTree *st=sc->addSapling(static_cast<float>(height), static_cast<int>(age), species->index()))
                    hits+=std::max(1., ru->resourceUnitSpecies(st->species_index)->species()->saplingGrowthParameters().representedStemNumberH(st->height));
                else
                    hits++;
            } else {
                hits++; // avoid an infinite loop
            }


        }
        total += pxcount;

    }

    // initialize grass cover
    if (init.columnIndex("grass_cover")>-1) {
        int grass_cover_value = init.value(low_index, "grass_cover").toInt();
        if (grass_cover_value<0 || grass_cover_value>100)
            throw IException(QString("The grass cover percentage (column 'grass_cover') for stand '%1' is '%2', which is invalid (expected: 0-100)").arg(stand_id).arg(grass_cover_value));
        GlobalSettings::instance()->model()->grassCover()->setInitialValues(lif_ptrs, grass_cover_value);
    }


    return total;
}
 

Rev 1190	Rev 1203
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	#include "species.h"	29	#include "species.h"
30		30
31	#include "helper.h"	31	#include "helper.h"
32	#include "random.h"	32	#include "random.h"
33	#include "expression.h"	33	#include "expression.h"
34	#include "expressionwrapper.h"	34	#include "expressionwrapper.h"
35	#include "environment.h"	35	#include "environment.h"
36	#include "csvfile.h"	36	#include "csvfile.h"
37	#include "mapgrid.h"	37	#include "mapgrid.h"
38	#include "snapshot.h"	38	#include "snapshot.h"
39	#include "grasscover.h"	39	#include "grasscover.h"
40		40
41	/** @class StandLoader	41	/** @class StandLoader
42	@ingroup tools	42	@ingroup tools
43	loads (initializes) trees for a "stand" from various sources.	43	loads (initializes) trees for a "stand" from various sources.
44	StandLoader initializes trees on the landscape. It reads (usually) from text files, creates the	44	StandLoader initializes trees on the landscape. It reads (usually) from text files, creates the
45	trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid).	45	trees and distributes the trees on the landscape (on the ResoureceUnit or on a stand defined by a grid).
46		46
47	See http://iland.boku.ac.at/initialize+trees	47	See http://iland.boku.ac.at/initialize+trees
48	*/	48	*/
49	// provide a mapping between "Picus"-style and "iLand"-style species Ids	49	// provide a mapping between "Picus"-style and "iLand"-style species Ids
50	static QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;	50	static QVector<int> picusSpeciesIds = QVector<int>() << 0 << 1 << 17;
51	static QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";	51	static QStringList iLandSpeciesIds = QStringList() << "piab" << "piab" << "fasy";
52		52
53	StandLoader::~StandLoader()	53	StandLoader::~StandLoader()
54	{	54	{
55	if (mRandom)	55	if (mRandom)
56	delete mRandom;	56	delete mRandom;
57	if (mHeightGridResponse)	57	if (mHeightGridResponse)
58	delete mHeightGridResponse;	58	delete mHeightGridResponse;
59	}	59	}
60		60
61		61
62	void StandLoader::copyTrees()	62	void StandLoader::copyTrees()
63	{	63	{
64	// we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards	64	// we assume that all stands are equal, so wie simply COPY the trees and modify them afterwards
65	const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();	65	const Grid<ResourceUnit*> &ruGrid=mModel->RUgrid();
66	ResourceUnit **p = ruGrid.begin();	66	ResourceUnit **p = ruGrid.begin();
67	if (!p)	67	if (!p)
68	throw IException("Standloader: invalid resource unit pointer!");	68	throw IException("Standloader: invalid resource unit pointer!");
69	++p; // skip the first...	69	++p; // skip the first...
70	const QVector<Tree> &tocopy = mModel->ru()->trees();	70	const QVector<Tree> &tocopy = mModel->ru()->trees();
71	for (; p!=ruGrid.end(); ++p) {	71	for (; p!=ruGrid.end(); ++p) {
72	QRectF rect = (*p)->boundingBox();	72	QRectF rect = (*p)->boundingBox();
73	foreach(const Tree& tree, tocopy) {	73	foreach(const Tree& tree, tocopy) {
74	Tree &newtree = (*p)->newTree();	74	Tree &newtree = (*p)->newTree();
75	newtree = tree; // copy tree data...	75	newtree = tree; // copy tree data...
76	newtree.setPosition(tree.position()+rect.topLeft());	76	newtree.setPosition(tree.position()+rect.topLeft());
77	newtree.setRU(*p);	77	newtree.setRU(*p);
78	newtree.setNewId();	78	newtree.setNewId();
79	}	79	}
80	}	80	}
81	if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";	81	if (logLevelInfo()) qDebug() << Tree::statCreated() << "trees loaded / copied.";
82	}	82	}
83		83
84	/** main routine of the stand setup.	84	/** main routine of the stand setup.
85	*/	85	*/
86	void StandLoader::processInit()	86	void StandLoader::processInit()
87	{	87	{
88	GlobalSettings *g = GlobalSettings::instance();	88	GlobalSettings *g = GlobalSettings::instance();
89	XmlHelper xml(g->settings().node("model.initialization"));	89	XmlHelper xml(g->settings().node("model.initialization"));
90		90
91	QString copy_mode = xml.value("mode", "copy");	91	QString copy_mode = xml.value("mode", "copy");
92	QString type = xml.value("type", "");	92	QString type = xml.value("type", "");
93	QString fileName = xml.value("file", "");	93	QString fileName = xml.value("file", "");
94		94
95	bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);	95	bool height_grid_enabled = xml.valueBool("heightGrid.enabled", false);
96	mHeightGridTries = xml.valueInt("heightGrid.maxTries", 10);	96	mHeightGridTries = xml.valueInt("heightGrid.maxTries", 10);
97	QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function	97	QScopedPointer<const MapGrid> height_grid; // use a QScopedPointer to guarantee that the resource is freed at the end of the processInit() function
98	if (height_grid_enabled) {	98	if (height_grid_enabled) {
99	QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init");	99	QString init_height_grid_file = GlobalSettings::instance()->path(xml.value("heightGrid.fileName"), "init");
100	qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file;	100	qDebug() << "initialization: using predefined tree heights map" << init_height_grid_file;
101		101
102	QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));	102	QScopedPointer<const MapGrid> p(new MapGrid(init_height_grid_file, false));
103	if (!p->isValid()) {	103	if (!p->isValid()) {
104	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));	104	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));
105	}	105	}
106	height_grid.swap(p);	106	height_grid.swap(p);
107	mInitHeightGrid = height_grid.data();	107	mInitHeightGrid = height_grid.data();
108		108
109	QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");	109	QString expr=xml.value("heightGrid.fitFormula", "polygon(x, 0,0, 0.8,1, 1.1, 1, 1.25,0)");
110	mHeightGridResponse = new Expression(expr);	110	mHeightGridResponse = new Expression(expr);
111	mHeightGridResponse->linearize(0., 2.);	111	mHeightGridResponse->linearize(0., 2.);
112	}	112	}
113		113
114	Tree::resetStatistics();	114	Tree::resetStatistics();
115		115
116	// one global init-file for the whole area:	116	// one global init-file for the whole area:
117	if (copy_mode=="single") {	117	if (copy_mode=="single") {
118	// useful for 1ha simulations only...	118	// useful for 1ha simulations only...
119	if (GlobalSettings::instance()->model()->ruList().size()>1)	119	if (GlobalSettings::instance()->model()->ruList().size()>1)
120	throw IException("Error initialization: 'mode' is 'single' but more than one resource unit is simulated (consider using another 'mode').");	120	throw IException("Error initialization: 'mode' is 'single' but more than one resource unit is simulated (consider using another 'mode').");
121		121
122	loadInitFile(fileName, type, 0, GlobalSettings::instance()->model()->ru()); // this is the first resource unit	122	loadInitFile(fileName, type, 0, GlobalSettings::instance()->model()->ru()); // this is the first resource unit
123	evaluateDebugTrees();	123	evaluateDebugTrees();
124	return;	124	return;
125	}	125	}
126		126
127		127
128	// call a single tree init for each resource unit	128	// call a single tree init for each resource unit
129	if (copy_mode=="unit") {	129	if (copy_mode=="unit") {
130	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {	130	foreach( const ResourceUnit *const_ru, g->model()->ruList()) {
131	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);	131	ResourceUnit ru = const_cast<ResourceUnit>(const_ru);
132	// set environment	132	// set environment
133	g->model()->environment()->setPosition(ru->boundingBox().center());	133	g->model()->environment()->setPosition(ru->boundingBox().center());
134	type = xml.value("type", "");	134	type = xml.value("type", "");
135	fileName = xml.value("file", "");	135	fileName = xml.value("file", "");
136	if (fileName.isEmpty())	136	if (fileName.isEmpty())
137	continue;	137	continue;
138	loadInitFile(fileName, type, 0, ru);	138	loadInitFile(fileName, type, 0, ru);
139	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";	139	if (logLevelInfo()) qDebug() << "loaded" << fileName << "on" << ru->boundingBox() << "," << ru->trees().count() << "trees.";
140	}	140	}
141	evaluateDebugTrees();	141	evaluateDebugTrees();
142	return;	142	return;
143	}	143	}
144		144
145	// map-modus: load a init file for each "polygon" in the standgrid	145	// map-modus: load a init file for each "polygon" in the standgrid
146	if (copy_mode=="map") {	146	if (copy_mode=="map") {
147	if (!g->model()->standGrid() \|\| !g->model()->standGrid()->isValid())	147	if (!g->model()->standGrid() \|\| !g->model()->standGrid()->isValid())
148	throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));	148	throw IException(QString("Stand-Initialization: model.initialization.mode is 'map' but there is no valid stand grid defined (model.world.standGrid)"));
149	QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");	149	QString map_file_name = GlobalSettings::instance()->path(xml.value("mapFileName"), "init");
150		150
151	CSVFile map_file(map_file_name);	151	CSVFile map_file(map_file_name);
152	if (map_file.rowCount()==0)	152	if (map_file.rowCount()==0)
153	throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));	153	throw IException(QString("Stand-Initialization: the map file %1 is empty or missing!").arg(map_file_name));
154	int ikey = map_file.columnIndex("id");	154	int ikey = map_file.columnIndex("id");
155	int ivalue = map_file.columnIndex("filename");	155	int ivalue = map_file.columnIndex("filename");
156	if (ikey<0 \|\| ivalue<0)	156	if (ikey<0 \|\| ivalue<0)
157	throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));	157	throw IException(QString("Stand-Initialization: the map file %1 does not contain the mandatory columns 'id' and 'filename'!").arg(map_file_name));
158	QString file_name;	158	QString file_name;
159	for (int i=0;i<map_file.rowCount();i++) {	159	for (int i=0;i<map_file.rowCount();i++) {
160	int key = map_file.value(i, ikey).toInt();	160	int key = map_file.value(i, ikey).toInt();
161	if (key>0) {	161	if (key>0) {
162	file_name = map_file.value(i, ivalue).toString();	162	file_name = map_file.value(i, ivalue).toString();
163	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;	163	if (logLevelInfo()) qDebug() << "loading" << file_name << "for grid id" << key;
164	if (!file_name.isEmpty())	164	if (!file_name.isEmpty())
165	loadInitFile(file_name, type, key, NULL);	165	loadInitFile(file_name, type, key, NULL);
166	}	166	}
167	}	167	}
168	mInitHeightGrid = 0;	168	mInitHeightGrid = 0;
169	evaluateDebugTrees();	169	evaluateDebugTrees();
170	return;	170	return;
171	}	171	}
172		172
173	// standgrid mode: load one large init file	173	// standgrid mode: load one large init file
174	if (copy_mode=="standgrid") {	174	if (copy_mode=="standgrid") {
175	fileName = GlobalSettings::instance()->path(fileName, "init");	175	fileName = GlobalSettings::instance()->path(fileName, "init");
176	if (!QFile::exists(fileName))	176	if (!QFile::exists(fileName))
177	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));	177	throw IException(QString("load-ini-file: file '%1' does not exist.").arg(fileName));
178	QString content = Helper::loadTextFile(fileName);	178	QString content = Helper::loadTextFile(fileName);
179	// this processes the init file (also does the checking) and	179	// this processes the init file (also does the checking) and
180	// stores in a QHash datastrucutre	180	// stores in a QHash datastrucutre
181	parseInitFile(content, fileName);	181	parseInitFile(content, fileName);
182		182
183	// setup the random distribution	183	// setup the random distribution
184	QString density_func = xml.value("model.initialization.randomFunction", "1-x^2");	184	QString density_func = xml.value("model.initialization.randomFunction", "1-x^2");
185	if (logLevelInfo()) qDebug() << "density function:" << density_func;	185	if (logLevelInfo()) qDebug() << "density function:" << density_func;
186	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {	186	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {
187	if (mRandom)	187	if (mRandom)
188	delete mRandom;	188	delete mRandom;
189	mRandom=new RandomCustomPDF(density_func);	189	mRandom=new RandomCustomPDF(density_func);
190	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;	190	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
191	}	191	}
192		192
193	if (mStandInitItems.isEmpty())	-
194	throw IException("StandLoader::processInit: 'mode' is 'standgrid' but the init file is either empty or contains no 'stand_id'-column.");	-
-		193	if (mStandInitItems.isEmpty()) {
-		194	qDebug() << "Initialize trees ('standgrid'-mode): no items to process (empty landscape).";
-		195	return;
-		196	//throw IException("StandLoader::processInit: 'mode' is 'standgrid' but the init file is either empty or contains no 'stand_id'-column.");
-		197	}
195	QHash<int, QVector<InitFileItem> >::const_iterator it = mStandInitItems.constBegin();	198	QHash<int, QVector<InitFileItem> >::const_iterator it = mStandInitItems.constBegin();
196	while (it!=mStandInitItems.constEnd()) {	199	while (it!=mStandInitItems.constEnd()) {
197	mInitItems = it.value(); // copy the items...	200	mInitItems = it.value(); // copy the items...
198	executeiLandInitStand(it.key());	201	executeiLandInitStand(it.key());
199	++it;	202	++it;
200	}	203	}
201	qDebug() << "finished setup of trees.";	204	qDebug() << "finished setup of trees.";
202	evaluateDebugTrees();	205	evaluateDebugTrees();
203	return;	206	return;
204		207
205	}	208	}
206	if (copy_mode=="snapshot") {	209	if (copy_mode=="snapshot") {
207	// load a snapshot from a file	210	// load a snapshot from a file
208	Snapshot shot;	211	Snapshot shot;
209		212
210	QString input_db = GlobalSettings::instance()->path(fileName);	213	QString input_db = GlobalSettings::instance()->path(fileName);
211	shot.loadSnapshot(input_db);	214	shot.loadSnapshot(input_db);
212	return;	215	return;
213	}	216	}
214	throw IException("StandLoader::processInit: invalid initalization.mode!");	217	throw IException("StandLoader::processInit: invalid initalization.mode!");
215	}	218	}
216		219
217	void StandLoader::processAfterInit()	220	void StandLoader::processAfterInit()
218	{	221	{
219	XmlHelper xml(GlobalSettings::instance()->settings().node("model.initialization"));	222	XmlHelper xml(GlobalSettings::instance()->settings().node("model.initialization"));
220		223
221	QString mode = xml.value("mode", "copy");	224	QString mode = xml.value("mode", "copy");
222	if (mode=="standgrid") {	225	if (mode=="standgrid") {
223	// load a file with saplings per polygon	226	// load a file with saplings per polygon
224	QString filename = xml.value("saplingFile", "");	227	QString filename = xml.value("saplingFile", "");
225	if (filename.isEmpty())	228	if (filename.isEmpty())
226	return;	229	return;
227	filename = GlobalSettings::instance()->path(filename, "init");	230	filename = GlobalSettings::instance()->path(filename, "init");
228	if (!QFile::exists(filename))	231	if (!QFile::exists(filename))
229	throw IException(QString("load-sapling-ini-file: file '%1' does not exist.").arg(filename));	232	throw IException(QString("load-sapling-ini-file: file '%1' does not exist.").arg(filename));
230	CSVFile init_file(filename);	233	CSVFile init_file(filename);
231	int istandid = init_file.columnIndex("stand_id");	234	int istandid = init_file.columnIndex("stand_id");
232	if (istandid==-1)	235	if (istandid==-1)
233	throw IException("Sapling-Init: the init file contains no 'stand_id' column (required in 'standgrid' mode).");	236	throw IException("Sapling-Init: the init file contains no 'stand_id' column (required in 'standgrid' mode).");
234		237
235	int stand_id = -99999;	238	int stand_id = -99999;
236	int ilow = -1, ihigh = 0;	239	int ilow = -1, ihigh = 0;
237	int total = 0;	240	int total = 0;
238	for (int i=0;i<init_file.rowCount();++i) {	241	for (int i=0;i<init_file.rowCount();++i) {
239	int row_stand = init_file.value(i, istandid).toInt();	242	int row_stand = init_file.value(i, istandid).toInt();
240	if (row_stand != stand_id) {	243	if (row_stand != stand_id) {
241	if (stand_id>=0) {	244	if (stand_id>=0) {
242	// process stand	245	// process stand
243	ihigh = i-1; // up to the last	246	ihigh = i-1; // up to the last
244	total += loadSaplingsLIF(stand_id, init_file, ilow, ihigh);	247	total += loadSaplingsLIF(stand_id, init_file, ilow, ihigh);
245	}	248	}
246	ilow = i; // mark beginning of new stand	249	ilow = i; // mark beginning of new stand
247	stand_id = row_stand;	250	stand_id = row_stand;
248	}	251	}
249	}	252	}
250	if (stand_id>=0)	253	if (stand_id>=0)
251	total += loadSaplingsLIF(stand_id, init_file, ilow, init_file.rowCount()-1); // the last stand	254	total += loadSaplingsLIF(stand_id, init_file, ilow, init_file.rowCount()-1); // the last stand
252	qDebug() << "initialization of sapling: total created:" << total;	255	qDebug() << "initialization of sapling: total created:" << total;
253		256
254	}	257	}
255		258
256	}	259	}
257		260
258	void StandLoader::evaluateDebugTrees()	261	void StandLoader::evaluateDebugTrees()
259	{	262	{
260	// evaluate debugging	263	// evaluate debugging
261	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");	264	QString dbg_str = GlobalSettings::instance()->settings().paramValueString("debug_tree");
262	int counter=0;	265	int counter=0;
263	if (!dbg_str.isEmpty()) {	266	if (!dbg_str.isEmpty()) {
264	if (dbg_str == "debugstamp") {	267	if (dbg_str == "debugstamp") {
265	qDebug() << "debug_tree = debugstamp: try touching all trees...";	268	qDebug() << "debug_tree = debugstamp: try touching all trees...";
266	// try to force an error if a stamp is invalid	269	// try to force an error if a stamp is invalid
267	AllTreeIterator at(GlobalSettings::instance()->model());	270	AllTreeIterator at(GlobalSettings::instance()->model());
268	double total_offset=0.;	271	double total_offset=0.;
269	while (Tree *t=at.next()) {	272	while (Tree *t=at.next()) {
270	total_offset += t->stamp()->offset();	273	total_offset += t->stamp()->offset();
271	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))	274	if (!GlobalSettings::instance()->model()->grid()->isIndexValid(t->positionIndex()))
272	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";	275	qDebug() << "evaluateDebugTrees: debugstamp: invalid position found!";
273	}	276	}
274	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";	277	qDebug() << "debug_tree = debugstamp: try touching all trees finished...";
275	return;	278	return;
276	}	279	}
277	TreeWrapper tw;	280	TreeWrapper tw;
278	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables	281	Expression dexp(dbg_str, &tw); // load expression dbg_str and enable external model variables
279	AllTreeIterator at(GlobalSettings::instance()->model());	282	AllTreeIterator at(GlobalSettings::instance()->model());
280	double result;	283	double result;
281	while (Tree *t = at.next()) {	284	while (Tree *t = at.next()) {
282	tw.setTree(t);	285	tw.setTree(t);
283	result = dexp.execute();	286	result = dexp.execute();
284	if (result) {	287	if (result) {
285	t->enableDebugging();	288	t->enableDebugging();
286	counter++;	289	counter++;
287	}	290	}
288	}	291	}
289	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";	292	qDebug() << "evaluateDebugTrees: enabled debugging for" << counter << "trees.";
290	}	293	}
291	}	294	}
292		295
293		296
294	/// load a single init file. Calls loadPicusFile() or loadiLandFile()	297	/// load a single init file. Calls loadPicusFile() or loadiLandFile()
295	/// @param fileName file to load	298	/// @param fileName file to load
296	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"	299	/// @param type init mode. allowed: "picus"/"single" or "iland"/"distribution"
297	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)	300	int StandLoader::loadInitFile(const QString &fileName, const QString &type, int stand_id, ResourceUnit *ru)
298	{	301	{
299	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");	302	QString pathFileName = GlobalSettings::instance()->path(fileName, "init");
300	if (!QFile::exists(pathFileName))	303	if (!QFile::exists(pathFileName))
301	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));	304	throw IException(QString("StandLoader::loadInitFile: File %1 does not exist!").arg(pathFileName));
302		305
303	if (type=="picus" \|\| type=="single") {	306	if (type=="picus" \|\| type=="single") {
304	if (stand_id>0)	307	if (stand_id>0)
305	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);	308	throw IException(QLatin1String("StandLoader::loadInitFile: initialization type %1 currently not supported for stand initilization mode!")+type);
306	return loadPicusFile(pathFileName, ru);	309	return loadPicusFile(pathFileName, ru);
307	}	310	}
308	if (type=="iland" \|\| type=="distribution")	311	if (type=="iland" \|\| type=="distribution")
309	return loadiLandFile(pathFileName, ru, stand_id);	312	return loadiLandFile(pathFileName, ru, stand_id);
310		313
311	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);	314	throw IException(QLatin1String("StandLoader::loadInitFile: unknown initalization.type:")+type);
312	}	315	}
313		316
314	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)	317	int StandLoader::loadPicusFile(const QString &fileName, ResourceUnit *ru)
315	{	318	{
316	QString content = Helper::loadTextFile(fileName);	319	QString content = Helper::loadTextFile(fileName);
317	if (content.isEmpty()) {	320	if (content.isEmpty()) {
318	qDebug() << "file not found: " + fileName;	321	qDebug() << "file not found: " + fileName;
319	return 0;	322	return 0;
320	}	323	}
321	return loadSingleTreeList(content, ru, fileName);	324	return loadSingleTreeList(content, ru, fileName);
322	}	325	}
323		326
324	/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.	327	/** load a list of trees (given by content) to a resource unit. Param fileName is just for error reporting.
325	returns the number of loaded trees.	328	returns the number of loaded trees.
326	*/	329	*/
327	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)	330	int StandLoader::loadSingleTreeList(const QString &content, ResourceUnit *ru, const QString &fileName)
328	{	331	{
329	if (!ru)	332	if (!ru)
330	ru = mModel->ru();	333	ru = mModel->ru();
331	Q_ASSERT(ru!=0);	334	Q_ASSERT(ru!=0);
332		335
333	QPointF offset = ru->boundingBox().topLeft();	336	QPointF offset = ru->boundingBox().topLeft();
334	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	337	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
335		338
336	QString my_content(content);	339	QString my_content(content);
337	// cut out the <trees> </trees> part if present	340	// cut out the <trees> </trees> part if present
338	if (content.contains("<trees>")) {	341	if (content.contains("<trees>")) {
339	QRegExp rx(".<trees>(.)</trees>.*");	342	QRegExp rx(".<trees>(.)</trees>.*");
340	rx.indexIn(content, 0);	343	rx.indexIn(content, 0);
341	if (rx.capturedTexts().count()<1)	344	if (rx.capturedTexts().count()<1)
342	return 0;	345	return 0;
343	my_content = rx.cap(1).trimmed();	346	my_content = rx.cap(1).trimmed();
344	}	347	}
345		348
346	CSVFile infile;	349	CSVFile infile;
347	infile.loadFromString(my_content);	350	infile.loadFromString(my_content);
348		351
349		352
350	int iID = infile.columnIndex("id");	353	int iID = infile.columnIndex("id");
351	int iX = infile.columnIndex("x");	354	int iX = infile.columnIndex("x");
352	int iY = infile.columnIndex("y");	355	int iY = infile.columnIndex("y");
353	int iBhd = infile.columnIndex("bhdfrom");	356	int iBhd = infile.columnIndex("bhdfrom");
354	if (iBhd<0)	357	if (iBhd<0)
355	iBhd = infile.columnIndex("dbh");	358	iBhd = infile.columnIndex("dbh");
356	double height_conversion = 100.;	359	double height_conversion = 100.;
357	int iHeight =infile.columnIndex("treeheight");	360	int iHeight =infile.columnIndex("treeheight");
358	if (iHeight<0) {	361	if (iHeight<0) {
359	iHeight = infile.columnIndex("height");	362	iHeight = infile.columnIndex("height");
360	height_conversion = 1.; // in meter	363	height_conversion = 1.; // in meter
361	}	364	}
362	int iSpecies = infile.columnIndex("species");	365	int iSpecies = infile.columnIndex("species");
363	int iAge = infile.columnIndex("age");	366	int iAge = infile.columnIndex("age");
364	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)	367	if (iX==-1 \|\| iY==-1 \|\| iBhd==-1 \|\| iSpecies==-1 \|\| iHeight==-1)
365	throw IException(QString("Initfile %1 is not valid!\nRequired columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));	368	throw IException(QString("Initfile %1 is not valid!\nRequired columns are: x,y, bhdfrom or dbh, species, treeheight or height.").arg(fileName));
366		369
367	double dbh;	370	double dbh;
368	bool ok;	371	bool ok;
369	int cnt=0;	372	int cnt=0;
370	QString speciesid;	373	QString speciesid;
371	for (int i=1;i<infile.rowCount();i++) {	374	for (int i=1;i<infile.rowCount();i++) {
372	dbh = infile.value(i, iBhd).toDouble();	375	dbh = infile.value(i, iBhd).toDouble();
373		376
374	//if (dbh<5.)	377	//if (dbh<5.)
375	// continue;	378	// continue;
376		379
377	QPointF f;	380	QPointF f;
378	if (iX>=0 && iY>=0) {	381	if (iX>=0 && iY>=0) {
379	f.setX( infile.value(i, iX).toDouble() );	382	f.setX( infile.value(i, iX).toDouble() );
380	f.setY( infile.value(i, iY).toDouble() );	383	f.setY( infile.value(i, iY).toDouble() );
381	f+=offset;	384	f+=offset;
382		385
383	}	386	}
384	// position valid?	387	// position valid?
385	if (!mModel->heightGrid()->valueAt(f).isValid())	388	if (!mModel->heightGrid()->valueAt(f).isValid())
386	continue;	389	continue;
387	Tree &tree = ru->newTree();	390	Tree &tree = ru->newTree();
388	tree.setPosition(f);	391	tree.setPosition(f);
389	if (iID>=0)	392	if (iID>=0)
390	tree.setId(infile.value(i, iID).toInt());	393	tree.setId(infile.value(i, iID).toInt());
391		394
392	tree.setDbh(dbh);	395	tree.setDbh(dbh);
393	tree.setHeight(infile.value(i,iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary	396	tree.setHeight(infile.value(i,iHeight).toDouble()/height_conversion); // convert from Picus-cm to m if necessary
394		397
395	speciesid = infile.value(i,iSpecies).toString();	398	speciesid = infile.value(i,iSpecies).toString();
396	int picusid = speciesid.toInt(&ok);	399	int picusid = speciesid.toInt(&ok);
397	if (ok) {	400	if (ok) {
398	int idx = picusSpeciesIds.indexOf(picusid);	401	int idx = picusSpeciesIds.indexOf(picusid);
399	if (idx==-1)	402	if (idx==-1)
400	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));	403	throw IException(QString("Loading init-file: invalid Picus-species-id. Species: %1").arg(picusid));
401	speciesid = iLandSpeciesIds[idx];	404	speciesid = iLandSpeciesIds[idx];
402	}	405	}
403	Species *s = speciesSet->species(speciesid);	406	Species *s = speciesSet->species(speciesid);
404	if (!ru \|\| !s)	407	if (!ru \|\| !s)
405	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));	408	throw IException(QString("Loading init-file: either resource unit or species invalid. Species: %1").arg(speciesid));
406	tree.setSpecies(s);	409	tree.setSpecies(s);
407		410
408	ok = true;	411	ok = true;
409	if (iAge>=0)	412	if (iAge>=0)
410	tree.setAge(infile.value(i, iAge).toInt(&ok), tree.height()); // this is a real age	413	tree.setAge(infile.value(i, iAge).toInt(&ok), tree.height()); // this is a real age
411	if (iAge<0 \|\| !ok \|\| tree.age()==0)	414	if (iAge<0 \|\| !ok \|\| tree.age()==0)
412	tree.setAge(0, tree.height()); // no real tree age available	415	tree.setAge(0, tree.height()); // no real tree age available
413		416
414	tree.setRU(ru);	417	tree.setRU(ru);
415	tree.setup();	418	tree.setup();
416	cnt++;	419	cnt++;
417	}	420	}
418	return cnt;	421	return cnt;
419	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";	422	//qDebug() << "loaded init-file contained" << lines.count() <<"lines.";
420	//qDebug() << "lines: " << lines;	423	//qDebug() << "lines: " << lines;
421	}	424	}
422		425
423	/** initialize trees on a resource unit based on dbh distributions.	426	/** initialize trees on a resource unit based on dbh distributions.
424	use a fairly clever algorithm to determine tree positions.	427	use a fairly clever algorithm to determine tree positions.
425	see http://iland.boku.ac.at/initialize+trees	428	see http://iland.boku.ac.at/initialize+trees
426	@param content tree init file (including headers) in a string	429	@param content tree init file (including headers) in a string
427	@param ru resource unit	430	@param ru resource unit
428	@param fileName source file name (for error reporting)	431	@param fileName source file name (for error reporting)
429	@return number of trees added	432	@return number of trees added
430	*/	433	*/
431	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)	434	int StandLoader::loadDistributionList(const QString &content, ResourceUnit *ru, int stand_id, const QString &fileName)
432	{	435	{
433	int total_count = parseInitFile(content, fileName, ru);	436	int total_count = parseInitFile(content, fileName, ru);
434	if (total_count==0)	437	if (total_count==0)
435	return 0;	438	return 0;
436		439
437		440
438	// setup the random distribution	441	// setup the random distribution
439	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");	442	QString density_func = GlobalSettings::instance()->settings().value("model.initialization.randomFunction", "1-x^2");
440	if (logLevelInfo()) qDebug() << "density function:" << density_func;	443	if (logLevelInfo()) qDebug() << "density function:" << density_func;
441	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {	444	if (!mRandom \|\| (mRandom->densityFunction()!= density_func)) {
442	if (mRandom)	445	if (mRandom)
443	delete mRandom;	446	delete mRandom;
444	mRandom=new RandomCustomPDF(density_func);	447	mRandom=new RandomCustomPDF(density_func);
445	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;	448	if (logLevelInfo()) qDebug() << "new probabilty density function:" << density_func;
446	}	449	}
447	if (stand_id>0) {	450	if (stand_id>0) {
448	// execute stand based initialization	451	// execute stand based initialization
449	executeiLandInitStand(stand_id);	452	executeiLandInitStand(stand_id);
450	} else {	453	} else {
451	// exeucte the initialization based on single resource units	454	// exeucte the initialization based on single resource units
452	executeiLandInit(ru);	455	executeiLandInit(ru);
453	ru->cleanTreeList();	456	ru->cleanTreeList();
454	}	457	}
455	return total_count;	458	return total_count;
456		459
457	}	460	}
458		461
459	int StandLoader::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru)	462	int StandLoader::parseInitFile(const QString &content, const QString &fileName, ResourceUnit* ru)
460	{	463	{
461	if (!ru)	464	if (!ru)
462	ru = mModel->ru();	465	ru = mModel->ru();
463	Q_ASSERT(ru!=0);	466	Q_ASSERT(ru!=0);
464	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU	467	SpeciesSet *speciesSet = ru->speciesSet(); // of default RU
465	Q_ASSERT(speciesSet!=0);	468	Q_ASSERT(speciesSet!=0);
466		469
467	//DebugTimer t("StandLoader::loadiLandFile");	470	//DebugTimer t("StandLoader::loadiLandFile");
468	CSVFile infile;	471	CSVFile infile;
469	infile.loadFromString(content);	472	infile.loadFromString(content);
470		473
471	int icount = infile.columnIndex("count");	474	int icount = infile.columnIndex("count");
472	int ispecies = infile.columnIndex("species");	475	int ispecies = infile.columnIndex("species");
473	int idbh_from = infile.columnIndex("dbh_from");	476	int idbh_from = infile.columnIndex("dbh_from");
474	int idbh_to = infile.columnIndex("dbh_to");	477	int idbh_to = infile.columnIndex("dbh_to");
475	int ihd = infile.columnIndex("hd");	478	int ihd = infile.columnIndex("hd");
476	int iage = infile.columnIndex("age");	479	int iage = infile.columnIndex("age");
477	int idensity = infile.columnIndex("density");	480	int idensity = infile.columnIndex("density");
478	if (icount<0 \|\| ispecies<0 \|\| idbh_from<0 \|\| idbh_to<0 \|\| ihd<0 \|\| iage<0)	481	if (icount<0 \|\| ispecies<0 \|\| idbh_from<0 \|\| idbh_to<0 \|\| ihd<0 \|\| iage<0)
479	throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));	482	throw IException(QString("load-ini-file: file '%1' containts not all required fields (count, species, dbh_from, dbh_to, hd, age).").arg(fileName));
480		483
481	int istandid = infile.columnIndex("stand_id");	484	int istandid = infile.columnIndex("stand_id");
482	mInitItems.clear();	485	mInitItems.clear();
483	mStandInitItems.clear();	486	mStandInitItems.clear();
484		487
485	InitFileItem item;	488	InitFileItem item;
486	bool ok;	489	bool ok;
487	int total_count = 0;	490	int total_count = 0;
488	for (int row=0;row<infile.rowCount();row++) {	491	for (int row=0;row<infile.rowCount();row++) {
489	item.count = infile.value(row, icount).toDouble();	492	item.count = infile.value(row, icount).toDouble();
490	total_count += item.count;	493	total_count += item.count;
491	item.dbh_from = infile.value(row, idbh_from).toDouble();	494	item.dbh_from = infile.value(row, idbh_from).toDouble();
492	item.dbh_to = infile.value(row, idbh_to).toDouble();	495	item.dbh_to = infile.value(row, idbh_to).toDouble();
493	item.hd = infile.value(row, ihd).toDouble();	496	item.hd = infile.value(row, ihd).toDouble();
494	if (item.hd==0. \|\| item.dbh_from / 100. * item.hd < 4.)	497	if (item.hd==0. \|\| item.dbh_from / 100. * item.hd < 4.)
495	qWarning() << QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ;	498	qWarning() << QString("load init file: file '%1' tries to init trees below 4m height. hd=%2, dbh=%3.").arg(fileName).arg(item.hd).arg(item.dbh_from) ;

Subversion Repositories public iLand

(root)/src/core/standloader.cpp @ 1222 - Rev 1190 → 1203