[SEDONA-432] Add RS_PixelAsCentroids

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

@@ -22,13 +22,15 @@
 import org.geotools.coverage.grid.GridCoordinates2D;
 import org.geotools.coverage.grid.GridCoverage2D;
 import org.geotools.geometry.DirectPosition2D;
+import org.geotools.referencing.operation.transform.AffineTransform2D;
 import org.locationtech.jts.geom.*;
 import org.opengis.coverage.PointOutsideCoverageException;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.FactoryException;
 import org.opengis.referencing.operation.TransformException;
 
 import java.awt.geom.Point2D;
+import java.awt.image.Raster;
 import java.util.ArrayList;
 import java.util.Collections;
 import java.util.List;
@@ -80,6 +82,43 @@ public static Geometry getPixelAsCentroid(GridCoverage2D raster, int colX, int r
         return  polygon.getCentroid();
     }
 
+    public static List<PixelRecord> getPixelAsCentroids(GridCoverage2D rasterGeom, int band) throws TransformException, FactoryException {
+        RasterUtils.ensureBand(rasterGeom, band);
+
+        int width = RasterAccessors.getWidth(rasterGeom);
+        int height = RasterAccessors.getHeight(rasterGeom);
+
+        Raster r = RasterUtils.getRaster(rasterGeom.getRenderedImage());
+        double[] pixels = r.getSamples(0, 0, width, height, band - 1, (double[]) null);
+
+        AffineTransform2D gridToCRS = RasterUtils.getGDALAffineTransform(rasterGeom);
+        double cellSizeX = gridToCRS.getScaleX();
+        double cellSizeY = gridToCRS.getScaleY();
+        double shearX = gridToCRS.getShearX();
+        double shearY = gridToCRS.getShearY();
+
+        int srid = RasterAccessors.srid(rasterGeom);
+        GeometryFactory geometryFactory = srid != 0 ? new GeometryFactory(new PrecisionModel(), srid) : GEOMETRY_FACTORY;
+
+        Point2D upperLeft = RasterUtils.getWorldCornerCoordinates(rasterGeom, 1, 1);
+        List<PixelRecord> pixelRecords = new ArrayList<>();
+
+        for (int y = 1; y <= height; y++) {
+            for (int x = 1; x <= width; x++) {
+                double pixelValue = pixels[(y - 1) * width + (x - 1)];
+
+                double worldX = upperLeft.getX() + (x - 0.5) * cellSizeX + (y - 0.5) * shearX;
+                double worldY = upperLeft.getY() + (y - 0.5) * cellSizeY + (x - 0.5) * shearY;
+
+                Coordinate centroidCoord = new Coordinate(worldX, worldY);
+                Geometry centroidGeom = geometryFactory.createPoint(centroidCoord);
+                pixelRecords.add(new PixelRecord(centroidGeom, pixelValue, x, y));
+            }
+        }
+
+        return pixelRecords;
+    }
+
     public static Geometry getPixelAsPoint(GridCoverage2D raster, int colX, int rowY) throws TransformException, FactoryException {
         int srid = RasterAccessors.srid(raster);
         Point2D point2D = RasterUtils.getWorldCornerCoordinatesWithRangeCheck(raster, colX, rowY);

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

@@ -0,0 +1,16 @@
+package org.apache.sedona.common.raster;
+
+import org.locationtech.jts.geom.Geometry;
+
+public class PixelRecord {
+    public final Geometry geom;
+    public final double value;
+    public final int colX, rowY;
+
+    public PixelRecord(Geometry geom, double value, int colX, int rowY) {
+        this.geom = geom;
+        this.value = value;
+        this.colX = colX;
+        this.rowY = rowY;
+    }
+}

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

@@ -122,6 +122,24 @@ public void testPixelAsCentroid() throws FactoryException, TransformException {
         assertEquals(expected, actual);
     }
 
+    @Test
+    public void testPixelAsCentroids() throws FactoryException, TransformException {
+        GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 12, 13, 134, -53, 9);
+        List<PixelRecord> points = PixelFunctions.getPixelAsCentroids(emptyRaster, 1);
+        String expected = "POINT (156.5 -75.5)";
+        PixelRecord point = points.get(26);
+        Geometry geom = point.geom;
+        assertEquals(expected, geom.toString());
+
+        // Testing with skewed raster
+        emptyRaster = RasterConstructors.makeEmptyRaster(1, 12, 13, 240, -193, 2, 1.5, 3, 2, 0);
+        points = PixelFunctions.getPixelAsCentroids(emptyRaster, 1);
+        expected = "POINT (252.5 -184.25)";
+        point = points.get(26);
+        geom = point.geom;
+        assertEquals(expected, geom.toString());
+    }
+
     @Test
     public void testPixelAsPointUpperLeft() throws FactoryException, TransformException {
         GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, 123, -230, 8);

docs/api/sql/Raster-operators.md

@@ -22,6 +22,54 @@ Output:
 POINT (156.5 -75.5)
 ```
 
+### RS_PixelAsCentroids
+Introduction: Returns a list of the centroid point geometry, the pixel value and its raster X and Y coordinates for each pixel in the raster at the specified band.
+Each centroid represents the geometric center of the corresponding pixel's area.
+
+Format: `RS_PixelAsCentroids(raster: Raster, band: Integer)`
+
+Since: `v1.5.1`
+
+Spark SQL Example:
+```sql
+SELECT ST_AsText(RS_PixelAsCentroids(raster, 1)) from rasters
+```
+
+Output:
+```
+[[POINT (-13065222 4021263.75),148.0,0,0], [POINT (-13065151 4021263.75),123.0,0,1], [POINT (-13065077 4021263.75),99.0,1,0], [POINT (-13065007 4021261.75),140.0,1,1]]
+```
+
+Spark SQL example for extracting Point, value, raster x and y coordinates:
+
+```scala
+val pointDf = sedona.read...
+val rasterDf = sedona.read.format("binaryFile").load("/some/path/*.tiff")
+var df = sedona.read.format("binaryFile").load("/some/path/*.tiff")
+df = df.selectExpr("RS_FromGeoTiff(content) as raster")
+
+df.selectExpr(
+  "explode(RS_PixelAsCentroids(raster, 1)) as exploded"
+).selectExpr(
+  "exploded.geom as geom",
+  "exploded.value as value",
+  "exploded.x as x",
+  "exploded.y as y"
+).show(3)
+```
+
+Output:
+
+```
++----------------------------------------------+-----+---+---+
+|geom                                          |value|x  |y  |
++----------------------------------------------+-----+---+---+
+|POINT (-13095781.835693639 4021226.5856936392)|0.0  |1  |1  |
+|POINT (-13095709.507080918 4021226.5856936392)|0.0  |2  |1  |
+|POINT (-13095637.178468198 4021226.5856936392)|0.0  |3  |1  |
++----------------------------------------------+-----+---+---+
+```
+
 ### RS_PixelAsPoint
 
 Introduction: Returns a point geometry of the specified pixel's upper-left corner. The pixel coordinates specified are 1-indexed.

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

@@ -231,6 +231,7 @@ object Catalog {
     function[RS_PixelAsPoint](),
     function[RS_PixelAsPolygon](),
     function[RS_PixelAsCentroid](),
+    function[RS_PixelAsCentroids](),
     function[RS_Count](),
     function[RS_Clip](),
     function[RS_Band](),

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

@@ -18,10 +18,19 @@
  */
 package org.apache.spark.sql.sedona_sql.expressions.raster
 
-import org.apache.sedona.common.raster.PixelFunctions
-import org.apache.spark.sql.catalyst.expressions.Expression
+import org.apache.sedona.common.raster.{PixelFunctions}
+import org.apache.sedona.sql.utils.GeometrySerializer
+import org.apache.spark.sql.catalyst.InternalRow
+import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression}
+import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback
+import org.apache.spark.sql.catalyst.util.GenericArrayData
+import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT}
 import org.apache.spark.sql.sedona_sql.expressions.InferrableFunctionConverter._
 import org.apache.spark.sql.sedona_sql.expressions.InferredExpression
+import org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer
+import org.apache.spark.sql.types.{AbstractDataType, ArrayType, DataType, DoubleType, IntegerType, StructType}
+
+import scala.collection.convert.ImplicitConversions.`collection AsScalaIterable`
 
 case class RS_Value(inputExpressions: Seq[Expression]) extends InferredExpression(
   inferrableFunction2(PixelFunctions.value), inferrableFunction3(PixelFunctions.value), inferrableFunction4(PixelFunctions.value)
@@ -49,6 +58,43 @@ case class RS_PixelAsCentroid(inputExpressions: Seq[Expression]) extends Inferre
   }
 }
 
+case class RS_PixelAsCentroids(inputExpressions: Seq[Expression])
+  extends Expression with CodegenFallback with ExpectsInputTypes {
+
+  override def nullable: Boolean = true
+
+  override def dataType: DataType = ArrayType(new StructType()
+    .add("geom", GeometryUDT)
+    .add("value", DoubleType)
+    .add("x", IntegerType)
+    .add("y", IntegerType))
+
+  override def eval(input: InternalRow): Any = {
+    val rasterGeom = inputExpressions(0).toRaster(input)
+    val band = inputExpressions(1).eval(input).asInstanceOf[Int]
+
+    if (rasterGeom == null) {
+      null
+    } else {
+      val pixelRecords = PixelFunctions.getPixelAsCentroids(rasterGeom, band)
+      val rows = pixelRecords.map { pixelRecord =>
+        val serializedGeom = GeometrySerializer.serialize(pixelRecord.geom)
+        val rowArray = Array[Any](serializedGeom, pixelRecord.value, pixelRecord.colX, pixelRecord.rowY)
+        InternalRow.fromSeq(rowArray)
+      }
+      new GenericArrayData(rows.toArray)
+    }
+  }
+
+  override def children: Seq[Expression] = inputExpressions
+
+  protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]): RS_PixelAsCentroids = {
+    copy(inputExpressions = newChildren)
+  }
+
+  override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType)
+}
+
 case class RS_Values(inputExpressions: Seq[Expression]) extends InferredExpression(
   inferrableFunction2(PixelFunctions.values), inferrableFunction3(PixelFunctions.values), inferrableFunction4(PixelFunctions.values)
 ) {

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

@@ -955,6 +955,34 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen
       assertEquals(expected, result)
     }
 
+    it("Passed RS_PixelAsCentroids with empty raster") {
+      val widthInPixel = 12
+      val heightInPixel = 13
+      val upperLeftX = 240
+      val upperLeftY = -193
+      val cellSize = 9
+      val numBands = 2
+      val result = sparkSession.sql(s"SELECT RS_PixelAsCentroids(RS_MakeEmptyRaster($numBands, $widthInPixel, $heightInPixel, $upperLeftX, $upperLeftY, $cellSize), 1)").first().getList(0);
+      val expected = "[POINT (253.5 -215.5),0.0,2,3]"
+      assertEquals(expected, result.get(25).toString)
+    }
+
+    it("Passed RS_PixelAsCentroids with raster") {
+      var df = sparkSession.read.format("binaryFile").load(resourceFolder + "raster/test1.tiff")
+      df = df.selectExpr("RS_FromGeoTiff(content) as raster")
+
+      var result = df.selectExpr(
+        "explode(RS_PixelAsCentroids(raster, 1)) as exploded"
+      ).selectExpr(
+        "exploded.geom as geom",
+        "exploded.value as value",
+        "exploded.x as x",
+        "exploded.y as y"
+      )
+
+      assert(result.count() == 264704)
+    }
+
     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")