Subversion Repositories public iLand

/********************************************************************************************

**    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 "dynamicstandout.h"

#include "debugtimer.h"

#include "statdata.h"

#include "model.h"

#include "resourceunit.h"

#include "species.h"

#include "expressionwrapper.h"

DynamicStandOut::DynamicStandOut()

{

    setName("dynamic stand output by species/RU", "dynamicstand");

    setDescription("Userdefined outputs for tree aggregates for each stand or species.\n"\

                   "Technically, each field is calculated 'live', i.e. it is looped over all trees, and eventually the statistics (percentiles) "\

                   "are calculated. The aggregated values are not scaled to any area unit.\n" \

                   "!!!Specifying the aggregation\n" \

                   "The ''by_species'' and ''by_ru'' option allow to define the aggregation level. When ''by_species'' is set to ''true'', " \

                   "a row for each species will be created, otherwise all trees of all species are aggregated to one row. " \

                   "Similarly, ''by_ru''=''true'' means outputs for each resource unit, while a value of ''false'' aggregates over the full project area.\n" \

                   "!!!Specifying filters\n" \

                   "You can use the 'rufilter' and 'treefilter' XML settings to reduce the limit the output to a subset of resource units / trees. " \

                   "Both filters are valid expressions (for resource unit level and tree level, respectively). For example, a ''treefilter'' of 'speciesindex=0' reduces the output to just one species.\n" \

                   "The ''condition'' filter is (when present) evaluated and the output is only executed when ''condition'' is true (variable='year') This can be used to constrain the output to specific years (e.g. 'in(year,100,200,300)' produces output only for the given year.\n" \

                   "!!!Specifying data columns\n"

                   "Each field is defined as: ''field.aggregatio''n (separated by a dot). A ''field'' is a valid [Expression]. ''Aggregation'' is one of the following:  " \

                   "mean, sum, min, max, p25, p50, p75, p5, 10, p90, p95 (pXX=XXth percentile), sd (std.dev.).\n" \

                   "Complex expression are allowed, e.g: if(dbh>50,1,0).sum (-> counts trees with dbh>50)");

    columns() << OutputColumn::year() << OutputColumn::ru()  << OutputColumn::id() << OutputColumn::species();

    // other colums are added during setup...

}

const QStringList aggList = QStringList() << "mean" << "sum" << "min" << "max" << "p25" << "p50" << "p75" << "p5"<< "p10" << "p90" << "p95" << "sd";

void DynamicStandOut::setup()

{

    QString filter = settings().value(".rufilter","");

    QString tree_filter = settings().value(".treefilter","");

    QString fieldList = settings().value(".columns", "");

    QString condition = settings().value(".condition", "");

    if (fieldList.isEmpty())

        return;

    mRUFilter.setExpression(filter);

    mTreeFilter.setExpression(tree_filter);

    mCondition.setExpression(condition);

    // clear columns

    columns().erase(columns().begin()+4, columns().end());

    mFieldList.clear();

    // setup fields

   if (!fieldList.isEmpty()) {

       QRegExp rx("([^\\.]+).(\\w+)[,\\s]*"); // two parts: before dot and after dot, and , + whitespace at the end

        int pos=0;

        QString field, aggregation;

        TreeWrapper tw;

        while ((pos = rx.indexIn(fieldList, pos)) != -1) {

            pos += rx.matchedLength();

            field = rx.cap(1); // field / expresssion

            aggregation = rx.cap(2);

            mFieldList.append(SDynamicField());

            // parse field

            if (field.count()>0 && !field.contains('(')) {

                // simple expression

                mFieldList.back().var_index = tw.variableIndex(field);

            } else {

                // complex expression

                mFieldList.back().var_index=-1;

                mFieldList.back().expression = field;

            }

            mFieldList.back().agg_index = aggList.indexOf(aggregation);

            if (mFieldList.back().agg_index==-1)

                 throw IException(QString("Invalid aggregate expression for dynamic output: %1\nallowed:%2")

                                          .arg(aggregation).arg(aggList.join(" ")));

             QString stripped_field=QString("%1_%2").arg(field, aggregation);

             stripped_field.replace(QRegExp("[\\[\\]\\,\\(\\)<>=!\\s]"), "_");

             stripped_field.replace("__", "_");

             columns() << OutputColumn(stripped_field, field, OutDouble);

        }

   }

}

void DynamicStandOut::exec()

{

    if (mFieldList.count()==0)

        return;

    if (!mCondition.isEmpty())

        if (!mCondition.calculate(GlobalSettings::instance()->currentYear()))

            return;

    DebugTimer t("dynamic stand output");

    bool per_species = GlobalSettings::instance()->settings().valueBool("output.dynamicstand.by_species", true);

    bool per_ru = GlobalSettings::instance()->settings().valueBool("output.dynamicstand.by_ru", true);

    if (per_ru) {

        // when looping over resource units, do it differently (old way)

        extractByResourceUnit(per_species);

        return;

    }

    Model *m = GlobalSettings::instance()->model();

    QVector<double> data; //statistics data

    TreeWrapper tw;

    Expression custom_expr;

    StatData stat; // statistcs helper class

    // grouping

    QVector<const Tree*> trees;

    for (QList<Species*>::const_iterator species = m->speciesSet()->activeSpecies().constBegin();species!=m->speciesSet()->activeSpecies().constEnd();++species) {

        trees.clear();

        AllTreeIterator all_trees(m);

        while (const Tree *t = all_trees.nextLiving()) {

            if (per_species && t->species() != *species)

                continue;

            trees.push_back(t);

        }

        if (trees.size()==0)

            continue;

        // dynamic calculations

        foreach (const SDynamicField &field, mFieldList) {

            if (!field.expression.isEmpty()) {

                // setup dynamic dynamic expression if present

                custom_expr.setExpression(field.expression);

                custom_expr.setModelObject(&tw);

            }

            // fetch data values from the trees

            data.clear();

            foreach(const Tree *t, trees) {

                tw.setTree(t);

                if (field.var_index>=0)

                    data.push_back(tw.value(field.var_index));

                else

                    data.push_back(custom_expr.execute());

            }

            // constant values (if not already present)

            if (isRowEmpty()) {

                *this << currentYear() << -1 << -1;

                if (per_species)

                    *this << (*species)->id();

                else

                    *this << "";

            }

            // calculate statistics

            stat.setData(data);

            // aggregate

            double value;

            switch (field.agg_index) {

            case 0: value = stat.mean(); break;

            case 1: value = stat.sum(); break;

            case 2: value = stat.min(); break;

            case 3: value = stat.max(); break;

            case 4: value = stat.percentile25(); break;

            case 5: value = stat.median(); break;

            case 6: value = stat.percentile75(); break;

            case 7: value = stat.percentile(5); break;

            case 8: value = stat.percentile(10); break;

            case 9: value = stat.percentile(90); break;

            case 10: value = stat.percentile(95); break;

            case 11: value = stat.standardDev(); break;

            default: value = 0.; break;

            }

            // add current value to output

            *this << value;

        }

        if (!isRowEmpty())

            writeRow();

        if (!per_species)

            break;

    }

}

void DynamicStandOut::extractByResourceUnit(const bool by_species)

{

    if (mFieldList.count()==0)

        return;

    Model *m = GlobalSettings::instance()->model();

    QVector<double> data; //statistics data

    StatData stat; // statistcs helper class

    TreeWrapper tw;

    RUWrapper ruwrapper;

    mRUFilter.setModelObject(&ruwrapper);

    Expression custom_expr;

    foreach(ResourceUnit *ru, m->ruList()) {

        if (ru->id()==-1)

            continue; // do not include if out of project area

        // test filter

        if (!mRUFilter.isEmpty()) {

            ruwrapper.setResourceUnit(ru);

            if (!mRUFilter.execute())

                continue;

        }

        foreach(const ResourceUnitSpecies *rus, ru->ruSpecies()) {

            if (by_species && rus->constStatistics().count()==0)

                continue;

            // dynamic calculations

            foreach (const SDynamicField &field, mFieldList) {

                if (!field.expression.isEmpty()) {

                    // setup dynamic dynamic expression if present

                    custom_expr.setExpression(field.expression);

                    custom_expr.setModelObject(&tw);

                }

                data.clear();

                bool has_trees = false;

                foreach(const Tree &tree, ru->trees()) {

                    if (by_species && tree.species()->index()!=rus->species()->index())

                        continue;

                    if (tree.isDead())

                        continue;

                    tw.setTree(&tree);

                    // apply treefilter

                    if (!mTreeFilter.isEmpty()) {

                        mTreeFilter.setModelObject(&tw);

                        if (!mTreeFilter.execute())

                            continue;

                    }

                    has_trees = true;

                    if (field.var_index>=0)

                        data.push_back(tw.value(field.var_index));

                    else

                        data.push_back(custom_expr.execute());

                }

                // do nothing if no trees are avaiable

                if (!has_trees)

                    continue;

                if (isRowEmpty()) {

                    *this << currentYear()  << ru->index() << ru->id();

                    if (by_species)

                        *this << rus->species()->id();

                    else

                        *this << "";

                 }

                // calculate statistics

                stat.setData(data);

                // aggregate

                double value;

                switch (field.agg_index) {

                case 0: value = stat.mean(); break;

                case 1: value = stat.sum(); break;

                case 2: value = stat.min(); break;

                case 3: value = stat.max(); break;

                case 4: value = stat.percentile25(); break;

                case 5: value = stat.median(); break;

                case 6: value = stat.percentile75(); break;

                case 7: value = stat.percentile(5); break;

                case 8: value = stat.percentile(10); break;

                case 9: value = stat.percentile(90); break;

                case 10: value = stat.percentile(95); break;

                case 11: value = stat.standardDev(); break;

                default: value = 0.; break;

                }

                // add current value to output

                *this << value;

            } // foreach (field)

            if (!isRowEmpty())

                writeRow();

            if (!by_species)

                break;

        } //foreach species

    } // foreach ressource unit

}