/*
* This file is part of JGrasstools (http://www.jgrasstools.org)
* (C) HydroloGIS - www.hydrologis.com
*
* JGrasstools 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/>.
*/
package org.jgrasstools.hortonmachine.modules.geomorphology.curvatures;
import static java.lang.Math.atan;
import static java.lang.Math.pow;
import static java.lang.Math.sqrt;
import static java.lang.Math.toDegrees;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_AUTHORCONTACTS;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_AUTHORNAMES;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_DRAFT;
import static org.jgrasstools.gears.i18n.GearsMessages.OMSHYDRO_LICENSE;
import static org.jgrasstools.gears.libs.modules.JGTConstants.doubleNovalue;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSCURVATURES_LABEL;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSCURVATURES_outPlan_DESCRIPTION;
import static org.jgrasstools.hortonmachine.i18n.HortonMessages.OMSCURVATURES_outProf_DESCRIPTION;
import java.awt.image.WritableRaster;
import javax.media.jai.iterator.RandomIter;
import oms3.annotations.Author;
import oms3.annotations.Bibliography;
import oms3.annotations.Description;
import oms3.annotations.Execute;
import oms3.annotations.In;
import oms3.annotations.Keywords;
import oms3.annotations.Label;
import oms3.annotations.License;
import oms3.annotations.Name;
import oms3.annotations.Out;
import oms3.annotations.Status;
import org.apache.commons.math3.linear.DecompositionSolver;
import org.apache.commons.math3.linear.MatrixUtils;
import org.apache.commons.math3.linear.RRQRDecomposition;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.geotools.coverage.grid.GridCoverage2D;
import org.jgrasstools.gears.libs.modules.GridNode;
import org.jgrasstools.gears.libs.modules.JGTConstants;
import org.jgrasstools.gears.libs.modules.JGTModel;
import org.jgrasstools.gears.utils.RegionMap;
import org.jgrasstools.gears.utils.coverage.CoverageUtilities;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
@Description("Estimates the longitudinal, normal and planar curvatures by means of a bivariate quadratic representation of the terrain.")
@Author(name = OMSHYDRO_AUTHORNAMES, contact = OMSHYDRO_AUTHORCONTACTS)
@Keywords("curvatures, bivariate, slope, aspect")
@Label(OMSCURVATURES_LABEL)
@Name("oms_curvaturesbivariate")
@Status(OMSHYDRO_DRAFT)
@License(OMSHYDRO_LICENSE)
@Bibliography("Multiscale Terrain Analysis of Multibeam Bathymetry Data for Habitat Mapping on the Continental Slope, Wilson M., 2007")
public class OmsCurvaturesBivariate extends JGTModel {
@Description("The map of the digital elevation model (DEM or pit).")
@In
public GridCoverage2D inElev = null;
@Description("The size of the analysis window in odd cells number.")
@In
public int pCells = 3;
// output
@Description(OMSCURVATURES_outProf_DESCRIPTION)
@Out
public GridCoverage2D outProf = null;
@Description(OMSCURVATURES_outPlan_DESCRIPTION)
@Out
public GridCoverage2D outPlan = null;
@Description("The map of slope.")
@Out
public GridCoverage2D outSlope = null;
@Description("The map of aspect")
@Out
public GridCoverage2D outAspect = null;
@Execute
public void process() throws Exception {
if (!concatOr(outProf == null, doReset)) {
return;
}
checkNull(inElev);
if (pCells < 3) {
pCells = 3;
}
if (pCells % 2 == 0) {
pCells = pCells + 1;
}
RegionMap regionMap = CoverageUtilities.getRegionParamsFromGridCoverage(inElev);
int nCols = regionMap.getCols();
int nRows = regionMap.getRows();
double xRes = regionMap.getXres();
double yRes = regionMap.getYres();
RandomIter elevationIter = CoverageUtilities.getRandomIterator(inElev);
WritableRaster profWR = CoverageUtilities.createDoubleWritableRaster(nCols, nRows, null, null, doubleNovalue);
WritableRaster planWR = CoverageUtilities.createDoubleWritableRaster(nCols, nRows, null, null, doubleNovalue);
WritableRaster slopeWR = CoverageUtilities.createDoubleWritableRaster(nCols, nRows, null, null, doubleNovalue);
WritableRaster aspectWR = CoverageUtilities.createDoubleWritableRaster(nCols, nRows, null, null, doubleNovalue);
final double[] planProfSlopeAspect = new double[4];
double disXX = Math.pow(xRes, 2.0);
double disYY = Math.pow(yRes, 2.0);
/*
* calculate curvatures
*/
pm.beginTask("Processing...", nRows - 2);
for( int r = 1; r < nRows - 1; r++ ) {
if (isCanceled(pm)) {
return;
}
for( int c = 1; c < nCols - 1; c++ ) {
calculateCurvatures(elevationIter, planProfSlopeAspect, nCols, nRows, c, r, xRes, yRes, disXX, disYY, pCells);
planWR.setSample(c, r, 0, planProfSlopeAspect[0]);
profWR.setSample(c, r, 0, planProfSlopeAspect[1]);
slopeWR.setSample(c, r, 0, planProfSlopeAspect[2]);
aspectWR.setSample(c, r, 0, planProfSlopeAspect[3]);
}
pm.worked(1);
}
pm.done();
if (isCanceled(pm)) {
return;
}
CoordinateReferenceSystem crs = inElev.getCoordinateReferenceSystem();
outProf = CoverageUtilities.buildCoverage("prof_curvature", profWR, regionMap, crs);
outPlan = CoverageUtilities.buildCoverage("plan_curvature", planWR, regionMap, crs);
outSlope = CoverageUtilities.buildCoverage("slope", slopeWR, regionMap, crs);
outAspect = CoverageUtilities.buildCoverage("aspect", aspectWR, regionMap, crs);
}
/**
* Calculate curvatures for a single cell.
*
* @param elevationIter the elevation map.
* @param planTangProf the array into which to insert the resulting [plan, tang, prof] curvatures.
* @param col the column to process.
* @param row the row to process.
* @param ncols the columns of the raster.
* @param nrows the rows of the raster.
* @param xRes
* @param yRes
* @param disXX the diagonal size of the cell, x component.
* @param disYY the diagonal size of the cell, y component.
*/
public static void calculateCurvatures( RandomIter elevationIter, final double[] planTangProf, int ncols, int nrows, int col,
int row, double xRes, double yRes, double disXX, double disYY, int windowSize ) {
GridNode node = new GridNode(elevationIter, ncols, nrows, xRes, yRes, col, row);
double[][] window = node.getWindow(windowSize, false);
if (!hasNovalues(window)) {
double[] parameters = calculateParameters(window);
double a = parameters[0];
double b = parameters[1];
double c = parameters[2];
double d = parameters[3];
double e = parameters[4];
double slope = atan(sqrt(d * d + e * e));
slope = toDegrees(slope);
double aspect = atan(e / d);
aspect = toDegrees(aspect);
double profcNumerator = -200.0 * (a * d * d + b * e * e + c * d * e);
double profcDenominator = (e * e + d * d) * pow((1 + e * e + d * d), 1.5);
double profc = profcNumerator / profcDenominator;
double plancNumerator = -200.0 * (b * d * d + a * e * e + c * d * e);
double plancDenominator = pow((e * e + d * d), 1.5);
double planc = plancNumerator / plancDenominator;
planTangProf[0] = planc;
planTangProf[1] = profc;
planTangProf[2] = slope;
planTangProf[3] = aspect;
} else {
planTangProf[0] = doubleNovalue;
planTangProf[1] = doubleNovalue;
planTangProf[2] = doubleNovalue;
planTangProf[3] = doubleNovalue;
}
}
private static boolean hasNovalues( double[][] window ) {
for( int i = 0; i < window.length; i++ ) {
for( int j = 0; j < window[0].length; j++ ) {
if (JGTConstants.isNovalue(window[i][j])) {
return true;
}
}
}
return false;
}
/**
* Calculates the parameters of a bivariate quadratic equation.
*
* @param elevationValues the window of points to use.
* @return the parameters of the bivariate quadratic equation as [a, b, c, d, e, f]
*/
private static double[] calculateParameters( final double[][] elevationValues ) {
int rows = elevationValues.length;
int cols = elevationValues[0].length;
int pointsNum = rows * cols;
final double[][] xyMatrix = new double[pointsNum][6];
final double[] valueArray = new double[pointsNum];
// TODO check on resolution
int index = 0;
for( int y = 0; y < rows; y++ ) {
for( int x = 0; x < cols; x++ ) {
xyMatrix[index][0] = x * x; // x^2
xyMatrix[index][1] = y * y; // y^2
xyMatrix[index][2] = x * y; // xy
xyMatrix[index][3] = x; // x
xyMatrix[index][4] = y; // y
xyMatrix[index][5] = 1;
valueArray[index] = elevationValues[y][x];
index++;
}
}
RealMatrix A = MatrixUtils.createRealMatrix(xyMatrix);
RealVector z = MatrixUtils.createRealVector(valueArray);
DecompositionSolver solver = new RRQRDecomposition(A).getSolver();
RealVector solution = solver.solve(z);
// start values for a, b, c, d, e, f, all set to 0.0
final double[] parameters = solution.toArray();
return parameters;
}
}