Skip to content

Commit

Permalink
[SEDONA-435] Add RS_PixelAsPolygons (#1131)
Browse files Browse the repository at this point in the history
  • Loading branch information
prantogg authored Nov 23, 2023
1 parent 6478129 commit 29b89f1
Show file tree
Hide file tree
Showing 6 changed files with 186 additions and 1 deletion.
Original file line number Diff line number Diff line change
Expand Up @@ -76,6 +76,56 @@ public static Geometry getPixelAsPolygon(GridCoverage2D raster, int colX, int ro
return GEOMETRY_FACTORY.createPolygon(coordinateArray);
}

public static List<PixelRecord> 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<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 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();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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<PixelRecord> 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);
Expand Down
51 changes: 50 additions & 1 deletion docs/api/sql/Raster-operators.md
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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](),
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down

0 comments on commit 29b89f1

Please sign in to comment.