WebSVN - public iLand - Diff - Rev 671 and 697


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"

/** @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;
}


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()+(*p)->boundingBox().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", "");

    Tree::resetStatistics();

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

    // copy trees from first unit to all other units:
    if (copy_mode=="copy") {
        loadInitFile(fileName, type);
        copyTrees();
        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);
            ru->setRandomGenerator(); //
            // 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
    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;
                loadInitFile(file_name, type, key, NULL);
            }
        }
        evaluateDebugTrees();
        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()) {
       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.";
}

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;
    QPoint pixelOffset;
    ResourceUnit *resource_unit;
    SInitPixel(): basal_area(0.), resource_unit(0) {}
};

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(ru->randomGenerator(), 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(ru->randomGenerator())*(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();

    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

    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(grid->grid().indexOf(i)));
       pixel_list.append(p);
    }
    double area_factor = grid->area(stand_id) / 10000.;

    int key;
    double rand_val, rand_fraction;
    int total_count = 0;
    foreach(const InitFileItem &item, mInitItems) {
        rand_fraction = fabs(double(item.density));
        int count = item.count * area_factor + 0.5; // round
        for (int i=0;i<count;i++) {


            // 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
            key = limit(int(pixel_list.count()*rand_val), 0, pixel_list.count()-1); // get from random number generator

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

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

Subversion Repositories public iLand

(root)/src/core/standloader.cpp @ 1222 - Rev 671 → 697