[SEDONA-404] Add RS_Resample

common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java

@@ -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
 {
@@ -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)));
+    }
+
 }

common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java

@@ -168,7 +168,7 @@ private static String createPaddedMatrixStringFromDouble(double[] values, int wi
         int maxDecimalPrecision = 0;
         for (double value : values) {
             String[] splitByDecimal = String.valueOf(value).split("\\.");
-            int preDecimal = splitByDecimal[0].length(), postDecimal = Math.min(decimalPrecision, splitByDecimal[1].length());
+            int preDecimal = splitByDecimal[0].length(), postDecimal = Math.min(decimalPrecision, splitByDecimal.length > 1 ? splitByDecimal[1].length() : 0);
             maxDecimalPrecision = Math.max(maxDecimalPrecision, postDecimal);
             int currWidth = preDecimal + 1; //add 1 for space occupied for decimal point
             maxPreDecimal = Math.max(maxPreDecimal, currWidth);

common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java

@@ -285,7 +285,7 @@ public static int[] getGridCoordinatesFromWorld(GridCoverage2D raster, double lo
         DirectPosition2D directPosition2D = new DirectPosition2D(raster.getCoordinateReferenceSystem2D(), longitude, latitude);
         DirectPosition worldCoord = raster.getGridGeometry().getCRSToGrid2D(PixelOrientation.UPPER_LEFT).transform((DirectPosition) directPosition2D, null);
         double[] coords = worldCoord.getCoordinate();
-        int[] gridCoords = new int[] {(int) coords[0], (int) coords[1]};
+        int[] gridCoords = new int[] {(int) Math.floor(coords[0]), (int) Math.floor(coords[1])};
         return gridCoords;
     }