org.jgrasstools.lesto.modules.raster.Las2BivariateRasterMosaic Maven / Gradle / Ivy
The newest version!
/*
* 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 tileLasPoints = lasData.getPointsInGeometry(roiPolygon, true);
if (tileLasPoints.size() > 100) {
final GridGeometry2D gridGeometry = outputCoverage.getGridGeometry();
int newCols = cols + 2 * pBuffer;
int newRows = rows + 2 * pBuffer;
GridGeometry2D bufferedGridGeometry = CoverageUtilities.gridGeometryFromRegionValues(env.getMaxY(), env.getMinY(),
env.getMaxX(), env.getMinX(), newCols, newRows, crs);
ArrayList[][] lasMatrix = new ArrayList[newCols][newRows];
for( int c = 0; c < newCols; c++ ) {
for( int r = 0; r < newRows; r++ ) {
ArrayList item = new ArrayList();
lasMatrix[c][r] = item;
}
}
// Splitting las into cells
final Point point = new Point();
for( LasRecord dot : tileLasPoints ) {
if (doIntensity) {
if (dot.intensity == NOINTENSITY) {
continue;
}
minValue = min(dot.intensity, minValue);
maxValue = max(dot.intensity, maxValue);
} else {
minValue = min(dot.z, minValue);
maxValue = max(dot.z, maxValue);
}
CoverageUtilities.colRowFromCoordinate(new Coordinate(dot.x, dot.y), bufferedGridGeometry, point);
lasMatrix[point.x][point.y].add(dot);
}
WritableRandomIter outWIter = null;
try {
outWIter = CoverageUtilities.getWritableRandomIterator(outWR);
for( int c = pBuffer; c < newCols - pBuffer; c++ ) {
if (c % 100 == 0) {
StringBuilder sb1 = new StringBuilder();
sb1.append(tileName);
sb1.append(": ");
sb1.append(c);
sb1.append(" of ");
sb1.append(newCols);
pm.message(sb1.toString());
}
for( int r = pBuffer; r < newRows - pBuffer; r++ ) {
Coordinate coordinate = CoverageUtilities.coordinateFromColRow(c, r, bufferedGridGeometry);
List currentPoints = new ArrayList();
for( int tmpC = c - pBuffer; tmpC <= c + pBuffer; tmpC++ ) {
for( int tmpR = r - pBuffer; tmpR <= r + pBuffer; tmpR++ ) {
currentPoints.addAll(lasMatrix[tmpC][tmpR]);
}
}
int size = currentPoints.size();
if (size >= pMinpoints) {
// need at least as many samples as parameters
try {
double[] parameters = calculateParameters(currentPoints);
double interpolatedValue = getInterpolatedValue(parameters, coordinate.x, coordinate.y);
// limit by min/max
if (interpolatedValue < minValue) {
interpolatedValue = minValue;
}
if (interpolatedValue > maxValue) {
interpolatedValue = maxValue;
}
CoverageUtilities.colRowFromCoordinate(coordinate, gridGeometry, point);
outWIter.setSample(point.x, point.y, 0, interpolatedValue);
} catch (SingularMatrixException ex) {
// we ignore the singular matrix leaving the cell undefined
}
}
}
}
} finally {
if (outWIter != null)
outWIter.done();
}
dumpRaster(outputCoverage, outTileFile.getAbsolutePath());
} else {
pm.message(tileName + " ignored because of no points there.");
}
}
private double[] calculateParameters( final List pointsInGeometry ) {
int pointsNum = pointsInGeometry.size();
final double[][] xyMatrix = new double[pointsNum][6];
final double[] valueArray = new double[pointsNum];
for( int i = 0; i < pointsNum; i++ ) {
LasRecord dot = pointsInGeometry.get(i);
xyMatrix[i][0] = dot.x * dot.x; // x^2
xyMatrix[i][1] = dot.y * dot.y; // y^2
xyMatrix[i][2] = dot.x * dot.y; // xy
xyMatrix[i][3] = dot.x; // x
xyMatrix[i][4] = dot.y; // y
xyMatrix[i][5] = 1;
if (doIntensity) {
valueArray[i] = dot.intensity;
} else {
valueArray[i] = dot.z;
}
}
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;
}
private double getInterpolatedValue( double[] parameters, double x, double y ) {
double z = parameters[0] * x * x + //
parameters[1] * y * y + //
parameters[2] * x * y + //
parameters[3] * x + //
parameters[4] * y + //
parameters[5];
return z;
}
}
© 2015 - 2025 Weber Informatics LLC | Privacy Policy