From a282ed0afbeaee0b3d18427ba3fb76b665ab63b2 Mon Sep 17 00:00:00 2001 From: Pranav Toggi Date: Wed, 22 Nov 2023 12:10:37 -0500 Subject: [PATCH] [SEDONA-431] Add RS_PixelAsPoints (#1125) --- .../sedona/common/raster/PixelFunctions.java | 40 ++++++++- .../sedona/common/raster/PixelRecord.java | 16 ++++ .../sedona/common/raster/FunctionsTest.java | 87 +++++++++++++++++++ docs/api/sql/Raster-operators.md | 48 ++++++++++ .../org/apache/sedona/sql/UDF/Catalog.scala | 1 + .../expressions/raster/PixelFunctions.scala | 52 ++++++++++- .../apache/sedona/sql/rasteralgebraTest.scala | 34 +++++++- 7 files changed, 272 insertions(+), 6 deletions(-) create mode 100644 common/src/main/java/org/apache/sedona/common/raster/PixelRecord.java 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 c4912d0ee3..c43537b2f8 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 @@ -22,6 +22,7 @@ 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; @@ -29,10 +30,10 @@ 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; -import java.util.stream.Collectors; public class PixelFunctions { @@ -91,6 +92,43 @@ public static Geometry getPixelAsPoint(GridCoverage2D raster, int colX, int rowY return GEOMETRY_FACTORY.createPoint(pointCoord); } + public static List getPixelAsPoints(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 pointRecords = 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 - 1) * cellSizeX + (y - 1) * shearX; + double worldY = upperLeft.getY() + (y - 1) * cellSizeY + (x - 1) * shearY; + + Coordinate pointCoord = new Coordinate(worldX, worldY); + Geometry pointGeom = geometryFactory.createPoint(pointCoord); + pointRecords.add(new PixelRecord(pointGeom, pixelValue, x, y)); + } + } + + return pointRecords; + } + public static List values(GridCoverage2D rasterGeom, int[] xCoordinates, int[] yCoordinates, int band) throws TransformException { RasterUtils.ensureBand(rasterGeom, band); // Check for invalid band index int numBands = rasterGeom.getNumSampleDimensions(); diff --git a/common/src/main/java/org/apache/sedona/common/raster/PixelRecord.java b/common/src/main/java/org/apache/sedona/common/raster/PixelRecord.java new file mode 100644 index 0000000000..58bdfae041 --- /dev/null +++ b/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; + } +} 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 d141944eb1..34504676df 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 @@ -183,6 +183,93 @@ public void testPixelAsPointFromRasterFile() throws IOException, TransformExcept assertEquals(expectedY, coordinate.getY(), 0.2d); } + @Test + public void testPixelAsPointSkewedRaster() throws FactoryException, TransformException { + // Testing with skewed raster + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 12, 13, 240, -193, 2, 1.5, 3, 2, 0); + String actual = Functions.asWKT(PixelFunctions.getPixelAsPoint(emptyRaster, 3, 3)); + String expected = "POINT (250 -186)"; + assertEquals(expected, actual); + } + + @Test + public void testPixelAsPointsOutputSize() throws FactoryException, TransformException { + GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 5, 10, 123, -230, 8); + List points = PixelFunctions.getPixelAsPoints(raster, 1); + assertEquals(50, points.size()); + } + + @Test + public void testPixelAsPointsValues() throws FactoryException, TransformException { + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, 123, -230, 8); + List points = PixelFunctions.getPixelAsPoints(emptyRaster, 1); + + PixelRecord point1 = points.get(0); + Geometry geom1 = (Geometry) point1.geom; + assertEquals(123, geom1.getCoordinate().x, FP_TOLERANCE); + assertEquals(-230, geom1.getCoordinate().y, FP_TOLERANCE); + assertEquals(0.0, point1.value, FP_TOLERANCE); + + PixelRecord point2 = points.get(22); + Geometry geom2 = (Geometry) point2.geom; + assertEquals(139, geom2.getCoordinate().x, FP_TOLERANCE); + assertEquals(-262, geom2.getCoordinate().y, FP_TOLERANCE); + assertEquals(0.0, point2.value, FP_TOLERANCE); + } + + @Test + public void testPixelAsPointsCustomSRIDPlanar() throws FactoryException, TransformException { + int srid = 3857; + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 5, -123, 54, 5, 5, 0, 0, srid); + List points = PixelFunctions.getPixelAsPoints(emptyRaster, 1); + PixelRecord point1 = points.get(0); + Geometry geom1 = (Geometry) point1.geom; + assertEquals(-123, geom1.getCoordinate().x, FP_TOLERANCE); + assertEquals(54, geom1.getCoordinate().y, FP_TOLERANCE); + assertEquals(srid, geom1.getSRID()); + } + + @Test + public void testPixelAsPointsSRIDSpherical() throws FactoryException, TransformException { + int srid = 4326; + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 5, 10, -123, 54, 5, -10, 0, 0, srid); + List points = PixelFunctions.getPixelAsPoints(emptyRaster, 1); + + PixelRecord point1 = points.get(11); + Geometry geom1 = (Geometry) point1.geom; + assertEquals(-118, geom1.getCoordinate().x, FP_TOLERANCE); + assertEquals(34, geom1.getCoordinate().y, FP_TOLERANCE); + assertEquals(srid, geom1.getSRID()); + } + + @Test + public void testPixelAsPointsFromRasterFile() throws IOException, TransformException, FactoryException { + GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster/test1.tiff"); + List points = PixelFunctions.getPixelAsPoints(raster, 1); + PixelRecord firstPoint = points.get(0); + Geometry firstGeom = (Geometry) firstPoint.geom; + + double expectedX = -1.3095818E7; + double expectedY = 4021262.75; + double val = 0.0; + + assertEquals(expectedX, firstGeom.getCoordinate().x, FP_TOLERANCE); + assertEquals(expectedY, firstGeom.getCoordinate().y, FP_TOLERANCE); + assertEquals(val, firstPoint.value, FP_TOLERANCE); + } + + @Test + public void testPixelAsPointsSkewedRaster() throws FactoryException, TransformException { + // Testing with skewed raster + GridCoverage2D emptyRaster = RasterConstructors.makeEmptyRaster(1, 12, 13, 240, -193, 2, 1.5, 3, 2, 0); + List points = PixelFunctions.getPixelAsPoints(emptyRaster, 1); + + PixelRecord point1 = points.get(26); + Geometry geom1 = (Geometry) point1.geom; + String expected = "POINT (250 -186)"; + assertEquals(expected, geom1.toString()); + } + @Test public void values() throws TransformException { // The function 'value' is implemented using 'values'. diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md index a002eb046f..2782b4d56b 100644 --- a/docs/api/sql/Raster-operators.md +++ b/docs/api/sql/Raster-operators.md @@ -56,6 +56,54 @@ Output: IndexOutOfBoundsException: Specified pixel coordinates (6, 2) do not lie in the raster ``` +### RS_PixelAsPoints +Introduction: Returns a list of the pixel's upper-left corner point geometry, the pixel value and its raster X and Y coordinates for each pixel in the raster at the specified band. + +Format: `RS_PixelAsPoints(raster: Raster, band: Integer)` + +Since: `v1.5.1` + +Spark SQL Example: +```sql +SELECT ST_AsText(RS_PixelAsPoints(raster, 1)) from rasters +``` + +Output: +``` +[[POINT (-13065223 4021262.75),148.0,0,0], [POINT (-13065150 4021262.75),123.0,0,1], [POINT (-13065078 4021262.75),99.0,1,0], [POINT (-13065006 4021262.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_PixelAsPoints(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 (-13095818 4021262.75) |0.0 |1 |1 | +|POINT (-13095745.67138728 4021262.75) |0.0 |2 |1 | +|POINT (-13095673.342774557 4021262.75)|0.0 |3 |1 | ++--------------------------------------+-----+---+---+ +``` + ### RS_PixelAsPolygon Introduction: Returns a polygon geometry that bounds the specified pixel. 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 718b014ffe..c219035ce1 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 @@ -229,6 +229,7 @@ object Catalog { function[RS_Rotation](), function[RS_GeoTransform](), function[RS_PixelAsPoint](), + function[RS_PixelAsPoints](), function[RS_PixelAsPolygon](), function[RS_PixelAsCentroid](), function[RS_Count](), 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 9cde85f710..4d4eb7f1d1 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 @@ -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) @@ -37,6 +46,45 @@ case class RS_PixelAsPoint(inputExpressions: Seq[Expression]) extends InferredEx } } +case class RS_PixelAsPoints(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.getPixelAsPoints(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_PixelAsPoints = { + copy(inputExpressions = newChildren) + } + + override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, IntegerType) +} + + + case class RS_PixelAsPolygon(inputExpressions: Seq[Expression]) extends InferredExpression(PixelFunctions.getPixelAsPolygon _) { protected def withNewChildrenInternal(newChilren: IndexedSeq[Expression]) = { copy(inputExpressions = newChilren) 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 2493bda2b2..1f5ad5fca3 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 @@ -20,11 +20,11 @@ package org.apache.sedona.sql import org.apache.sedona.common.raster.MapAlgebra import org.apache.sedona.common.utils.RasterUtils -import org.apache.spark.sql.SaveMode -import org.apache.spark.sql.functions.{collect_list, expr} +import org.apache.spark.sql.{Row, SaveMode} +import org.apache.spark.sql.functions.{col, collect_list, expr} import org.geotools.coverage.grid.GridCoverage2D import org.junit.Assert.{assertEquals, assertNull, assertTrue} -import org.locationtech.jts.geom.{Coordinate, Geometry} +import org.locationtech.jts.geom.{Coordinate, Geometry, Point} import org.scalatest.{BeforeAndAfter, GivenWhenThen} import java.awt.image.DataBuffer @@ -985,6 +985,34 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen assertEquals(expectedY, actualCoordinates.y, 1e-5) } + it("Passed RS_PixelAsPoints 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_PixelAsPoints(RS_MakeEmptyRaster($numBands, $widthInPixel, $heightInPixel, $upperLeftX, $upperLeftY, $cellSize), 1)").first().getList(0); + val expected = "[POINT (127.19000244140625 -12),0.0,2,1]" + assertEquals(expected, result.get(1).toString) + } + + it("Passed RS_PixelAsPoints 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_PixelAsPoints(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_PixelAsPolygon with raster") { val widthInPixel = 5 val heightInPixel = 10