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 "spatialanalysis.h"

#include "globalsettings.h"

#include "model.h"

#include "tree.h"

#include "stamp.h"

#include "helper.h"

#include "resourceunit.h"

#include "scriptgrid.h"

#include <QJSEngine>

#include <QJSValue>

void SpatialAnalysis::addToScriptEngine()

{

    SpatialAnalysis *spati = new SpatialAnalysis;

    QJSValue v = GlobalSettings::instance()->scriptEngine()->newQObject(spati);

    GlobalSettings::instance()->scriptEngine()->globalObject().setProperty("SpatialAnalysis", v);

}

SpatialAnalysis::~SpatialAnalysis()

{

    if (mRumple)

        delete mRumple;

}

double SpatialAnalysis::rumpleIndexFullArea()

{

    if (!mRumple)

        mRumple = new RumpleIndex;

    double rum = mRumple->value();

    return rum;

}

/// extract patches (clumps) from the grid 'src'.

/// Patches are defined as adjacent pixels (8-neighborhood)

/// Return: vector with number of pixels per patch (first element: patch 1, second element: patch 2, ...)

QList<int> SpatialAnalysis::extractPatches(Grid<double> &src, int min_size, QString fileName)

{

    mClumpGrid.setup(src.metricRect(), src.cellsize());

    mClumpGrid.wipe();

    // now loop over all pixels and run a floodfill algorithm

    QPoint start;

    QQueue<QPoint> pqueue; // for the flood fill algorithm

    QList<int> counts;

    int patch_index = 0;

    int total_size = 0;

    int patches_skipped = 0;

    for (int i=0;i<src.count();++i) {

        if (src[i]>0. && mClumpGrid[i]==0) {

            start = src.indexOf(i);

            pqueue.clear();

            patch_index++;

            // quick and dirty implementation of the flood fill algroithm.

            // based on: http://en.wikipedia.org/wiki/Flood_fill

            // returns the number of pixels colored

            pqueue.enqueue(start);

            int found = 0;

            while (!pqueue.isEmpty()) {

                QPoint p = pqueue.dequeue();

                if (!src.isIndexValid(p))

                    continue;

                if (src.valueAtIndex(p)>0. && mClumpGrid.valueAtIndex(p) == 0) {

                    mClumpGrid.valueAtIndex(p) = patch_index;

                    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;

                }

            }

            if (found<min_size) {

                // delete the patch again

                pqueue.enqueue(start);

                while (!pqueue.isEmpty()) {

                    QPoint p = pqueue.dequeue();

                    if (!src.isIndexValid(p))

                        continue;

                    if (mClumpGrid.valueAtIndex(p) == patch_index) {

                        mClumpGrid.valueAtIndex(p) = -1;

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

                    }

                }

                --patch_index;

                patches_skipped++;

            } else {

                // save the patch in the result

                counts.push_back(found);

                total_size+=found;

            }

        }

    }

    // remove the -1 again...

    mClumpGrid.limit(0,999999);

    qDebug() << "extractPatches: found" << patch_index << "patches, total valid pixels:" << total_size << "skipped" << patches_skipped;

    if (!fileName.isEmpty()) {

        qDebug() << "extractPatches: save to file:" << GlobalSettings::instance()->path(fileName);

        Helper::saveToTextFile(GlobalSettings::instance()->path(fileName), gridToESRIRaster(mClumpGrid) );

    }

    return counts;

}

void SpatialAnalysis::saveRumpleGrid(QString fileName)

{

    if (!mRumple)

        mRumple = new RumpleIndex;

    Helper::saveToTextFile(GlobalSettings::instance()->path(fileName), gridToESRIRaster(mRumple->rumpleGrid()) );

}

void SpatialAnalysis::saveCrownCoverGrid(QString fileName)

{

    calculateCrownCover();

    Helper::saveToTextFile(GlobalSettings::instance()->path(fileName), gridToESRIRaster(mCrownCoverGrid) );

}

QJSValue SpatialAnalysis::patches(QJSValue grid, int min_size)

{

    ScriptGrid *sg = qobject_cast<ScriptGrid*>(grid.toQObject());

    if (sg) {

        // extract patches (keep patches with a size >= min_size

        mLastPatches = extractPatches(*sg->grid(), min_size, QString());

        // create a (double) copy of the internal clump grid, and return this grid

        // as a JS value

        QJSValue v = ScriptGrid::createGrid(mClumpGrid.toDouble(),"patch");

        return v;

    }

    return QJSValue();

}

void SpatialAnalysis::calculateCrownCover()

{

    mCrownCoverGrid.setup(GlobalSettings::instance()->model()->RUgrid().metricRect(),

                          GlobalSettings::instance()->model()->RUgrid().cellsize());

    // calculate the crown cover per resource unit. We use the "reader"-stamps of the individual trees

    // as they represent the crown (size). We also simply hijack the LIF grid for our calculations.

    FloatGrid *grid = GlobalSettings::instance()->model()->grid();

    grid->initialize(0.f);

    // we simply iterate over all trees of all resource units (not bothering about multithreading here)

    int x,y;

    AllTreeIterator ati(GlobalSettings::instance()->model());

    while (Tree *t = ati.nextLiving()) {

        // apply the reader-stamp

        const Stamp *reader = t->stamp()->reader();

        QPoint pos_reader = t->positionIndex(); // tree position

        pos_reader-=QPoint(reader->offset(), reader->offset());

        int reader_size = reader->size();

        int rx = pos_reader.x();

        int ry = pos_reader.y();

        // the reader stamps are stored such as to have a sum of 1.0 over all pixels

        // (i.e.: they express the percentage for each cell contributing to the full crown).

        // we thus calculate a the factor to "blow up" cell values; a fully covered cell has then a value of 1,

        // and values between 0-1 are cells that are partially covered by the crown.

        double crown_factor = reader->crownArea()/double(cPxSize*cPxSize);

        // add the reader-stamp values: multiple (partial) crowns can add up to being fully covered

        for (y=0;y<reader_size; ++y) {

            for (x=0;x<reader_size;++x) {

                 grid->valueAtIndex(rx+x, ry+y) += (*reader)(x,y)*crown_factor;

            }

        }

    }

    // now aggregate values for each resource unit

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

    for (float *rg = mCrownCoverGrid.begin(); rg!=mCrownCoverGrid.end();++rg) {

        ResourceUnit *ru =  model->RUgrid().constValueAtIndex(mCrownCoverGrid.indexOf(rg));

        if (!ru) {

            *rg=0.f;

            continue;

        }

        float cc_sum = 0.f;

        GridRunner<float> runner(grid, mCrownCoverGrid.cellRect(mCrownCoverGrid.indexOf(rg)));

        while (float *gv = runner.next()) {

            if (model->heightGridValue(runner.currentIndex().x(), runner.currentIndex().y()).isValid())

                if (*gv >= 0.5f) // 0.5: half of a 2m cell is covered by a tree crown; is a bit pragmatic but seems reasonable (and works)

                    cc_sum++;

        }

        if (ru->stockableArea()>0.) {

            double value = cPxSize*cPxSize*cc_sum/ru->stockableArea();

            *rg = limit(value, 0., 1.);

        }

    }

}

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

/********************************* RumpleIndex class ************************************/

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

void RumpleIndex::setup()

{

    mRumpleGrid.clear();

    if (!GlobalSettings::instance()->model()) return;

    // the rumple grid hast the same dimensions as the resource unit grid (i.e. 100 meters)

    mRumpleGrid.setup(GlobalSettings::instance()->model()->RUgrid().metricRect(),

                    GlobalSettings::instance()->model()->RUgrid().cellsize());

}

void RumpleIndex::calculate()

{

    if (mRumpleGrid.isEmpty())

        setup();

    mRumpleGrid.initialize(0.f);

    HeightGrid *hg = GlobalSettings::instance()->model()->heightGrid();

    // iterate over the resource units and calculate the rumple index / surface area for each resource unit

    HeightGridValue* hgv_8[8]; // array holding pointers to height grid values (neighborhood)

    float heights[9];  // array holding heights (8er neighborhood + center pixel)

    int total_valid_pixels = 0;

    float total_surface_area = 0.f;

    for (float *rg = mRumpleGrid.begin(); rg!=mRumpleGrid.end();++rg) {

        int valid_pixels = 0;

        float surface_area_sum = 0.f;

        GridRunner<HeightGridValue> runner(hg, mRumpleGrid.cellRect(mRumpleGrid.indexOf(rg)));

        while (runner.next()) {

            if (runner.current()->isValid()) {

                runner.neighbors8(hgv_8);

                bool valid = true;

                float *hp = heights;

                *hp++ = runner.current()->height;

                // retrieve height values from the grid

                for (int i=0;i<8;++i) {

                    *hp++ = hgv_8[i] ? hgv_8[i]->height: 0;

                    if (hgv_8[i] && !hgv_8[i]->isValid()) valid = false;

                    if (!hgv_8[i]) valid = false;

                }

                // calculate surface area only for cells which are (a) within the project area, and (b) all neighboring pixels are inside the forest area

                if (valid) {

                   valid_pixels++;

                   float surface_area = calculateSurfaceArea(heights, hg->cellsize());

                   surface_area_sum += surface_area;

                }

            }

        }

        if (valid_pixels>0) {

            float rumple_index = surface_area_sum / ( (float) valid_pixels * hg->cellsize()*hg->cellsize());

            *rg = rumple_index;

            total_valid_pixels += valid_pixels;

            total_surface_area += surface_area_sum;

        }

    }

    mRumpleIndex = 0.;

    if (total_valid_pixels>0) {

        float rumple_index = total_surface_area / ( (float) total_valid_pixels * hg->cellsize()*hg->cellsize());

        mRumpleIndex = rumple_index;

    }

    mLastYear = GlobalSettings::instance()->currentYear();

}

double RumpleIndex::value(const bool force_recalculate)

{

    if (force_recalculate ||  mLastYear != GlobalSettings::instance()->currentYear())

        calculate();

    return mRumpleIndex;

}

double RumpleIndex::test_triangle_area()

{

    // check calculation: numbers for Jenness paper

    float hs[9]={165, 170, 145, 160, 183, 155,122,175,190};

    double area = calculateSurfaceArea(hs, 100);

    return area;

}

inline double surface_length(float h1, float h2, float l)

{

    return sqrt((h1-h2)*(h1-h2) + l*l);

}

inline double heron_triangle_area(float a, float b, float c)

{

    float s=(a+b+c)/2.f;

    return sqrt(s * (s-a) * (s-b) * (s-c) );

}

/// calculate the surface area of a pixel given its height value, the height of the 8 neigboring pixels, and the cellsize

/// the algorithm is based on http://www.jennessent.com/downloads/WSB_32_3_Jenness.pdf

double RumpleIndex::calculateSurfaceArea(const float *heights, const float cellsize)

{

    // values in the height array [0..8]: own height / north/east/west/south/ NE/NW/SE/SW

    // step 1: calculate length on 3d surface between all edges

    //   8(A) * 1(B) * 5(C)       <- 0: center cell, indices in the "heights" grid, A..I: codes used by Jenness

    //   4(D) * 0(E) * 2(F)

    //   7(G) * 3(H) * 6(I)

    float slen[16]; // surface lengths (divided by 2)

    // horizontal

    slen[0] = surface_length(heights[8], heights[1], cellsize)/2.f;

    slen[1] = surface_length(heights[1], heights[5], cellsize)/2.f;

    slen[2] = surface_length(heights[4], heights[0], cellsize)/2.f;

    slen[3] = surface_length(heights[0], heights[2], cellsize)/2.f;

    slen[4] = surface_length(heights[7], heights[3], cellsize)/2.f;

    slen[5] = surface_length(heights[3], heights[6], cellsize)/2.f;

    // vertical

    slen[6] = surface_length(heights[8], heights[4], cellsize)/2.f;

    slen[7] = surface_length(heights[1], heights[0], cellsize)/2.f;

    slen[8] = surface_length(heights[5], heights[2], cellsize)/2.f;

    slen[9] = surface_length(heights[4], heights[7], cellsize)/2.f;

    slen[10] = surface_length(heights[0], heights[3], cellsize)/2.f;

    slen[11] = surface_length(heights[2], heights[6], cellsize)/2.f;

    // diagonal

    float cellsize_diag = cellsize * M_SQRT2;

    slen[12] = surface_length(heights[0], heights[8], cellsize_diag)/2.f;

    slen[13] = surface_length(heights[0], heights[5], cellsize_diag)/2.f;

    slen[14] = surface_length(heights[0], heights[7], cellsize_diag)/2.f;

    slen[15] = surface_length(heights[0], heights[6], cellsize_diag)/2.f;

    // step 2: combine the three sides of all the 8 sub triangles using Heron's formula

    double surface_area = 0.;

    surface_area += heron_triangle_area(slen[12], slen[0], slen[7]); // i

    surface_area += heron_triangle_area(slen[7], slen[1], slen[13]); // ii

    surface_area += heron_triangle_area(slen[6], slen[2], slen[12]); // iii

    surface_area += heron_triangle_area(slen[13], slen[8], slen[3]); // iv

    surface_area += heron_triangle_area(slen[2], slen[9], slen[14]); // v

    surface_area += heron_triangle_area(slen[3], slen[11], slen[15]); // vi

    surface_area += heron_triangle_area(slen[14], slen[10], slen[4]); // vii

    surface_area += heron_triangle_area(slen[10], slen[15], slen[5]); // viii

    return surface_area;

}

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

/* ******************************* Spatial Layers **************************************** */

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

void SpatialLayeredGrid::setup()

{

    addGrid("rumple", 0);

}

void SpatialLayeredGrid::createGrid(const int grid_index)

{

    Q_UNUSED(grid_index); // TODO: what should happen here?

}

int SpatialLayeredGrid::addGrid(const QString name, FloatGrid *grid)

{

    mGridNames.append(name);

    mGrids.append(grid);

    return mGrids.size();

}