WebSVN - public iLand - Diff - Rev 697 and 699


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

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


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

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

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

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


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

 
 
 
 
 
 
 

    // 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;
    }
    if (copy_mode=="snapshot") {
        // load a snapshot from a file
        Snapshot shot;

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

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

Subversion Repositories public iLand

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