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 c43537b2f8..d46ddcf9bb 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 @@ -76,6 +76,56 @@ public static Geometry getPixelAsPolygon(GridCoverage2D raster, int colX, int ro return GEOMETRY_FACTORY.createPolygon(coordinateArray); } + public static List getPixelAsPolygons(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 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 worldX1 = upperLeft.getX() + (x - 1) * cellSizeX + (y - 1) * shearX; + double worldY1 = upperLeft.getY() + (y - 1) * cellSizeY + (x - 1) * shearY; + double worldX2 = worldX1 + cellSizeX; + double worldY2 = worldY1 + shearY; + double worldX3 = worldX2 + shearX; + double worldY3 = worldY2 + cellSizeY; + double worldX4 = worldX1 + shearX; + double worldY4 = worldY1 + cellSizeY; + + Coordinate[] coordinates = new Coordinate[] { + new Coordinate(worldX1, worldY1), + new Coordinate(worldX2, worldY2), + new Coordinate(worldX3, worldY3), + new Coordinate(worldX4, worldY4), + new Coordinate(worldX1, worldY1) + }; + + Geometry polygon = geometryFactory.createPolygon(coordinates); + pixelRecords.add(new PixelRecord(polygon, pixelValue, x, y)); + } + } + + return pixelRecords; + } + public static Geometry getPixelAsCentroid(GridCoverage2D raster, int colX, int rowY) throws FactoryException, TransformException { Geometry polygon = PixelFunctions.getPixelAsPolygon(raster, colX, rowY); return polygon.getCentroid(); 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 34504676df..2df5e8f03f 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 @@ -108,6 +108,25 @@ public void testPixelAsPolygon() throws FactoryException, TransformException { assertEquals(expected,actual); } + @Test + public void testPixelAsPointsPolygons() throws FactoryException, TransformException { + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, 123, -230, 8); + List points = PixelFunctions.getPixelAsPolygons(emptyRaster, 1); + + PixelRecord point = points.get(11); + Geometry geom = (Geometry) point.geom; + String expected = "POLYGON ((131 -246, 139 -246, 139 -254, 131 -254, 131 -246))"; + assertEquals(expected,geom.toString()); + + // Testing with skewed raster + emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, 234, -43, 3, 4, 2,3,0); + points = PixelFunctions.getPixelAsPolygons(emptyRaster, 1); + point = points.get(11); + geom = (Geometry) point.geom; + expected = "POLYGON ((241 -32, 244 -29, 246 -25, 243 -28, 241 -32))"; + assertEquals(expected,geom.toString()); + } + @Test public void testPixelAsCentroid() throws FactoryException, TransformException { GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 12, 13, 134, -53, 9); diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md index 2782b4d56b..0e5b3dcbd7 100644 --- a/docs/api/sql/Raster-operators.md +++ b/docs/api/sql/Raster-operators.md @@ -13,7 +13,7 @@ Since: `v1.5.0` Spark SQL Example: ```sql -SELECT ST_AsText(RS_PixelAsPolygon(RS_MakeEmptyRaster(1, 12, 13, 134, -53, 9), 3, 3)) +SELECT ST_AsText(RS_PixelAsCentroid(RS_MakeEmptyRaster(1, 12, 13, 134, -53, 9), 3, 3)) ``` Output: @@ -126,6 +126,55 @@ Output: POLYGON ((131 -246, 139 -246, 139 -254, 131 -254, 131 -246)) ``` +### RS_PixelAsPolygons +Introduction: Returns a list of the polygon geometry, the pixel value and its raster X and Y coordinates for each pixel in the raster at the specified band. + +Format: `RS_PixelAsPolygons(raster: Raster, band: Integer)` + +Since: `v1.5.1` + +Spark SQL Example: +```sql +SELECT ST_AsText(RS_PixelAsPolygons(raster, 1)) from rasters +``` + +Output: +``` +[[POLYGON ((123.19000244140625 -12, 127.19000244140625 -12, 127.19000244140625 -16, 123.19000244140625 -16, 123.19000244140625 -12)),0.0,1,1], +[POLYGON ((127.19000244140625 -12, 131.19000244140625 -12, 131.19000244140625 -16, 127.19000244140625 -16, 127.19000244140625 -12)),0.0,2,1], +[POLYGON ((131.19000244140625 -12, 135.19000244140625 -12, 135.19000244140625 -16, 131.19000244140625 -16, 131.19000244140625 -12)),0.0,3,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_PixelAsPolygons(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| ++--------------------+-----+---+---+ +|POLYGON ((-130958...| 0.0| 1| 1| +|POLYGON ((-130957...| 0.0| 2| 1| +|POLYGON ((-130956...| 0.0| 3| 1| ++--------------------+-----+---+---+ +``` + ## Geometry Functions ### RS_Envelope 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 c219035ce1..d930c68181 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 @@ -231,6 +231,7 @@ object Catalog { function[RS_PixelAsPoint](), function[RS_PixelAsPoints](), function[RS_PixelAsPolygon](), + function[RS_PixelAsPolygons](), function[RS_PixelAsCentroid](), function[RS_Count](), function[RS_Clip](), diff --git a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala index 4d4eb7f1d1..562bd1f467 100644 --- a/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala +++ b/spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/PixelFunctions.scala @@ -91,6 +91,43 @@ case class RS_PixelAsPolygon(inputExpressions: Seq[Expression]) extends Inferred } } +case class RS_PixelAsPolygons(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.getPixelAsPolygons(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_PixelAsPolygons = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType) +} + case class RS_PixelAsCentroid(inputExpressions: Seq[Expression]) extends InferredExpression(PixelFunctions.getPixelAsCentroid _) { protected def withNewChildrenInternal(newChildren: IndexedSeq[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 1f5ad5fca3..4fd5592bfb 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 @@ -1025,6 +1025,35 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen assertEquals(expected, result) } + it("Passed RS_PixelAsPolygons with empty raster") { + val widthInPixel = 5 + val heightInPixel = 10 + val upperLeftX = 123.19 + val upperLeftY = -12 + val cellSize = 4 + val numBands = 2 + val result = sparkSession.sql(s"SELECT RS_PixelAsPolygons(RS_MakeEmptyRaster($numBands, $widthInPixel, $heightInPixel, $upperLeftX, $upperLeftY, $cellSize), 1)").first().getList(0); + val expected = "[POLYGON ((127.19000244140625 -20, 131.19000244140625 -20, 131.19000244140625 -24, 127.19000244140625 -24, 127.19000244140625 -20)),0.0,2,3]" + assertEquals(expected, result.get(11).toString) + } + + + it("Passed RS_PixelAsPolygons 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_PixelAsPolygons(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_PixelAsCentroid with raster") { val widthInPixel = 12 val heightInPixel = 13