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 "global.h"
#include "dem.h"

#include "globalsettings.h"
#include "model.h"

#include "gisgrid.h"

// from here: http://www.scratchapixel.com/lessons/3d-advanced-lessons/interpolation/bilinear-interpolation/
template<typename T>
T bilinear(
   const T &tx,
   const T &ty,
   const T &c00,
   const T &c10,
   const T &c01,
   const T &c11)
{
#if 1
    T a = c00 * (T(1) - tx) + c10 * tx;
    T b = c01 * (T(1) - tx) + c11 * tx;
    return a * (T(1) - ty) + b * ty;
#else
    return (T(1) - tx) * (T(1) - ty) * c00 +
        tx * (T(1) - ty) * c10 +
        (T(1) - tx) * ty * c01 +
        tx * ty * c11;
#endif
}

/// loads a DEM from a ESRI style text file.
/// internally, the DEM has always a resolution of 10m
bool DEM::loadFromFile(const QString &fileName)
{
    if (!GlobalSettings::instance()->model())
        throw IException("DEM::create10mGrid: no valid model to retrieve height grid.");

    HeightGrid *h_grid = GlobalSettings::instance()->model()->heightGrid();
    if (!h_grid || h_grid->isEmpty())
        throw IException("GisGrid::create10mGrid: no valid height grid to copy grid size.");

    GisGrid gis_grid;
    if (!gis_grid.loadFromFile(fileName))
        throw IException(QString("Unable to load DEM file %1").arg(fileName));
    // create a grid with the same size as the height grid
    // (height-grid: 10m size, covering the full extent)
    clear();
    aspect_grid.clear();
    slope_grid.clear();
    view_grid.clear();

    setup(h_grid->metricRect(),h_grid->cellsize());

    const QRectF &world = GlobalSettings::instance()->model()->extent();

    if (fmod(gis_grid.cellSize(), cellsize()) != 0) {
        QPointF p;
        // simple copy of the data
        for (int i=0;i<count();i++) {
            p = cellCenterPoint(indexOf(i));
            if (gis_grid.value(p) != gis_grid.noDataValue() && world.contains(p) )
                valueAtIndex(i) = gis_grid.value(p);
            else
                valueAtIndex(i) = -1;
        }
    } else {
        // bilinear approximation approach
        qDebug() << "DEM: built-in bilinear interpolation from cell size" << gis_grid.cellSize();
        int f = gis_grid.cellSize() / cellsize(); // size-factor
        initialize(-1.f);
        int ixmin = 10000000, iymin = 1000000, ixmax = -1, iymax = -1;
        for (int y=0;y<gis_grid.rows();++y)
            for (int x=0;x<gis_grid.cols(); ++x){
                Vector3D p3d = gis_grid.coord(x,y);
                if (world.contains(p3d.x(), p3d.y())) {
                    QPoint pt=indexAt(QPointF(p3d.x(), p3d.y()));
                    valueAt(p3d.x(), p3d.y()) = gis_grid.value(x,y);
                    ixmin = std::min(ixmin, pt.x()); ixmax=std::max(ixmax, pt.x());
                    iymin = std::min(iymin, pt.y()); iymax=std::max(iymax, pt.y());
                }
            }
        for (int y=iymin;y<=iymax-f;y+=f)
            for (int x=ixmin;x<=ixmax-f;x+=f){
                float c00 = valueAtIndex(x, y);
                float c10 = valueAtIndex(x+f, y);
                float c01 = valueAtIndex(x, y+f);
                float c11 = valueAtIndex(x+f, y+f);
                for (int my=0;my<f;++my)
                    for (int mx=0;mx<f;++mx)
                        valueAtIndex(x+mx, y+my) = bilinear<float>(mx/float(f), my/float(f), c00, c10, c01, c11);
            }
    }

    return true;
}

/// calculate slope and aspect at a given point.
/// results are params per reference.
/// returns the height at point (x/y)
/// calculation follows: Burrough, P. A. and McDonell, R.A., 1998.Principles of Geographical Information Systems.(Oxford University Press, New York), p. 190.
/// http://uqu.edu.sa/files2/tiny_mce/plugins/filemanager/files/4280125/Principles%20of%20Geographical%20Information%20Systems.pdf
/// @param point metric coordinates of point to derive orientation
/// @param rslope_angle RESULTING (passed by reference) slope angle as percentage (i.e: 1:=45 degrees)
/// @param rslope_aspect RESULTING slope direction in degrees (0: North, 90: east, 180: south, 270: west)
float DEM::orientation(const QPointF &point, float &rslope_angle, float &rslope_aspect) const
{
    QPoint pt = indexAt(point);
    if (pt.x()>0 && pt.x()<sizeX()+1 && pt.y()>0 && pt.y()<sizeY()-1) {
        float *p = const_cast<DEM*>(this)->ptr(pt.x(), pt.y());
        float z2 = *(p-sizeX());
        float z4 = *(p-1);
        float z6 = *(p+1);
        float z8 = *(p+sizeX());
        float g = (-z4 + z6) / (2*cellsize());
        float h = (z2 - z8) / (2*cellsize());

        if (z2<=0. || z4<=0. || z6<=0. || z8<=0) {
            rslope_angle = 0.;
            rslope_aspect = 0.;
            return *p;
        }
        rslope_angle = sqrt(g*g + h*h);
        // atan2: returns -pi : +pi
        // North: -pi/2, east: 0, south: +pi/2, west: -pi/+pi
        float aspect = atan2(-h, -g);
        // transform to degree:
        // north: 0, east: 90, south: 180, west: 270
        aspect = aspect * 180. / M_PI + 360. + 90.;
        aspect = fmod(aspect, 360.f);

        rslope_aspect = aspect;
        return *p;
    } else {
        rslope_angle = 0.;
        rslope_aspect = 0.;
        return 0.;
    }
}

void DEM::createSlopeGrid() const
{
    if (slope_grid.isEmpty()) {
        // setup custom grids with the same size as this DEM
        slope_grid.setup(*this);
        view_grid.setup(*this);
        aspect_grid.setup(*this);
    } else {
        return;
    }
    float *slope = slope_grid.begin();
    float *view = view_grid.begin();
    float *aspect = aspect_grid.begin();

    // use fixed values for azimuth (315) and angle (45 deg) and calculate
    // norm vectors
    float sun_x = cos(315. * M_PI/180.) * cos(45.*M_PI/180.);
    float sun_y = sin(315. * M_PI/180.) * cos(45.*M_PI/180.);
    float sun_z = sin(45.*M_PI/180.);

    float a_x, a_y, a_z;
    for (float *p = begin(); p!=end(); ++p, ++slope, ++view, ++aspect) {
        QPointF pt = cellCenterPoint(indexOf(p));
        float height = orientation(pt, *slope, *aspect);
        // calculate the view value:
        if (height>0) {
            float h = atan(*slope);
            a_x = cos(*aspect * M_PI/180.) * cos(h);
            a_y = sin(*aspect * M_PI/180.) * cos(h);
            a_z = sin(h);

            // use the scalar product to calculate the angle, and then
            // transform from [-1,1] to [0,1]
            *view = (a_x*sun_x + a_y*sun_y + a_z*sun_z + 1.)/2.;
        } else {
            *view = 0.;
        }
    }

}