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 "abe_global.h"
#include "actsalvage.h"

#include "fmstp.h"
#include "fmstand.h"
#include "fomescript.h"
#include "scheduler.h"
#include "forestmanagementengine.h"
#include "fmtreelist.h"

#include "tree.h"
#include "expression.h"
#include "expressionwrapper.h"
#include "mapgrid.h"
#include "helper.h"

namespace ABE {

/** @class ActSalvage
    @ingroup abe
    The ActSalvage class handles salvage logging after disturbances.

  */

ActSalvage::ActSalvage(FMSTP *parent): Activity(parent)
{
    mCondition = 0;
    mMaxPreponeActivity = 0;

    mBaseActivity.setIsSalvage(true);
    mBaseActivity.setIsRepeating(true);
    mBaseActivity.setExecuteImmediate(true);

}

ActSalvage::~ActSalvage()
{
    if (mCondition)
        delete mCondition;
}

void ActSalvage::setup(QJSValue value)
{
    Activity::setup(value); // setup base events
    events().setup(value, QStringList() << "onBarkBeetleAttack");

    QString condition = FMSTP::valueFromJs(value, "disturbanceCondition").toString();
    if (!condition.isEmpty() && condition!="undefined") {
        mCondition = new Expression(condition);
    }
    mMaxPreponeActivity = FMSTP::valueFromJs(value, "maxPrepone", "0").toInt();
    mThresholdSplit = FMSTP::valueFromJs(value, "thresholdSplitStand", "0.1").toNumber();
    mThresholdClear = FMSTP::valueFromJs(value, "thresholdClearStand", "0.9").toNumber();
    mThresholdMinimal = FMSTP::valueFromJs(value, "thresholdIgnoreDamage", "5").toNumber();
    mDebugSplit = FMSTP::boolValueFromJs(value, "debugSplit", false);

}

bool ActSalvage::execute(FMStand *stand)
{


    if (stand->property("_run_salvage").toBool()) {
        // 2nd phase: do the after disturbance cleanup of a stand.
        bool simu = stand->currentFlags().isDoSimulate();
        stand->currentFlags().setDoSimulate(false);
        stand->currentFlags().setFinalHarvest(true); // this should be accounted as "final harvest"
        // execute the "onExecute" event
        bool result =  Activity::execute(stand);
        stand->currentFlags().setDoSimulate(simu);
        stand->currentFlags().setFinalHarvest(false);
        stand->setProperty("_run_salvage", false);
        return result;
    }

    // the salvaged timber is already accounted for - so nothing needs to be done here.
    // however, we check if there is a planned activity for the stand which could be executed sooner
    // than planned.
    bool preponed = const_cast<FMUnit*>(stand->unit())->scheduler()->forceHarvest(stand, mMaxPreponeActivity);
    if (stand->trace())
        qCDebug(abe) << "Salvage activity executed. Changed scheduled activites (preponed): " << preponed;

    const_cast<FMUnit*>(stand->unit())->scheduler()->addExtraHarvest(stand, stand->totalHarvest(), Scheduler::Salvage);
    // check if we should re-assess the stand grid (after large disturbances)
    // as a preliminary check we only look closer, if we have more than  x m3/ha of damage.
    if (stand->disturbedTimber()/stand->area() > mThresholdMinimal)
        checkStandAfterDisturbance(stand);

    // the harvest happen(ed) anyways.
    //stand->resetHarvestCounter(); // set back to zero...
    return true;
}

QStringList ActSalvage::info()
{
    QStringList lines = Activity::info();
    lines << QString("condition: %1").arg(mCondition?mCondition->expression():"-");
    lines << QString("maxPrepone: %1").arg(mMaxPreponeActivity);
    return lines;
}

bool ActSalvage::evaluateRemove(Tree *tree) const
{
    if (!mCondition)
        return true; // default: remove all trees
    TreeWrapper tw(tree);
    bool result = static_cast<bool>( mCondition->execute(0, &tw) );
    return result;
}

bool ActSalvage::barkbeetleAttack(FMStand *stand, double generations, int infested_px_ha)
{

    //QJSValue params;
    QJSValueList params=QJSValueList() << QJSValue(generations) << QJSValue(infested_px_ha);

    QJSValue result = events().run(QStringLiteral("onBarkBeetleAttack"), stand, &params);
    if (!result.isBool())
        qCDebug(abe) << "Salvage-Activity:onBarkBeetleAttack: expecting a boolean return";
    return result.toBool();
}

void ActSalvage::checkStandAfterDisturbance(FMStand *stand)
{
    //
    FMTreeList *trees = ForestManagementEngine::instance()->scriptBridge()->treesObj();
    //trees->runGrid();
    trees->prepareStandGrid(QStringLiteral("height"), QString());

    const int min_split_size = 50; // min size (100=1ha)
    FloatGrid &grid = trees->standGrid();
    static int no_split = 0;
    if (mDebugSplit)
        trees->exportStandGrid(QString("temp/height_%1.txt").arg(++no_split));

    float h_max = grid.max();

    double r_low;
    int h_lower = 0, h_higher=0;
    if (h_max==0.f) {
        // total disturbance...
        r_low = 1.;
    } else {
        // check coverage of disturbed area.
        for (float *p=grid.begin(); p!=grid.end(); ++p)
            if (*p>=0.f) {
                if (*p < h_max*0.33)
                    ++h_lower;
                else
                    ++h_higher;
            }
        if (h_lower==0 && h_higher==0)
            return;
        r_low = h_lower / double(h_lower+h_higher);
    }

    if (r_low < mThresholdSplit || (r_low<0.5 && h_lower<min_split_size)) {
        // no big damage: return and do nothing
        return;
    }

    // restart if a large fraction is cleared, or if the remaining forest is <0.25ha
    if (r_low > mThresholdClear || (r_low>0.5 && h_higher<min_split_size)) {
        // total disturbance: restart rotation...
        qCDebug(abe) << "ActSalvage: total damage for stand" << stand->id() << "Restarting rotation.";
        stand->setProperty("_run_salvage", true);
        stand->reset(stand->stp());
        return;
    }
    // medium disturbance: check if need to split the stand area:
    Grid<int> my_map(grid.cellsize(), grid.sizeX(), grid.sizeY());
    GridRunner<float> runner(&grid);
    GridRunner<int> id_runner(&my_map);
    float *neighbors[8];
    int n_empty=0;
    while (runner.next() && id_runner.next()) {
        if (*runner.current()==-1.) {
            *id_runner.current() = -1;
            continue;
        }
        runner.neighbors8(neighbors);
        double empty = 0.;
        int valid = 0;
        for (int i=0;i<8;++i) {
            if (neighbors[i] && *neighbors[i]<h_max*0.33)
                empty++;
            if (neighbors[i])
                valid++;
        }
        if (valid)
            empty /= double(valid);
        // empty cells are marked with 0; areas covered by forest set to stand_id; -1: out-of-stand areas
        // if a cell is empty, some neighbors (i.e. >50%) need to be empty too;
        // if a cell is *not* empty, it has to be surrounded by a larger fraction of empty points (75%)
        if ( (*runner.current()<h_max*0.33 && empty>0.5)
             || (empty>=0.75) ) {
            *id_runner.current() = 0;
            n_empty++;
        } else {
            *id_runner.current() = stand->id();
        }
    }
    if (mDebugSplit)
        Helper::saveToTextFile(GlobalSettings::instance()->path(QString("temp/split_before_%1.txt").arg(no_split)), gridToESRIRaster(my_map) );


    // now flood-fill 0ed areas....
    // if the "new" areas are too small (<0.25ha), then nothing happens.
    QVector<QPair<int,int> > cleared_small_areas; // areas of cleared "patches"
    QVector<QPair<int,int> > stand_areas; // areas of remaining forested "patches"
    int fill_color = -1;
    int stand_fill_color = stand->id() + 1000;
    id_runner.reset();
    while (id_runner.next()) {
        if (*id_runner.current()==0) {
            int s = floodFillHelper(my_map, id_runner.currentIndex(), 0, --fill_color);
            cleared_small_areas.push_back(QPair<int,int>(fill_color, s)); // patch size
        } else if (*id_runner.current()==stand->id()) {
            int s=floodFillHelper(my_map, id_runner.currentIndex(), stand->id(), stand_fill_color);
            stand_areas.push_back(QPair<int,int>(stand_fill_color,s));
            stand_fill_color++;
        }
    }
    if (mDebugSplit)
        Helper::saveToTextFile(GlobalSettings::instance()->path(QString("temp/split_stands_%1.txt").arg(no_split)), gridToESRIRaster(my_map) );


    // special case: remainnig forest are only small patches
    int max_size=0;
    for (int i=0;i<stand_areas.size();++i)
        max_size=std::max(max_size, stand_areas[i].second);
    if (max_size<min_split_size) {
        // total disturbance: restart rotation...
        qCDebug(abe) << "ActSalvage: total damage for stand" << stand->id() << "(remaining patches too small). Restarting rotation.";
        stand->setProperty("_run_salvage", true);
        stand->reset(stand->stp());
        return;
    }

    // clear small areas
    QVector<int> neighbor_ids;
    bool finished=false;
    int iter=100;
    while (!finished && cleared_small_areas.size()>0 && --iter>0) {

        // find smallest area....
        int i_min=-1;
        for (int i=0;i<cleared_small_areas.size();++i) {
            if (cleared_small_areas[i].second<min_split_size) {
                if (i_min==-1 || (i_min>-1 && cleared_small_areas[i].second<cleared_small_areas[i_min].second))
                    i_min = i;
            }
        }
        if (i_min==-1) {
            finished = true;
            continue;
        }

        // loook for neighbors of the area
        // attach to largest "cleared" neighbor (if such a neighbor exists)
        neighborFinderHelper(my_map,neighbor_ids, cleared_small_areas[i_min].first);
        if (neighbor_ids.size()==0) {
            // patch fully surrounded by "out of project area". We'll add it to the *first* stand map entry
            neighbor_ids.append(stand_areas.first().first);
        }
        // look for "empty patches" first
        int i_empty=-1; int max_size=0;
        for (int i=0;i<cleared_small_areas.size();++i) {
            if (neighbor_ids.contains(cleared_small_areas[i].first))
                if (cleared_small_areas[i].second>max_size) {
                    i_empty=i;
                    max_size=cleared_small_areas[i].second;
                }
        }
        if (i_empty>-1) {
            // replace "i_min" with "i_empty"
            int r = replaceValueHelper(my_map, cleared_small_areas[i_min].first, cleared_small_areas[i_empty].first);
            cleared_small_areas[i_empty].second += r;
            cleared_small_areas.remove(i_min);
            continue;
        }


        if (stand_areas.size()>0) {
            // attach to largest stand part which is a neighbor

            int i_empty=-1; int max_size=0;
            for (int i=0;i<stand_areas.size();++i) {
                if (neighbor_ids.contains(stand_areas[i].first))
                    if (stand_areas[i].second>max_size) {
                        i_empty=i;
                        max_size=stand_areas[i].second;
                    }
            }
            if (i_empty>-1) {
                // replace "i_min" with "i_empty"
                int r = replaceValueHelper(my_map, cleared_small_areas[i_min].first, stand_areas[i_empty].first);
                stand_areas[i_empty].second += r;
                cleared_small_areas.remove(i_min);
            }
        }
        if (iter==3)
            qCDebug(abe) << "ActSalvage:Loop1: no solution.";
    }


    // clear small forested stands
    finished = false;
    iter = 100;
    while (!finished && --iter>0) {
        finished=true;
        for (int i=0;i<stand_areas.size();++i) {
            if (stand_areas[i].second < min_split_size) {
                neighborFinderHelper(my_map, neighbor_ids, stand_areas[i].first);

                if (neighbor_ids.size()>0) {
                    int r = replaceValueHelper(my_map, stand_areas[i].first, neighbor_ids[0]);
                    if (neighbor_ids[0]>0) {
                        // another stand
                        for (int j=0;j<stand_areas.size();++j)
                            if (stand_areas[j].first == neighbor_ids[0]) {
                                stand_areas[j].second += r;
                            }
                    } else {
                        // clearing
                        for (int j=0;j<cleared_small_areas.size();++j)
                            if (cleared_small_areas[j].first == neighbor_ids[0]) {
                                cleared_small_areas[j].second += r;
                            }
                    }

                    stand_areas.remove(i);
                    finished=false;
                    break;
                }
            }
        }
    }
    if (iter==0)
        qCDebug(abe) << "ActSalvage:Loop2: no solution.";
    if (mDebugSplit)
        Helper::saveToTextFile(GlobalSettings::instance()->path(QString("temp/split_final_%1.txt").arg(no_split)), gridToESRIRaster(my_map) );


    // determine final new stands....
    QVector<int> new_stands; // internal ids that should become new stands

    for (int i=0;i<cleared_small_areas.size();++i)
        new_stands.push_back(cleared_small_areas[i].first);

    // only add new stands - keep the old stand as is
    //    if (new_stands.size()>0) {
//        // if there are no "cleared" parts, we keep the stand as is.
//        for (int i=0;i<stand_areas.size();++i)
//            if (stand_areas[i].first != stand->id()+1000)
//                new_stands.push_back(stand_areas[i].first);
//    }


    for (int i=0;i<new_stands.size(); ++i) {
        // ok: we have new stands. Now do the actual splitting
        FMStand *new_stand = ForestManagementEngine::instance()->splitExistingStand(stand);
        // copy back to the stand grid
        GridRunner<int> sgrid(ForestManagementEngine::instance()->standGrid()->grid(), grid.metricRect());
        id_runner.reset();
        int n_px=0;
        while (sgrid.next() && id_runner.next()) {
            if (*id_runner.current() == new_stands[i]) {
                *sgrid.current() = new_stand->id();
                ++n_px;
            }
        }

        // the new stand  is prepared.
        // at the end of this years execution, the stand will be re-evaluated.
        new_stand->setInitialId(stand->id());
        // year of splitting: all the area of the stand is still accounted for for the "old" stand
        // in the next year (after the update of the stand grid), the old stand shrinks and the new
        // stands get their correct size.
        // new_stand->setArea(n_px / (cHeightSize*cHeightSize));
        new_stand->setProperty("_run_salvage", true);
        new_stand->reset(stand->stp());
        qCDebug(abe) << "ActSalvage: new stand" << new_stand->id() << "parent stand" << stand->id() << "#split:" << no_split;
    }
}

// quick and dirty implementation of the flood fill algroithm.
// based on: http://en.wikipedia.org/wiki/Flood_fill
// returns the number of pixels colored
int ActSalvage::floodFillHelper(Grid<int> &grid, QPoint start, int old_color, int color)
{
    QQueue<QPoint> pqueue;
    pqueue.enqueue(start);
    int found = 0;
    while (!pqueue.isEmpty()) {
        QPoint p = pqueue.dequeue();
        if (!grid.isIndexValid(p))
            continue;
        if (grid.valueAtIndex(p)==old_color) {
            grid.valueAtIndex(p) = color;
            pqueue.enqueue(p+QPoint(-1,0));
            pqueue.enqueue(p+QPoint(1,0));
            pqueue.enqueue(p+QPoint(0,-1));
            pqueue.enqueue(p+QPoint(0,1));
            pqueue.enqueue(p+QPoint(1,1));
            pqueue.enqueue(p+QPoint(1,-1));
            pqueue.enqueue(p+QPoint(-1,1));
            pqueue.enqueue(p+QPoint(-1,-1));
            ++found;
        }
    }
    return found;
}

// find all neigbors of color 'stand_id' and save in the 'neighbors' vector
int ActSalvage::neighborFinderHelper(Grid<int> &grid, QVector<int> &neighbors, int stand_id)
{
    GridRunner<int> id_runner(&grid);
    neighbors.clear();
    int *nb[8];
    while (id_runner.next()) {
        if (*id_runner.current()==stand_id) {
            id_runner.neighbors8(nb);
            for (int i=0;i<8;++i)
                if (nb[i] && *nb[i]!=-1 && *nb[i]!=stand_id) {
                    if (!neighbors.contains(*nb[i]))
                        neighbors.push_back(*nb[i]);
                }
        }
    }
    return neighbors.size();
}

int ActSalvage::replaceValueHelper(Grid<int> &grid, int old_value, int new_value)
{
    int n=0;
    for (int *p=grid.begin(); p!=grid.end(); ++p)
        if (*p == old_value) {
            *p = new_value;
            ++n;
        }
    return n;
}


} // namespace