apache · jiayuasu · Oct 31, 2023 · Oct 25, 2023 · Oct 25, 2023 · Oct 25, 2023
@@ -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,70 @@ 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;
+
+        RenderedImage renderedImage = raster.getRenderedImage();
+
+        // x ordinate of the upper-left corner of the upper-left pixel
+        double offsetX = renderedImage.getTileGridXOffset();
+
+        // y ordinate of the upper-left corner of the upper-left pixel
+        double offsetY = renderedImage.getTileGridYOffset();
+
+        double[] metadata = metadata(raster);
+        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();

@@ -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, 0.0, 0.0};
+        assertArrayEquals(expected, actual, 1e-9);
+    }
+
     @Test
     public void testUpperLeftX() throws FactoryException {
         GridCoverage2D gridCoverage2D = RasterConstructors.makeEmptyRaster(1, 3, 4, 1,2, 5);

@@ -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: Grid offset X for upper-left X coordinate
+- 5: Grid offset Y for upper-left Y 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, 0.0, 0.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.

@@ -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](),

@@ -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)

@@ -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, 0.0, 0.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)