[SEDONA-437] Add implicit CRS transformation

R/tests/testthat/test-data-interface-raster.R

@@ -165,7 +165,7 @@ test_that("Passed RS_Value with raster", {
     binary_sdf %>% 
     mutate(
       raster = RS_FromGeoTiff(content),
-      val = RS_Value(raster, ST_Point(-13077301.685, 4002565.802))
+      val = RS_Value(raster, ST_SetSRID(ST_Point(-13077301.685, 4002565.802), RS_SRID(raster)))
     ) %>% 
     select(val) %>% 
     collect()
@@ -191,7 +191,7 @@ test_that("Passed RS_Values with raster", {
     binary_sdf %>% 
     mutate(
       raster = RS_FromGeoTiff(content),
-      val = RS_Values(raster, array(ST_Point(-13077301.685, 4002565.802), NULL))
+      val = RS_Values(raster, array(ST_SetSRID(ST_Point(-13077301.685, 4002565.802), RS_SRID(raster)), NULL))
     ) %>% 
     select(val) %>%
     collect()

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

@@ -18,6 +18,8 @@
  */
 package org.apache.sedona.common.raster;
 
+import org.apache.commons.lang3.tuple.Pair;
+import org.apache.sedona.common.Functions;
 import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.grid.GridCoverage2D;
 import org.geotools.geometry.DirectPosition2D;
@@ -108,6 +110,15 @@ public static GridCoverage2D setValues(GridCoverage2D raster, int band, int colX
     public static GridCoverage2D setValues(GridCoverage2D raster, int band, Geometry geom, double value, boolean keepNoData) throws FactoryException, TransformException {
         RasterUtils.ensureBand(raster, band);
 
+        Pair<GridCoverage2D, Geometry> pair = RasterUtils.setDefaultCRSAndTransform(raster, geom);
+        raster = pair.getLeft();
+        geom = pair.getRight();
+
+        // checking if the raster contains the geometry
+        if (!RasterPredicates.rsIntersects(raster, geom)) {
+            throw new IllegalArgumentException("The provided geometry is not intersecting the raster. Please provide a geometry that is in the raster's extent.");
+        }
+
         String bandDataType = RasterBandAccessors.getBandType(raster, band);
 
         GridCoverage2D rasterizedGeom;
@@ -140,7 +151,6 @@ public static GridCoverage2D setValues(GridCoverage2D raster, int band, Geometry
         }
         // Converting geometry to raster and then iterating through them
         else {
-            // Starting pixel location on the given raster
             int[] pixelLocation = RasterUtils.getGridCoordinatesFromWorld(raster, colX, rowY);
             int x = pixelLocation[0], y = pixelLocation[1];
 

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

@@ -18,6 +18,8 @@
  */
 package org.apache.sedona.common.raster;
 
+import org.apache.commons.lang3.tuple.Pair;
+import org.apache.sedona.common.Functions;
 import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.grid.GridCoordinates2D;
 import org.geotools.coverage.grid.GridCoverage2D;
@@ -32,20 +34,19 @@
 import java.awt.geom.Point2D;
 import java.awt.image.Raster;
 import java.util.ArrayList;
+import java.util.Arrays;
 import java.util.Collections;
 import java.util.List;
 
 public class PixelFunctions
 {
     private static GeometryFactory GEOMETRY_FACTORY = new GeometryFactory();
-    public static Double value(GridCoverage2D rasterGeom, Geometry geometry, int band) throws TransformException
-    {
-        return values(rasterGeom, Collections.singletonList(geometry), band).get(0);
+    public static Double value(GridCoverage2D rasterGeom, Geometry geometry, int band) throws TransformException, FactoryException {
+        return values(rasterGeom, Arrays.asList(geometry), band).get(0);
     }
 
-    public static Double value(GridCoverage2D rasterGeom, Geometry geometry) throws TransformException
-    {
-        return values(rasterGeom, Collections.singletonList(geometry), 1).get(0);
+    public static Double value(GridCoverage2D rasterGeom, Geometry geometry) throws TransformException, FactoryException {
+        return values(rasterGeom, Arrays.asList(geometry), 1).get(0);
     }
 
     public static Double value(GridCoverage2D rasterGeom, int colX, int rowY, int band) throws TransformException
@@ -247,10 +248,26 @@ public static List<Double> values(GridCoverage2D rasterGeom, int[] xCoordinates,
         return result;
     }
 
-    public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries, int band) throws TransformException {
+    public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries, int band) throws TransformException, FactoryException {
         RasterUtils.ensureBand(rasterGeom, band); // Check for invalid band index
-        int numBands = rasterGeom.getNumSampleDimensions();
 
+        for (int i = 0; i < geometries.size(); i++) {
+
+            if (geometries.get(i) == null) {
+                continue;
+            }
+
+
+            Pair<GridCoverage2D, Geometry> pair = RasterUtils.setDefaultCRSAndTransform(rasterGeom, geometries.get(i));
+
+            geometries.set(i, pair.getRight());
+
+            if (i == 0) {
+                rasterGeom = pair.getLeft();
+            }
+        }
+
+        int numBands = rasterGeom.getNumSampleDimensions();
         double noDataValue = RasterUtils.getNoDataValue(rasterGeom.getSampleDimension(band - 1));
         double[] pixelBuffer = new double[numBands];
 
@@ -278,7 +295,7 @@ public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geom
         return result;
     }
 
-    public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries) throws TransformException {
+    public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries) throws TransformException, FactoryException {
         return values(rasterGeom, geometries, 1);
     }
 

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

@@ -19,6 +19,8 @@
 package org.apache.sedona.common.raster;
 
 import org.apache.commons.lang3.ArrayUtils;
+import org.apache.commons.lang3.tuple.Pair;
+import org.apache.sedona.common.Functions;
 import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.GridSampleDimension;
 import org.geotools.coverage.grid.GridCoverage2D;
@@ -189,6 +191,10 @@ public static GridCoverage2D clip(GridCoverage2D raster, int band, Geometry geom
         RasterUtils.ensureBand(raster, band);
         GridCoverage2D singleBandRaster = RasterBandAccessors.getBand(raster, new int[]{band});
 
+        Pair<GridCoverage2D, Geometry> pair = RasterUtils.setDefaultCRSAndTransform(singleBandRaster, geometry);
+        singleBandRaster = pair.getLeft();
+        geometry = pair.getRight();
+
         // Crop the raster
         // this will shrink the extent of the raster to the geometry
         Crop cropObject = new Crop();

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

@@ -19,8 +19,11 @@
 package org.apache.sedona.common.utils;
 
 import com.sun.media.imageioimpl.common.BogusColorSpace;
+import org.apache.commons.lang3.tuple.Pair;
+import org.apache.sedona.common.Functions;
 import org.apache.sedona.common.FunctionsGeoTools;
 import org.apache.sedona.common.raster.RasterAccessors;
+import org.apache.sedona.common.raster.RasterEditors;
 import org.geotools.coverage.Category;
 import org.geotools.coverage.CoverageFactoryFinder;
 import org.geotools.coverage.GridSampleDimension;
@@ -436,6 +439,38 @@ public static Geometry convertCRSIfNeeded(Geometry geometry, CoordinateReference
         return geometry;
     }
 
+    /**
+     * If the raster has a CRS, then it transforms the geom to the raster's CRS. 
+     * If any of the inputs, raster or geom doesn't have a CRS, it defaults to 4326.
+     * @param raster
+     * @param geom
+     * @return
+     * @throws FactoryException
+     */
+    public static Pair<GridCoverage2D, Geometry> setDefaultCRSAndTransform(GridCoverage2D raster, Geometry geom) throws FactoryException {
+        int rasterSRID = RasterAccessors.srid(raster);
+        int geomSRID = Functions.getSRID(geom);
+
+        if (rasterSRID == 0) {
+            raster = RasterEditors.setSrid(raster, 4326);
+            rasterSRID = RasterAccessors.srid(raster);
+        }
+
+        if (geomSRID == 0) {
+            geom = Functions.setSRID(geom, 4326);
+            geomSRID = Functions.getSRID(geom);
+        }
+
+
+        if (rasterSRID != geomSRID) {
+            // implicitly converting roi geometry CRS to raster CRS
+            geom = convertCRSIfNeeded(geom, raster.getCoordinateReferenceSystem());
+            // have to set the SRID as RasterUtils.convertCRSIfNeeded doesn't set it even though the geometry is in raster's CRS
+            geom = Functions.setSRID(geom, RasterAccessors.srid(raster));
+        }
+        return Pair.of(raster, geom);
+    }
+
     /**
      * Converts data types to data type codes
      * @param s pixel data type possible values {D, I, B, F, S, US} <br><br>

common/src/test/java/org/apache/sedona/common/FunctionsTest.java

@@ -1167,23 +1167,23 @@ public void nRingsMultiPolygonMixed() throws Exception {
     @Test
     public void testBuffer() {
         Polygon polygon = GEOMETRY_FACTORY.createPolygon(coordArray(50, 50, 50, 150, 150, 150, 150, 50, 50, 50));
-        String actual = Functions.asWKT(Functions.buffer(polygon, 15));
-        String expected = "POLYGON ((50 35, 47.07364516975807 35.288220793951545, 44.25974851452364 36.1418070123307, 41.666446504705966 37.52795581546182, 39.39339828220179 39.39339828220179, 37.527955815461816 41.66644650470597, 36.141807012330695 44.25974851452366, 35.288220793951545 47.07364516975807, 35 50, 35 150, 35.288220793951545 152.92635483024193, 36.1418070123307 155.74025148547634, 37.52795581546182 158.33355349529404, 39.39339828220179 160.6066017177982, 41.66644650470597 162.4720441845382, 44.25974851452365 163.8581929876693, 47.07364516975808 164.71177920604845, 50 165, 150 165, 152.92635483024193 164.71177920604845, 155.74025148547634 163.8581929876693, 158.33355349529404 162.4720441845382, 160.6066017177982 160.6066017177982, 162.4720441845382 158.33355349529404, 163.8581929876693 155.74025148547634, 164.71177920604845 152.92635483024193, 165 150, 165 50, 164.71177920604845 47.07364516975807, 163.8581929876693 44.25974851452365, 162.4720441845382 41.666446504705966, 160.6066017177982 39.39339828220179, 158.33355349529404 37.52795581546182, 155.74025148547634 36.1418070123307, 152.92635483024193 35.288220793951545, 150 35, 50 35))";
+        String actual = Functions.asWKT(Functions.reducePrecision(Functions.buffer(polygon, 15), 4));
+        String expected = "POLYGON ((47.0736 35.2882, 44.2597 36.1418, 41.6664 37.528, 39.3934 39.3934, 37.528 41.6664, 36.1418 44.2597, 35.2882 47.0736, 35 50, 35 150, 35.2882 152.9264, 36.1418 155.7403, 37.528 158.3336, 39.3934 160.6066, 41.6664 162.472, 44.2597 163.8582, 47.0736 164.7118, 50 165, 150 165, 152.9264 164.7118, 155.7403 163.8582, 158.3336 162.472, 160.6066 160.6066, 162.472 158.3336, 163.8582 155.7403, 164.7118 152.9264, 165 150, 165 50, 164.7118 47.0736, 163.8582 44.2597, 162.472 41.6664, 160.6066 39.3934, 158.3336 37.528, 155.7403 36.1418, 152.9264 35.2882, 150 35, 50 35, 47.0736 35.2882))";
         assertEquals(expected, actual);
 
         LineString lineString = GEOMETRY_FACTORY.createLineString(coordArray(0, 0, 50, 70, 100, 100));
-        actual = Functions.asWKT(Functions.buffer(lineString, 10, "side=left"));
-        expected = "POLYGON ((100 100, 50 70, 0 0, -8.137334712067348 5.812381937190963, 41.86266528793265 75.81238193719096, 43.21673095875923 77.34760240582902, 44.855042445724735 78.57492925712545, 94.85504244572473 108.57492925712545, 100 100))";
+        actual = Functions.asWKT(Functions.reducePrecision(Functions.buffer(lineString, 10, "side=left"), 4));
+        expected = "POLYGON ((50 70, 0 0, -8.1373 5.8124, 41.8627 75.8124, 43.2167 77.3476, 44.855 78.5749, 94.855 108.5749, 100 100, 50 70))";
         assertEquals(expected, actual);
 
         lineString = GEOMETRY_FACTORY.createLineString(coordArray(0, 0, 50, 70, 70, -3));
-        actual = Functions.asWKT(Functions.buffer(lineString, 10, "endcap=square"));
-        expected = "POLYGON ((41.86266528793265 75.81238193719096, 43.21555008457904 77.3465120530184, 44.85228625762473 78.57327494173381, 46.70439518001618 79.44134465372912, 48.69438734657371 79.914402432785, 50.73900442057982 79.9726562392556, 52.75270263976913 79.6136688198111, 54.65123184115194 78.8524596785218, 56.355160363552315 77.72087668296376, 57.79319835113832 76.26626359641972, 58.90518041582699 74.54947928466231, 59.64458286836891 72.642351470786, 79.64458286836891 -0.3576485292139977, 82.28693433915491 -10.002231397582907, 62.997768602417096 -15.28693433915491, 45.912786454208465 47.07325050180659, 8.137334712067348 -5.812381937190963, 2.324952774876386 -13.949716649258315, -13.94971664925831 -2.3249527748763885, 41.86266528793265 75.81238193719096))";
+        actual = Functions.asWKT(Functions.reducePrecision(Functions.buffer(lineString, 10, "endcap=square"), 4));
+        expected = "POLYGON ((43.2156 77.3465, 44.8523 78.5733, 46.7044 79.4413, 48.6944 79.9144, 50.739 79.9727, 52.7527 79.6137, 54.6512 78.8525, 56.3552 77.7209, 57.7932 76.2663, 58.9052 74.5495, 59.6446 72.6424, 79.6446 -0.3576, 82.2869 -10.0022, 62.9978 -15.2869, 45.9128 47.0733, 8.1373 -5.8124, 2.325 -13.9497, -13.9497 -2.325, 41.8627 75.8124, 43.2156 77.3465))";
         assertEquals(expected, actual);
 
         Point point = GEOMETRY_FACTORY.createPoint(new Coordinate(100, 90));
-        actual = Functions.asWKT(Functions.buffer(point, 10, "quad_segs=2"));
-        expected = "POLYGON ((110 90, 107.07106781186548 82.92893218813452, 100 80, 92.92893218813452 82.92893218813452, 90 90, 92.92893218813452 97.07106781186548, 100 100, 107.07106781186548 97.07106781186548, 110 90))";
+        actual = Functions.asWKT(Functions.reducePrecision(Functions.buffer(point, 10, "quad_segs=2"), 4));
+        expected = "POLYGON ((107.0711 82.9289, 100 80, 92.9289 82.9289, 90 90, 92.9289 97.0711, 100 100, 107.0711 97.0711, 110 90, 107.0711 82.9289))";
         assertEquals(expected, actual);
     }
 

common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java

@@ -19,6 +19,9 @@
 package org.apache.sedona.common.raster;
 
 import org.apache.sedona.common.Constructors;
+import org.apache.sedona.common.Functions;
+import org.apache.sedona.common.FunctionsGeoTools;
+import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.grid.GridCoverage2D;
 import org.junit.Test;
 import org.locationtech.jts.geom.Geometry;
@@ -72,6 +75,54 @@ public void testSetValuesWithGeomInRaster() throws IOException, ParseException,
         assertEquals(expected, actual, 0d);
     }
 
+    @Test
+    public void testSetValuesWithRasterNoSRID() throws IOException, ParseException, FactoryException, TransformException {
+        GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif");
+        String polygon = "POLYGON ((-77.9148 37.9545, -77.9123 37.8898, -77.9938 37.8878, -77.9964 37.9524, -77.9148 37.9545))";
+        Geometry geom = Constructors.geomFromWKT(polygon, 0);
+
+        GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false);
+
+        Geometry point = Constructors.geomFromWKT("POINT (-77.9146 37.8916)", 0);
+        double actual = PixelFunctions.value(result, point, 1);
+        double expected = 10.0;
+        assertEquals(expected, actual, 0d);
+
+        point = Constructors.geomFromWKT("POINT (-77.9549 37.9357)", 0);
+        actual = PixelFunctions.value(result, point, 1);
+        assertEquals(expected, actual, 0d);
+    }
+
+    @Test
+    public void testSetValuesWithRaster() throws IOException, FactoryException, ParseException, TransformException {
+        GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif");
+        String polygon = "POLYGON ((-8682522.873537656 4572703.890837922, -8673439.664183248 4572993.532747675, -8673155.57366801 4563873.2099182755, -8701890.325907696 4562931.7093397, -8682522.873537656 4572703.890837922))";
+        Geometry geom = Constructors.geomFromWKT(polygon, 3857);
+
+        GridCoverage2D result = PixelFunctionEditors.setValues(raster, 1, geom, 10, false);
+
+        Geometry point = Constructors.geomFromWKT("POINT (243700 4197797)", 26918);
+        double actual = PixelFunctions.value(result, point, 1);
+        double expected = 10.0;
+        assertEquals(expected, actual, 0d);
+
+        point = Constructors.geomFromWKT("POINT (235749.0869 4200557.7397)", 26918);
+        actual = PixelFunctions.value(result, point, 1);
+        assertEquals(expected, actual, 0d);
+
+        point = Constructors.geomFromWKT("POINT (240311 4202806)", 26918);
+        actual = PixelFunctions.value(result, point, 1);
+        assertEquals(expected, actual, 0d);
+
+        point = Constructors.geomFromWKT("POINT (223670 4197650)", 26918);
+        actual = PixelFunctions.value(result, point, 1);
+        assertEquals(expected, actual, 0d);
+
+        point = Constructors.geomFromWKT("POINT (243850 4197640)", 26918);
+        actual = PixelFunctions.value(result, point, 1);
+        assertEquals(expected, actual, 0d);
+    }
+
     @Test
     public void testSetValuesGeomVariant() throws FactoryException, ParseException, TransformException {
         GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 4, 6, 1, -1, 1, -1, 0, 0, 0);