From 5bd0172a8865b162d93320d1229ec59f57a79088 Mon Sep 17 00:00:00 2001 From: Furqaan Khan <46216254+furqaankhan@users.noreply.github.com> Date: Wed, 1 Nov 2023 01:22:29 +0530 Subject: [PATCH] [SEDONA-411] Add RS_Rotation (#1061) --- .../sedona/common/raster/RasterAccessors.java | 64 +++++++++++++++++++ .../common/raster/RasterAccessorsTest.java | 19 +++++- docs/api/sql/Raster-operators.md | 54 ++++++++++++++++ .../org/apache/sedona/sql/UDF/Catalog.scala | 2 + .../expressions/raster/RasterAccessors.scala | 12 ++++ .../apache/sedona/sql/rasteralgebraTest.scala | 14 ++++ 6 files changed, 163 insertions(+), 2 deletions(-) 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 5d1d82cc2b..20062d4254 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 @@ -32,6 +32,7 @@ import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.opengis.referencing.operation.TransformException; +import java.awt.image.RenderedImage; import java.util.Arrays; import java.util.Set; @@ -127,6 +128,69 @@ public static String getGeoReference(GridCoverage2D raster, String format) { } } + + public static double[] getGeoTransform(GridCoverage2D raster) throws FactoryException { + // size of a pixel along the transformed i axis + double magnitudeI; + + // size of a pixel along the transformed j axis + double magnitudeJ; + + // angle by which the raster is rotated (Radians positive clockwise) + double thetaI; + + // angle from transformed i axis to transformed j axis (Radians positive counter-clockwise) + double thetaIJ; + + double[] metadata = metadata(raster); + + // x ordinate of the upper-left corner of the upper-left pixel + double offsetX = metadata[0]; + + // y ordinate of the upper-left corner of the upper-left pixel + double offsetY = metadata[1]; + + double scaleX = metadata[4]; + double scaleY = metadata[5]; + double skewX = metadata[6]; + double skewY = metadata[7]; + + // pixel size in i direction - west-east + magnitudeI = Math.sqrt(scaleX * scaleX + skewY * skewY); + + // pixel size in j direction - north-south + magnitudeJ = Math.sqrt(scaleY * scaleY + skewX * skewX); + + // Rotation + thetaI = Math.acos(scaleX / magnitudeI); + double thetaTest = Math.acos(skewY / magnitudeI); + if (thetaTest < Math.PI / 2) { + thetaI = -thetaI; + } + + // Angular separation + thetaIJ = Math.acos((((scaleX * skewX) + (skewY * scaleY)) / (magnitudeI * magnitudeJ))); + thetaTest = Math.acos(((-skewY * skewX) + (scaleX * scaleY)) / (magnitudeI * magnitudeJ)); + if (thetaTest > Math.PI / 2) { + thetaIJ = -thetaIJ; + } + + double[] result = new double[6]; + result[0] = magnitudeI; + result[1] = magnitudeJ; + result[2] = thetaI; + result[3] = thetaIJ; + result[4] = offsetX; + result[5] = offsetY; + + return result; + } + + public static double getRotation(GridCoverage2D raster) throws FactoryException { + double rotation = getGeoTransform(raster)[2]; + return rotation; + } + public static Geometry getGridCoord(GridCoverage2D raster, double x, double y) throws TransformException { int[] coords = RasterUtils.getGridCoordinatesFromWorld(raster, x, y); coords = Arrays.stream(coords).map(number -> number + 1).toArray(); 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 541b7fa81c..2193225f0d 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 @@ -29,8 +29,7 @@ import java.io.IOException; -import static org.junit.Assert.assertEquals; -import static org.junit.Assert.assertThrows; +import static org.junit.Assert.*; public class RasterAccessorsTest extends RasterTestBase { @@ -60,6 +59,22 @@ public void testSrid() throws FactoryException { assertEquals(4326, RasterAccessors.srid(multiBandRaster)); } + @Test + public void testRotation() throws IOException, FactoryException { + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(2, 10, 15, 1, 2, 1, -2, 10, 10, 0); + double actual = RasterAccessors.getRotation(emptyRaster); + double expected = -1.4711276743037347; + assertEquals(expected, actual, 1e-9); + } + + @Test + public void testGeoTransform() throws FactoryException { + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 10, 15, 1, 2, 1, -1, 10, 10, 0); + double[] actual = RasterAccessors.getGeoTransform(emptyRaster); + double[] expected = new double[] {10.04987562112089, 10.04987562112089, -1.4711276743037347, -1.5707963267948966, 1.0, 2.0}; + assertArrayEquals(expected, actual, 1e-9); + } + @Test public void testUpperLeftX() throws FactoryException { GridCoverage2D gridCoverage2D = RasterConstructors.makeEmptyRaster(1, 3, 4, 1,2, 5); diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md index 4bda003051..3bafeb4480 100644 --- a/docs/api/sql/Raster-operators.md +++ b/docs/api/sql/Raster-operators.md @@ -292,6 +292,38 @@ SELECT RS_GeoReferrence(ST_MakeEmptyRaster(1, 3, 4, 100.0, 200.0,2.0, -3.0, 0.1, 198.500000 ``` +### RS_GeoTransform + +Introduction: Returns an array of parameters that represent the GeoTranformation of the raster. The array contains the following values: + +- 0: pixel width along west-east axis (x axis) +- 1: pixel height along north-south axis (y axis) +- 2: Rotation of the raster +- 3: Angular separation between x axis and y axis +- 4: X ordinate of upper-left coordinate +- 5: Y ordinate of upper-left coordinate + +!!!note + Refer to [this image](https://www.researchgate.net/figure/Relation-between-the-cartesian-axes-x-y-and-i-j-axes-of-the-pixels_fig3_313860913) for a clear understanding between i & j axis and x & y axis. + +Format: `RS_GeoTransform(raster: Raster)` + +Since: `v1.5.1` + +Spark SQL Example: + +```sql +SELECT RS_GeoTransform( + RS_MakeEmptyRaster(2, 10, 15, 1, 2, 1, -2, 1, 2, 0) + ) +``` + +Output: + +``` +[2.23606797749979, 2.23606797749979, -1.1071487177940904, -2.214297435588181, 1.0, 2.0] +``` + ### RS_Height Introduction: Returns the height of the raster. @@ -352,6 +384,28 @@ Output: 54 ``` +### RS_Rotation + +Introduction: Returns the uniform rotation of the raster in radian. + +Format: `RS_Rotation(raster: Raster)` + +Since: `v1.5.1` + +Spark SQL Example: + +```sql +SELECT RS_Rotation( + RS_MakeEmptyRaster(2, 10, 15, 1, 2, 1, -2, 1, 2, 0) + ) +``` + +Output: + +``` +-1.1071487177940904 +``` + ### RS_ScaleX Introduction: Returns the pixel width of the raster in CRS units. diff --git a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala index 98e258584d..92238b7427 100644 --- a/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala +++ b/spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala @@ -226,6 +226,8 @@ object Catalog { function[RS_SkewX](), function[RS_SkewY](), function[RS_GeoReference](), + function[RS_Rotation](), + function[RS_GeoTransform](), function[RS_PixelAsPoint](), function[RS_PixelAsPolygon](), function[RS_PixelAsCentroid](), diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala index 77aa82eccc..e6f13eb3fe 100644 --- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterAccessors.scala @@ -86,6 +86,18 @@ case class RS_GeoReference(inputExpressions: Seq[Expression]) extends InferredEx } } +case class RS_Rotation(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getRotation _) { + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } +} + +case class RS_GeoTransform(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getGeoTransform _) { + protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = { + copy(inputExpressions = newChildren) + } +} + case class RS_SkewX(inputExpressions: Seq[Expression]) extends InferredExpression(RasterAccessors.getSkewX _) { protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]): Expression = { copy(inputExpressions = newChildren) diff --git a/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala b/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala index e533b8fda7..ba4767c4c2 100644 --- a/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala +++ b/spark/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala @@ -609,6 +609,20 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen assertEquals(expected, result) } + it("Passed RS_Rotation") { + val df = sparkSession.sql("SELECT RS_MakeEmptyRaster(2, 10, 15, 1, 2, 1, -2, 1, 2, 0) as raster") + val actual = df.selectExpr("RS_Rotation(raster)").first().get(0) + val expected = -1.1071487177940904 + assertEquals(expected, actual) + } + + it("Passed RS_GeoTransform") { + val df = sparkSession.sql("SELECT RS_MakeEmptyRaster(2, 10, 15, 1, 2, 1, -2, 1, 2, 0) as raster") + val actual = df.selectExpr("RS_GeoTransform(raster)").first().getSeq(0) + val expected = mutable.ArraySeq(2.23606797749979, 2.23606797749979, -1.1071487177940904, -2.214297435588181, 1.0, 2.0) + assertTrue(expected.equals(actual)) + } + it("Passed RS_SkewX") { val df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff") val result = df.selectExpr("RS_SkewX(RS_FromGeoTiff(content))").first().getDouble(0)