[SEDONA-344] Add RS_RasterToWorldCoordX and RS_RasterToWorldCoordY (#947

)

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

@@ -18,16 +18,13 @@
  */
 package org.apache.sedona.common.raster;
 
-import org.geotools.coverage.grid.GridCoordinates2D;
+import org.apache.sedona.common.utils.RasterUtils;
 import org.geotools.coverage.grid.GridCoverage2D;
-import org.geotools.coverage.grid.GridEnvelope2D;
-import org.geotools.coverage.grid.GridGeometry2D;
 import org.geotools.geometry.DirectPosition2D;
-import org.geotools.referencing.CRS;
 import org.locationtech.jts.geom.*;
 import org.opengis.coverage.PointOutsideCoverageException;
 import org.opengis.geometry.DirectPosition;
-import org.opengis.metadata.spatial.PixelOrientation;
+import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.operation.TransformException;
 
 import java.awt.geom.Point2D;
@@ -46,29 +43,16 @@ public static Double value(GridCoverage2D rasterGeom, Geometry geometry, int ban
         return values(rasterGeom, Collections.singletonList(geometry), band).get(0);
     }
 
-    public static Geometry getPixelAsPoint(GridCoverage2D raster, int colX, int rowY) throws TransformException {
-        GridGeometry2D gridGeometry2D = raster.getGridGeometry();
-        GridEnvelope2D gridEnvelope2D = gridGeometry2D.getGridRange2D();
-        GridCoordinates2D gridCoordinates2D = new GridCoordinates2D(colX - 1, rowY - 1);
-        int srid = 0;
-        String srs = CRS.toSRS(raster.getCoordinateReferenceSystem2D(), true);
-        if (!"Generic cartesian 2D".equalsIgnoreCase(srs)) {
-           srid = Integer.parseInt(srs);
-        }
-        if (gridEnvelope2D.contains(gridCoordinates2D)) {
-            Point2D point2D = gridGeometry2D.getGridToCRS2D(PixelOrientation.UPPER_LEFT).transform(gridCoordinates2D, null);
-            DirectPosition2D directPosition2D = new DirectPosition2D(gridGeometry2D.getCoordinateReferenceSystem2D(), point2D.getX(), point2D.getY());
-            Coordinate pointCoord = new Coordinate(directPosition2D.getX(), directPosition2D.getY());
-            if (srid != 0) {
-                GeometryFactory factory = new GeometryFactory(new PrecisionModel(), srid);
-                return factory.createPoint(pointCoord);
-            }
-            return GEOMETRY_FACTORY.createPoint(pointCoord);
-        }else {
-            throw new IndexOutOfBoundsException(String.format("Specified pixel coordinates (%d, %d) do not lie in the raster", colX, rowY));
+    public static Geometry getPixelAsPoint(GridCoverage2D raster, int colX, int rowY) throws TransformException, FactoryException {
+        int srid = RasterAccessors.srid(raster);
+        Point2D point2D = RasterUtils.getCornerCoordinatesWithRangeCheck(raster, colX, rowY);
+        Coordinate pointCoord = new Coordinate(point2D.getX(), point2D.getY());
+        if (srid != 0) {
+            GeometryFactory factory = new GeometryFactory(new PrecisionModel(), srid);
+            return factory.createPoint(pointCoord);
         }
+        return GEOMETRY_FACTORY.createPoint(pointCoord);
     }
-
     public static List<Double> values(GridCoverage2D rasterGeom, List<Geometry> geometries, int band) throws TransformException {
         int numBands = rasterGeom.getNumSampleDimensions();
         if (band < 1 || band > numBands) {

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

@@ -26,6 +26,7 @@
 import org.geotools.referencing.operation.transform.AffineTransform2D;
 import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.TransformException;
 
 import java.util.Optional;
 
@@ -74,6 +75,14 @@ public static double getScaleY(GridCoverage2D raster) {
         return RasterUtils.getGDALAffineTransform(raster).getScaleY();
     }
 
+    public static double getWorldCoordX(GridCoverage2D raster, int colX, int rowY) throws TransformException {
+        return RasterUtils.getCornerCoordinates(raster, colX, rowY).getX();
+    }
+
+    public static double getWorldCoordY(GridCoverage2D raster, int colX, int rowY) throws TransformException {
+        return RasterUtils.getCornerCoordinates(raster, colX, rowY).getY();
+    }
+
     /**
      * Returns the metadata of a raster as an array of doubles.
      * @param raster the raster

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

@@ -27,6 +27,8 @@
 import org.geotools.coverage.grid.GridGeometry2D;
 import org.geotools.geometry.DirectPosition2D;
 import org.geotools.referencing.operation.transform.AffineTransform2D;
+import org.hsqldb.index.Index;
+import org.opengis.coverage.grid.GridCoordinates;
 import org.opengis.metadata.spatial.PixelOrientation;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
@@ -96,4 +98,19 @@ public static AffineTransform2D getAffineTransform(GridCoverage2D raster, PixelO
     public static Point2D getCornerCoordinates(GridCoverage2D raster, int colX, int rowY) throws TransformException {
         return raster.getGridGeometry().getGridToCRS2D(PixelOrientation.UPPER_LEFT).transform(new GridCoordinates2D(colX - 1, rowY - 1), null);
     }
+
+    /***
+     * Returns the world coordinates of the given grid coordinate. The expected grid coordinates are 1 indexed. The function also enforces a range check to make sure given grid coordinates are actually inside the grid.
+     * @param raster
+     * @param colX
+     * @param rowY
+     * @return
+     * @throws IndexOutOfBoundsException
+     * @throws TransformException
+     */
+    public static Point2D getCornerCoordinatesWithRangeCheck(GridCoverage2D raster, int colX, int rowY) throws IndexOutOfBoundsException, TransformException {
+        GridCoordinates2D gridCoordinates2D = new GridCoordinates2D(colX - 1, rowY - 1);
+        if (!(raster.getGridGeometry().getGridRange2D().contains(gridCoordinates2D))) throw new IndexOutOfBoundsException(String.format("Specified pixel coordinates (%d, %d) do not lie in the raster", colX, rowY));
+        return raster.getGridGeometry().getGridToCRS2D(PixelOrientation.UPPER_LEFT).transform(gridCoordinates2D, null);
+    }
 }

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

@@ -116,7 +116,7 @@ public void testPixelAsPointOutOfBounds() throws FactoryException {
     }
 
     @Test
-    public void testPixelAsPointFromRasterFile() throws IOException, TransformException {
+    public void testPixelAsPointFromRasterFile() throws IOException, TransformException, FactoryException {
         GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster/test1.tiff");
         Geometry actualPoint = PixelFunctions.getPixelAsPoint(raster, 1, 1);
         Coordinate coordinate = actualPoint.getCoordinate();

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

@@ -21,6 +21,7 @@
 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;
 
@@ -88,6 +89,52 @@ public void testScaleY() throws UnsupportedOperationException, FactoryException
         GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(2, 10, 15, 0, 0, 1, -2, 0, 0, 0);
         assertEquals(-2, RasterAccessors.getScaleY(emptyRaster), 1e-9);
     }
+
+    @Test
+    public void testWorldCoordX() throws FactoryException, TransformException {
+        int colX = 1, rowY = 1;
+        GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326);
+        double actualX = RasterAccessors.getWorldCoordX(emptyRaster, colX, rowY);
+        double expectedX = -123;
+        assertEquals(expectedX,actualX, 0.1d);
+        colX = 2;
+        expectedX = -118;
+        actualX = RasterAccessors.getWorldCoordX(emptyRaster, colX, rowY);
+        assertEquals(expectedX, actualX, 0.1d);
+    }
+
+    @Test
+    public void testWorldCoordXOutOfBounds() throws FactoryException, TransformException {
+        int colX = 6;
+        int rowY = 5;
+        GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326);
+        double actualX = RasterAccessors.getWorldCoordX(emptyRaster, colX, rowY);
+        double expectedX = -98;
+        assertEquals(expectedX, actualX, 0.1d);
+    }
+
+    @Test
+    public void testWorldCoordY() throws FactoryException, TransformException {
+        int colX = 1, rowY = 1;
+        GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326);
+        double actualY = RasterAccessors.getWorldCoordY(emptyRaster, colX, rowY);
+        double expectedY = 54;
+        assertEquals(expectedY,actualY, 0.1d);
+        rowY = 2;
+        expectedY = 44;
+        actualY = RasterAccessors.getWorldCoordY(emptyRaster, colX, rowY);
+        assertEquals(expectedY, actualY, 0.1d);
+    }
+
+    @Test
+    public void testWorldCoordYOutOfBounds() throws FactoryException, TransformException{
+        int colX = 4;
+        int rowY = 11;
+        GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326);
+        double actualY = RasterAccessors.getWorldCoordY(emptyRaster, colX, rowY);
+        double expectedY = -46;
+        assertEquals(expectedY, actualY, 0.1d);
+    }
 
     @Test
     public void testMetaData()

docs/api/sql/Raster-operators.md

@@ -90,6 +90,42 @@ Output:
 512
 ```
 
+### RS_RasterToWorldCoordX
+
+Introduction: Returns the upper left X coordinate of the given row and column of the given raster geometric units of the geo-referenced raster. If any out of bounds values are given, the X coordinate of the assumed point considering existing raster pixel size and skew values will be returned.
+
+Format: `RS_RasterToWorldCoordX(raster: Raster, colX: int, rowY: int)`
+
+Since: `1.5.0`
+
+Spark SQL example:
+```sql
+SELECT RS_RasterToWorldCoordX(ST_MakeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326), 1, 1) from rasters
+```
+
+Output:
+```
+-123
+```
+
+### RS_RasterToWorldCoordY
+
+Introduction: Returns the upper left Y coordinate of the given row and column of the given raster geometric units of the geo-referenced raster. If any out of bounds values are given, the Y coordinate of the assumed point considering existing raster pixel size and skew values will be returned.
+
+Format: `RS_RasterToWorldCoordY(raster: Raster, colX: int, rowY: int)`
+
+Since: `1.5.0`
+
+Spark SQL example:
+```sql
+SELECT RS_RasterToWorldCoordY(ST_MakeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, 4326), 1, 1) from rasters
+```
+
+Output:
+```
+54
+```
+
 ### RS_ScaleX
 
 Introduction: Returns the pixel width of the raster in CRS units.

sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala

@@ -215,6 +215,8 @@ object Catalog {
     function[RS_ScaleY](),
     function[RS_PixelAsPoint](),
     function[RS_ConvexHull](),
+    function[RS_RasterToWorldCoordX](),
+    function[RS_RasterToWorldCoordY](),
     function[RS_Within](),
     function[RS_Contains]()
   )

sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala

@@ -78,3 +78,15 @@ case class RS_ScaleY(inputExpressions: Seq[Expression]) extends InferredExpressi
   }
 }
 
+case class RS_RasterToWorldCoordX(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getWorldCoordX _) {
+  protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+    copy(inputExpressions = newChildren)
+  }
+}
+
+case class RS_RasterToWorldCoordY(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getWorldCoordY _) {
+  protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+    copy(inputExpressions = newChildren)
+  }
+}
+

sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala

@@ -559,6 +559,20 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen
       assertEquals(expectedY, actualCoordinates.y, 1e-5)
     }
 
+    it("Passed RS_RasterToWorldCoordX with raster") {
+      var df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff")
+      df = df.selectExpr("RS_FromGeoTiff(content) as raster")
+      val result = df.selectExpr("RS_RasterToWorldCoordX(raster, 1, 1)").first().getDouble(0)
+      assertEquals(-13095817.809482181, result, 0.5d)
+    }
+
+    it("Passed RS_RasterToWorldCoordY with raster") {
+      var df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff")
+      df = df.selectExpr("RS_FromGeoTiff(content) as raster")
+      val result = df.selectExpr("RS_RasterToWorldCoordY(raster, 1, 1)").first().getDouble(0)
+      assertEquals(4021262.7487925636, result, 0.5d)
+    }
+
     it("Passed RS_Contains") {
       assert(sparkSession.sql("SELECT RS_Contains(RS_MakeEmptyRaster(1, 20, 20, 2, 22, 1), ST_GeomFromWKT('POLYGON ((5 5, 5 10, 10 10, 10 5, 5 5))'))").first().getBoolean(0))
       assert(!sparkSession.sql("SELECT RS_Contains(RS_MakeEmptyRaster(1, 20, 20, 2, 22, 1), ST_GeomFromWKT('POLYGON ((2 2, 2 25, 20 25, 20 2, 2 2))'))").first().getBoolean(0))