Skip to content

Commit

Permalink
[SEDONA-404] Add RS_Resample (#1038)
Browse files Browse the repository at this point in the history
  • Loading branch information
@iGN5117
iGN5117 authored Sep 29, 2023
1 parent db07953 commit b135e8e
Show file tree
Hide file tree
Showing 7 changed files with 477 additions and 2 deletions.
133 changes: 133 additions & 0 deletions common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -25,14 +25,23 @@
import org.geotools.coverage.grid.GridCoverageFactory;
import org.geotools.coverage.grid.GridEnvelope2D;
import org.geotools.coverage.grid.GridGeometry2D;
import org.geotools.coverage.processing.Operations;
import org.geotools.geometry.Envelope2D;
import org.geotools.referencing.crs.DefaultEngineeringCRS;
import org.geotools.referencing.operation.transform.AffineTransform2D;
import org.opengis.coverage.grid.GridCoverage;
import org.opengis.metadata.spatial.PixelOrientation;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.MathTransform2D;
import org.opengis.referencing.operation.TransformException;

import javax.media.jai.Interpolation;
import java.awt.geom.Point2D;
import java.util.Map;
import java.util.Objects;

public class RasterEditors
{
Expand Down Expand Up @@ -97,4 +106,128 @@ public static GridCoverage2D setGeoReference(GridCoverage2D raster, double upper
String geoRedCoord = String.format("%f %f %f %f %f %f", scaleX, skewY, skewX, scaleY, upperLeftX, upperLeftY);
return setGeoReference(raster, geoRedCoord, "GDAL");
}

public static GridCoverage2D resample(GridCoverage2D raster, double widthOrScale, double heightOrScale, double gridX, double gridY, boolean useScale,String algorithm) throws TransformException {
/*
* Old Parameters
*/
AffineTransform2D affine = RasterUtils.getGDALAffineTransform(raster);
int originalWidth = RasterAccessors.getWidth(raster), originalHeight = RasterAccessors.getHeight(raster);
double upperLeftX = affine.getTranslateX(), upperLeftY = affine.getTranslateY();
double originalSkewX = affine.getShearX(), originalSkewY = affine.getShearY();
double originalScaleX = affine.getScaleX(), originalScaleY = affine.getScaleY();
CoordinateReferenceSystem crs = raster.getCoordinateReferenceSystem2D();

/*
* New Parameters
*/
int newWidth = useScale ? originalWidth : (int)Math.floor(widthOrScale);
int newHeight = useScale ? originalHeight : (int)Math.floor(heightOrScale);
double newScaleX = useScale ? widthOrScale : originalScaleX;
double newScaleY = useScale ? heightOrScale : originalScaleY;
double newUpperLeftX = upperLeftX, newUpperLeftY = upperLeftY;

if (noConfigChange(originalWidth, originalHeight, upperLeftX, upperLeftY, originalScaleX, originalScaleY, newWidth, newHeight, gridX, gridY, newScaleX, newScaleY, useScale)) {
// no reconfiguration parameters provided
return raster;
}


Envelope2D envelope2D = raster.getEnvelope2D();
//process scale changes due to changes in widthOrScale and heightOrScale
if (!useScale) {
newScaleX = (Math.abs(envelope2D.getMaxX() - envelope2D.getMinX())) / newWidth;
newScaleY = Math.signum(originalScaleY) * Math.abs(envelope2D.getMaxY() - envelope2D.getMinY()) / newHeight;
}else {
//height and width cannot have floating point, ceil them to next greatest integer in that case.
newWidth = (int) Math.ceil(Math.abs(envelope2D.getMaxX() - envelope2D.getMinX()) / Math.abs(newScaleX));
newHeight = (int) Math.ceil(Math.abs(envelope2D.getMaxY() - envelope2D.getMinY()) / Math.abs(newScaleY));
}

if (!approximateEquals(upperLeftX, gridX) || !approximateEquals(upperLeftY, gridY)) {
//change upperLefts to gridX/Y to check if any warping is needed
GridCoverage2D tempRaster = setGeoReference(raster, gridX, gridY, newScaleX, newScaleY, originalSkewX, originalSkewY);

//check expected grid coordinates for old upperLefts
int[] expectedCellCoordinates = RasterUtils.getGridCoordinatesFromWorld(tempRaster, upperLeftX, upperLeftY);

//get expected world coordinates at the expected grid coordinates
Point2D expectedGeoPoint = RasterUtils.getWorldCornerCoordinates(tempRaster, expectedCellCoordinates[0] + 1, expectedCellCoordinates[1] + 1);

//check for shift
if (!approximateEquals(expectedGeoPoint.getX(), upperLeftX)) {
if (!useScale) {
newScaleX = Math.abs(envelope2D.getMaxX() - expectedGeoPoint.getX()) / newWidth;
}else {
//width cannot have floating point, ceil it to next greatest integer in that case.
newWidth = (int) Math.ceil(Math.abs(envelope2D.getMaxX() - expectedGeoPoint.getX()) / Math.abs(newScaleX));
}
newUpperLeftX = expectedGeoPoint.getX();
}

if (!approximateEquals(expectedGeoPoint.getY(), upperLeftY)) {
if (!useScale) {
newScaleY = Math.signum(newScaleY) * Math.abs(envelope2D.getMinY() - expectedGeoPoint.getY()) / newHeight;
}else {
//height cannot have floating point, ceil it to next greatest integer in that case.
newHeight = (int) Math.ceil(Math.abs(envelope2D.getMinY() - expectedGeoPoint.getY()) / Math.abs(newScaleY));
}
newUpperLeftY = expectedGeoPoint.getY();
}
}

MathTransform transform = new AffineTransform2D(newScaleX, originalSkewY, originalSkewX, newScaleY, newUpperLeftX, newUpperLeftY);
GridGeometry2D gridGeometry = new GridGeometry2D(
new GridEnvelope2D(0, 0, newWidth, newHeight),
PixelInCell.CELL_CORNER,
transform, crs, null);
Interpolation resamplingAlgorithm = Interpolation.getInstance(0);
if (!Objects.isNull(algorithm) && !algorithm.isEmpty()) {
if (algorithm.equalsIgnoreCase("nearestneighbor")) {
resamplingAlgorithm = Interpolation.getInstance(0);
}else if (algorithm.equalsIgnoreCase("bilinear")) {
resamplingAlgorithm = Interpolation.getInstance(1);
}else if (algorithm.equalsIgnoreCase("bicubic")) {
resamplingAlgorithm = Interpolation.getInstance(2);
}
}
GridCoverage2D newRaster = (GridCoverage2D) Operations.DEFAULT.resample(raster, null, gridGeometry, resamplingAlgorithm);
return newRaster;
}
public static GridCoverage2D resample(GridCoverage2D raster, double widthOrScale, double heightOrScale, boolean useScale, String algorithm) throws TransformException {
return resample(raster, widthOrScale, heightOrScale, RasterAccessors.getUpperLeftX(raster), RasterAccessors.getUpperLeftY(raster), useScale, algorithm);

}

public static GridCoverage2D resample(GridCoverage2D raster, GridCoverage2D referenceRaster, boolean useScale, String algorithm) throws FactoryException, TransformException {
int srcSRID = RasterAccessors.srid(raster);
int destSRID = RasterAccessors.srid(referenceRaster);
if (srcSRID != destSRID) {
throw new IllegalArgumentException("Provided input raster and reference raster have different SRIDs");
}
double[] refRasterMetadata = RasterAccessors.metadata(referenceRaster);
int newWidth = (int) refRasterMetadata[2];
int newHeight = (int) refRasterMetadata[3];
double gridX = refRasterMetadata[0];
double gridY = refRasterMetadata[1];
double newScaleX = refRasterMetadata[4];
double newScaleY = refRasterMetadata[5];
if (useScale) {
return resample(raster, newScaleX, newScaleY, gridX, gridY, useScale, algorithm);
}
return resample(raster, newWidth, newHeight, gridX, gridY, useScale, algorithm);
}

private static boolean approximateEquals(double a, double b) {
double tolerance = 1E-6;
return Math.abs(a - b) <= tolerance;
}

private static boolean noConfigChange(int oldWidth, int oldHeight, double oldUpperX, double oldUpperY, double originalScaleX, double originalScaleY,
int newWidth, int newHeight, double newUpperX, double newUpperY, double newScaleX, double newScaleY, boolean useScale) {
if (!useScale)
return ((oldWidth == newWidth) && (oldHeight == newHeight) && (approximateEquals(oldUpperX, newUpperX)) && (approximateEquals(oldUpperY, newUpperY)));
return ((approximateEquals(originalScaleX, newScaleX)) && (approximateEquals(originalScaleY, newScaleY)) && (approximateEquals(oldUpperX, newUpperX)) && (approximateEquals(oldUpperY, newUpperY)));
}

}
179 changes: 179 additions & 0 deletions common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,11 +21,13 @@
import org.geotools.coverage.grid.GridCoverage2D;
import org.junit.Test;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.operation.TransformException;

import java.io.IOException;
import java.util.Arrays;

import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertThrows;

public class RasterEditorsTest extends RasterTestBase {
@Test
Expand Down Expand Up @@ -87,4 +89,181 @@ public void testSetGeoReferenceWithEmptyRasterSRID() throws FactoryException {
expected = "10.000000 \n3.000000 \n3.000000 \n-10.000000 \n15.000000 \n-7.000000";
assertEquals(expected, actual);
}


@Test
public void testResample() throws FactoryException, TransformException {
double[] values = {1, 2, 3, 5, 4, 5, 6, 9, 7, 8, 9, 10};
GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "d", 4, 3, 0, 0, 2, -2, 0, 0, 0);
raster = MapAlgebra.addBandFromArray(raster, values, 1, null);

//test with height and width
GridCoverage2D newRaster = RasterEditors.resample(raster, 6, 5, 1, -1, false, "nearestneighbor");
String res = RasterOutputs.asMatrix(newRaster);
String expectedRes = "| 1.0 1.0 2.0 3.0 3.0 5.0|\n" +
"| 1.0 1.0 2.0 3.0 3.0 5.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
double[] metadata = RasterAccessors.metadata(newRaster);
double[] expectedMetadata = {-0.33333333333333326,0.19999999999999996,6,5,1.388888888888889,-1.24,0,0,0,1};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}

//test with scaleX and scaleY
newRaster = RasterEditors.resample(raster, 1.2, -1.4, 1, -1, true, null);
res = RasterOutputs.asMatrix(newRaster);
expectedRes = "| 1.0 1.0 2.0 3.0 3.0 5.0 5.0|\n" +
"| 1.0 1.0 2.0 3.0 3.0 5.0 5.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0 9.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0 10.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0 10.0|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
metadata = RasterAccessors.metadata(newRaster);
expectedMetadata = new double[]{-0.20000000298023224, 0.4000000059604645, 7.0, 5.0, 1.2, -1.4, 0.0, 0.0, 0.0, 1.0};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}
}

@Test
public void testResampleResizeFlavor() throws FactoryException, TransformException {
double[] values = {1, 2, 3, 5, 4, 5, 6, 9, 7, 8, 9, 10};
GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "d", 4, 3, 0, 0, 2, -2, 0, 0, 0);
raster = MapAlgebra.addBandFromArray(raster, values, 1, null);
GridCoverage2D newRaster = RasterEditors.resample(raster, 6, 5, false, "nearestneighbor");
String res = RasterOutputs.asMatrix(newRaster);
String expectedRes = "| 1.0 2.0 2.0 3.0 5.0 5.0|\n" +
"| 1.0 2.0 2.0 3.0 5.0 5.0|\n" +
"| 4.0 5.0 5.0 6.0 9.0 9.0|\n" +
"| 7.0 8.0 8.0 9.0 10.0 10.0|\n" +
"| 7.0 8.0 8.0 9.0 10.0 10.0|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
double[] metadata = RasterAccessors.metadata(newRaster);
double[] expectedMetadata = {0,0,6,5,1.3333333333333333,-1.2,0,0,0,1};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}

//check with scaleX and scaleY
newRaster = RasterEditors.resample(raster, 1.2, -1.4, true, null);
res = RasterOutputs.asMatrix(newRaster);
expectedRes = "| 1.0 1.0 2.0 3.0 3.0 5.0 5.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0 9.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0 9.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0 10.0|\n" +
"| NaN NaN NaN NaN NaN NaN NaN|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
metadata = RasterAccessors.metadata(newRaster);
expectedMetadata = new double[]{0,0,7,5,1.2,-1.4,0,0,0,1};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}
}


@Test
public void testResampleRefRaster() throws FactoryException, TransformException {
double[] values = {1, 2, 3, 5, 4, 5, 6, 9, 7, 8, 9, 10};
GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "d", 4, 3, 0, 0, 2, -2, 0, 0, 0);
GridCoverage2D refRaster = RasterConstructors.makeEmptyRaster(2, "d", 6, 5, 1, -1, 1.2, -1.4, 0, 0, 0);
raster = MapAlgebra.addBandFromArray(raster, values, 1, null);

//test with height and width
GridCoverage2D newRaster = RasterEditors.resample(raster, refRaster, false, null);
String res = RasterOutputs.asMatrix(newRaster);
String expectedRes = "| 1.0 1.0 2.0 3.0 3.0 5.0|\n" +
"| 1.0 1.0 2.0 3.0 3.0 5.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
double[] metadata = RasterAccessors.metadata(newRaster);
double[] expectedMetadata = {-0.33333333333333326,0.19999999999999996,6,5,1.388888888888889,-1.24,0,0,0,1};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}

//test with scaleX and scaleY
newRaster = RasterEditors.resample(raster, refRaster, true, null);
res = RasterOutputs.asMatrix(newRaster);
expectedRes = "| 1.0 1.0 2.0 3.0 3.0 5.0 5.0|\n" +
"| 1.0 1.0 2.0 3.0 3.0 5.0 5.0|\n" +
"| 4.0 4.0 5.0 6.0 6.0 9.0 9.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0 10.0|\n" +
"| 7.0 7.0 8.0 9.0 9.0 10.0 10.0|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
metadata = RasterAccessors.metadata(newRaster);
expectedMetadata = new double[]{-0.20000000298023224, 0.4000000059604645, 7.0, 5.0, 1.2, -1.4, 0.0, 0.0, 0.0, 1.0};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}
}

@Test
public void testResampleDiffAlgorithms() throws FactoryException, TransformException {
/*
Even though we cannot match interpolation with that of PostGIS for other algorithms, this is a sanity test case to detect potentially invalid changes to the function
*/
double[] values = {1, 2, 3, 4, 5, 6, 7, 8, 9};
GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "d", 3, 3, 0, 0, 2, -2, 0, 0, 0);
raster = MapAlgebra.addBandFromArray(raster, values, 1, null);

//test bilinear
GridCoverage2D newRaster = RasterEditors.resample(raster, 5, 5, 0, 0, false, "bilinear");
String res = RasterOutputs.asMatrix(newRaster);
String expectedRes = "| NaN NaN NaN NaN NaN|\n" +
"| NaN 2.600000 3.200000 3.800000 4.200000|\n" +
"| NaN 4.400000 5.000000 5.600000 6.000000|\n" +
"| NaN 6.200000 6.800000 7.400000 7.800000|\n" +
"| NaN 7.400000 8.000000 8.600000 9.000000|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
double[] metadata = RasterAccessors.metadata(newRaster);
double[] expectedMetadata = {0, 0, 5, 5, 1.2, -1.2, 0, 0, 0, 1};
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}



//test bicubic
newRaster = RasterEditors.resample(raster, 5, 5, 0, 0, false, "bicubic");
res = RasterOutputs.asMatrix(newRaster);
expectedRes = "| NaN NaN NaN NaN NaN|\n" +
"| NaN 2.305379 2.979034 3.648548 4.042909|\n" +
"| NaN 4.326345 5.000000 5.669513 6.063874|\n" +
"| NaN 6.334885 7.008540 7.678053 8.072415|\n" +
"| NaN 7.517968 8.191623 8.861137 9.255498|\n";
//verify correct interpolation
assertEquals(expectedRes, res);
metadata = RasterAccessors.metadata(newRaster);
//verify correct raster geometry
for (int i = 0; i < metadata.length; i++) {
assertEquals(expectedMetadata[i], metadata[i], 1e-6);
}
}

@Test
public void testResampleRefRasterDiffSRID() throws FactoryException {
GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "d", 3, 3, 0, 0, 2, -2, 0, 0, 0);
GridCoverage2D refRaster = RasterConstructors.makeEmptyRaster(2, "d", 5, 5, 1, -1, 1.2, -1.2, 0, 0, 4326);
assertThrows("Provided input raster and reference raster have different SRIDs", IllegalArgumentException.class, () -> RasterEditors.resample(raster, refRaster, false, null));
}

}
Loading

0 comments on commit b135e8e

Please sign in to comment.