You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ji...@apache.org on 2023/08/07 17:59:47 UTC
[sedona] branch master updated: [SEDONA-344] Add RS_RasterToWorldCoordX and RS_RasterToWorldCoordY (#947)
This is an automated email from the ASF dual-hosted git repository.
jiayu pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sedona.git
The following commit(s) were added to refs/heads/master by this push:
new af67bd77 [SEDONA-344] Add RS_RasterToWorldCoordX and RS_RasterToWorldCoordY (#947)
af67bd77 is described below
commit af67bd77ab5c0b5b6b07cc517b8aaa2640655045
Author: Nilesh Gajwani <ni...@gmail.com>
AuthorDate: Mon Aug 7 13:59:40 2023 -0400
[SEDONA-344] Add RS_RasterToWorldCoordX and RS_RasterToWorldCoordY (#947)
---
.../sedona/common/raster/PixelFunctions.java | 36 +++++------------
.../sedona/common/raster/RasterAccessors.java | 9 +++++
.../apache/sedona/common/utils/RasterUtils.java | 17 ++++++++
.../apache/sedona/common/raster/FunctionsTest.java | 2 +-
.../sedona/common/raster/RasterAccessorsTest.java | 47 ++++++++++++++++++++++
docs/api/sql/Raster-operators.md | 36 +++++++++++++++++
.../scala/org/apache/sedona/sql/UDF/Catalog.scala | 2 +
.../expressions/raster/RasterAccessors.scala | 12 ++++++
.../org/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 5ac5561e..39bffe01 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 class PixelFunctions
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) {
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 3ab1d27f..ea830410 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.crs.DefaultEngineeringCRS;
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 class RasterAccessors
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 b60b3cd5..d42a4973 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.GridCoverageFactory;
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 class RasterUtils {
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 9a1f2033..60a8c00b 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 class FunctionsTest extends RasterTestBase {
}
@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 060c4838..4dcd60dd 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 @@ package org.apache.sedona.common.raster;
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 class RasterAccessorsTest extends RasterTestBase
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 0e8617a9..f669d041 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 f10883ea..c02ebf1b 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 e1e9bd11..cc121e0f 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 3b94916f..86324230 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))