From af67bd77ab5c0b5b6b07cc517b8aaa2640655045 Mon Sep 17 00:00:00 2001 From: Nilesh Gajwani Date: Mon, 7 Aug 2023 13:59:40 -0400 Subject: [PATCH] [SEDONA-344] Add RS_RasterToWorldCoordX and RS_RasterToWorldCoordY (#947) --- .../sedona/common/raster/PixelFunctions.java | 36 ++++---------- .../sedona/common/raster/RasterAccessors.java | 9 ++++ .../sedona/common/utils/RasterUtils.java | 17 +++++++ .../sedona/common/raster/FunctionsTest.java | 2 +- .../common/raster/RasterAccessorsTest.java | 47 +++++++++++++++++++ docs/api/sql/Raster-operators.md | 36 ++++++++++++++ .../org/apache/sedona/sql/UDF/Catalog.scala | 2 + .../expressions/raster/RasterAccessors.scala | 12 +++++ .../apache/sedona/sql/rasteralgebraTest.scala | 14 ++++++ 9 files changed, 148 insertions(+), 27 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java b/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java index 5ac5561ebf..39bffe010e 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/PixelFunctions.java +++ b/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 values(GridCoverage2D rasterGeom, List geometries, int band) throws TransformException { int numBands = rasterGeom.getNumSampleDimensions(); if (band < 1 || band > numBands) { diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java index 3ab1d27fbb..ea83041067 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterAccessors.java +++ b/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 diff --git a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java b/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java index b60b3cd554..d42a497310 100644 --- a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java +++ b/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); + } } diff --git a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java b/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java index 9a1f20331f..60a8c00b5c 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/FunctionsTest.java +++ b/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(); diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java index 060c48389d..4dcd60dde8 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterAccessorsTest.java +++ b/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() diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md index 0e8617a995..f669d041cd 100644 --- a/docs/api/sql/Raster-operators.md +++ b/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. diff --git a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala index f10883eaaf..c02ebf1b16 100644 --- a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala +++ b/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]() ) diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala index e1e9bd11fa..cc121e0f9a 100644 --- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala +++ b/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) + } +} + diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala b/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala index 3b94916fe2..86324230aa 100644 --- a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala +++ b/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))