From a616b0f31011c943a4928fbabf75f47c1ec21d02 Mon Sep 17 00:00:00 2001 From: Pranav Toggi Date: Fri, 15 May 2026 17:59:08 -0700 Subject: [PATCH] feat: make rasterization orientation-agnostic by supporting positive scaleY --- .../sedona/common/raster/Rasterization.java | 88 ++++++--- .../sedona/common/utils/RasterUtils.java | 166 ++++++++++++++++ .../raster/RasterBandAccessorsTest.java | 21 ++ .../common/raster/RasterBandEditorsTest.java | 143 +++++++++++++- .../common/raster/RasterConstructorsTest.java | 181 ++++++++++++++++++ .../common/raster/RasterizationTests.java | 160 +++++++++++++++- .../sedona/common/utils/RasterUtilsTest.java | 129 +++++++++++++ docs/api/sql/Raster-Operators/RS_AsRaster.md | 1 - .../rasterization/test_rasterization.tiff | Bin 524672 -> 1218 bytes 9 files changed, 860 insertions(+), 29 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java index f5f32aecffb..b70b8e9665d 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java +++ b/common/src/main/java/org/apache/sedona/common/raster/Rasterization.java @@ -18,6 +18,9 @@ */ package org.apache.sedona.common.raster; +import static org.apache.sedona.common.utils.RasterUtils.getRasterScaleY; +import static org.apache.sedona.common.utils.RasterUtils.getRasterUpperLeftY; + import java.awt.image.WritableRaster; import java.math.BigDecimal; import java.math.RoundingMode; @@ -69,6 +72,11 @@ protected static List rasterize( rasterized = coverageFactory.create("rasterized", params.writableRaster, rasterExtent); } + // If original raster is bottom-up, flip the rasterized result vertically + if (params.bottomUp) { + rasterized = RasterUtils.flipVerticallyGridSpace(rasterized); + } + // Return results compatible with the original function List objects = new ArrayList<>(); objects.add(params.writableRaster); @@ -136,16 +144,22 @@ private static void rasterizePoint( int x = startX; int y = startY; - for (double worldY = geomExtent.getMinY(); - worldY < geomExtent.getMaxY(); - worldY += params.scaleY, y++) { + for (double worldY = geomExtent.getMaxY(); + worldY > geomExtent.getMinY(); + worldY += params.scaleY, y--) { x = startX; for (double worldX = geomExtent.getMinX(); worldX < geomExtent.getMaxX(); worldX += params.scaleX, x++) { - // Flip y-axis (since raster Y starts from top-left) - int yIndex = -y - 1; + // Make zero indexed + int yIndex = y - 1; + + // Adjust for bottom-up rasters + // Reverse the y index + if (params.bottomUp) { + yIndex = params.writableRaster.getHeight() - 1 - yIndex; + } // Create envelope for this pixel double cellMaxX = worldX + params.scaleX; @@ -182,9 +196,9 @@ private static void rasterizeLineString( } double x0 = (start.x - params.upperLeftX) / params.scaleX; - double y0 = (params.upperLeftY - start.y) / params.scaleY; + double y0 = (start.y - params.upperLeftY) / params.scaleY; double x1 = (end.x - params.upperLeftX) / params.scaleX; - double y1 = (params.upperLeftY - end.y) / params.scaleY; + double y1 = (end.y - params.upperLeftY) / params.scaleY; // Apply Bresenham for this segment drawLineBresenham(params, x0, y0, x1, y1, value, 0.2); @@ -221,6 +235,12 @@ private static void drawLineBresenham( int rasterX = (int) (Math.floor(x)); int rasterY = (int) (Math.floor(y)); + // Adjust for bottom-up rasters + // Reverse the y index + if (params.bottomUp) { + rasterY = params.writableRaster.getHeight() - 1 - rasterY; + } + // Only write if within raster bounds if (rasterX >= 0 && rasterX < params.writableRaster.getWidth() @@ -338,9 +358,9 @@ static ReferencedEnvelope rasterizeGeomExtent( // Using BigDecimal to avoid floating point errors double upperLeftX = metadata[0]; - double upperLeftY = metadata[1]; + double upperLeftY = getRasterUpperLeftY(metadata); double scaleX = metadata[4]; - double scaleY = metadata[5]; + double scaleY = getRasterScaleY(metadata); // Compute the aligned min/max values double alignedMinX = @@ -464,13 +484,14 @@ private static RasterizationParams calculateRasterizationParams( upperLeftY = geomExtent.getMaxY(); } else { upperLeftX = metadata[0]; - upperLeftY = metadata[1]; + upperLeftY = getRasterUpperLeftY(metadata); } WritableRaster writableRaster; if (useGeometryExtent) { int geomExtentWidth = (int) (Math.round(geomExtent.getWidth() / metadata[4])); - int geomExtentHeight = (int) (Math.round(geomExtent.getHeight() / -metadata[5])); + int geomExtentHeight = + (int) (Math.round(geomExtent.getHeight() / -getRasterScaleY(metadata))); writableRaster = RasterFactory.createBandedRaster( @@ -483,19 +504,30 @@ private static RasterizationParams calculateRasterizationParams( RasterUtils.getDataTypeCode(pixelType), rasterWidth, rasterHeight, 1, null); } + boolean bottomUp = metadata[5] > 0; + return new RasterizationParams( - writableRaster, pixelType, metadata[4], -metadata[5], upperLeftX, upperLeftY); + writableRaster, + pixelType, + metadata[4], + getRasterScaleY(metadata), + upperLeftX, + upperLeftY, + bottomUp); } private static void validateRasterMetadata(double[] rasterMetadata) { if (rasterMetadata[4] < 0) { - throw new IllegalArgumentException("ScaleX cannot be negative"); - } - if (rasterMetadata[5] > 0) { - throw new IllegalArgumentException("ScaleY must be negative"); + throw new IllegalArgumentException( + String.format( + "Invalid raster metadata: scaleX is negative (%.2f). Right-to-left rasters are not supported.", + rasterMetadata[4])); } if (rasterMetadata[6] != 0 || rasterMetadata[7] != 0) { - throw new IllegalArgumentException("SkewX and SkewY must be zero"); + throw new IllegalArgumentException( + String.format( + "Invalid raster metadata: skewX is %.2f and skewY is %.2f. Both values must be zero.", + rasterMetadata[6], rasterMetadata[7])); } } @@ -506,6 +538,7 @@ private static class RasterizationParams { double scaleY; double upperLeftX; double upperLeftY; + boolean bottomUp; RasterizationParams( WritableRaster writableRaster, @@ -513,13 +546,15 @@ private static class RasterizationParams { double scaleX, double scaleY, double upperLeftX, - double upperLeftY) { + double upperLeftY, + boolean bottomUp) { this.writableRaster = writableRaster; this.pixelType = pixelType; this.scaleX = scaleX; this.scaleY = scaleY; this.upperLeftX = upperLeftX; this.upperLeftY = upperLeftY; + this.bottomUp = bottomUp; } } @@ -602,15 +637,15 @@ private static Map> computeScanlineIntersections( // Using BigDecimal to avoid floating point errors double yStart = Math.round( - (BigDecimal.valueOf(params.upperLeftY) - .subtract(BigDecimal.valueOf(worldP1.y)) + (BigDecimal.valueOf(worldP1.y) + .subtract(BigDecimal.valueOf(params.upperLeftY)) .divide(BigDecimal.valueOf(params.scaleY), RoundingMode.CEILING)) .doubleValue()); double yEnd = Math.round( - (BigDecimal.valueOf(params.upperLeftY) - .subtract(BigDecimal.valueOf(worldP2.y)) + (BigDecimal.valueOf(worldP2.y) + .subtract(BigDecimal.valueOf(params.upperLeftY)) .divide(BigDecimal.valueOf(params.scaleY), RoundingMode.FLOOR)) .doubleValue()); @@ -619,7 +654,7 @@ private static Map> computeScanlineIntersections( yStart = Math.min((params.writableRaster.getHeight() - 0.5), Math.abs(yStart) - 0.5); double p1X = (worldP1.x - params.upperLeftX) / params.scaleX; - double p1Y = (params.upperLeftY - worldP1.y) / params.scaleY; + double p1Y = (worldP1.y - params.upperLeftY) / params.scaleY; if (worldP1.x == worldP2.x) { // Vertical line case: directly set xIntercept. Avoid divide by zero error when @@ -695,6 +730,13 @@ private static void fillPolygon( for (Map.Entry> entry : scanlineFillRanges.entrySet()) { int y = entry.getKey(); + + // Adjust for bottom-up rasters + // Reverse the y index + if (params.bottomUp) { + y = params.writableRaster.getHeight() - 1 - y; + } + for (int[] range : entry.getValue()) { if (range.length == 1) { params.writableRaster.setSample(range[0], y, 0, value); diff --git a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java b/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java index f962bb92605..6113ee943d5 100644 --- a/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java +++ b/common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java @@ -984,4 +984,170 @@ public static GridCoverage2D applyRasterMask(GridCoverage2D raster, GridCoverage false); return modifiedRaster; } + + /** + * This function calculates the upper-left Y world coordinate of the raster. It accounts for both + * bottom-up and top-down raster conventions. + * + *

For bottom-up rasters (metadata[5] > 0), the upper-left Y world coordinate is directly taken + * from metadata[1]. - For top-down rasters (metadata[5] < 0), the upper-left Y coordinate is + * adjusted by adding the product of the raster height and the scale factor (metadata[5]). + * + * @param metadata The metadata array containing raster properties. + * @return The upper-left Y world coordinate of the raster. + */ + public static double getRasterUpperLeftY(double[] metadata) { + if (metadata[5] < 0) { + return metadata[1]; + } else { + int rasterHeight = (int) metadata[3]; + return metadata[1] + (rasterHeight * metadata[5]); + } + } + + /** + * This function calculates the Y scale factor of the raster. It ensures that the scale factor is + * always negative, regardless of the raster convention. + * + *

For bottom-up rasters (metadata[5] > 0), the scale factor is returned as-is. - For top-down + * rasters (metadata[5] < 0), the scale factor is negated to maintain consistency. + * + * @param metadata The metadata array containing raster properties. + * @return The Y scale factor of the raster. + */ + public static double getRasterScaleY(double[] metadata) { + if (metadata[5] < 0) { + return metadata[5]; + } else { + return -metadata[5]; + } + } + + /** + * GRID SPACE: Grid space refers to the raster's index space (col, row). Pixel values remain at + * the same array indices, but the grid-to-CRS affine transform is rewritten. + * + *

EFFECT: - The raster is reinterpreted spatially (top-down ↔ bottom-up). - Pixel values are + * NOT moved in memory. - Pixel world coordinates DO change. - The world envelope (bounding box) + * is preserved. + * + *

Use this when you want to change how the raster is oriented in CRS space without touching + * the pixel data. + */ + public static GridCoverage2D flipVerticallyGridSpace(GridCoverage2D raster) { + + GridGeometry2D gridGeometry = raster.getGridGeometry(); + + MathTransform mt = gridGeometry.getGridToCRS2D(); + if (!(mt instanceof AffineTransform2D)) { + throw new UnsupportedOperationException( + "flipVertically only supports AffineTransform2D grid-to-CRS transforms"); + } + + AffineTransform2D affine = (AffineTransform2D) mt; + + RenderedImage image = raster.getRenderedImage(); + int height = image.getHeight(); + + /* + * Affine transform definition (GeoTools / GDAL compatible): + * + * worldX = scaleX * col + shearX * row + translateX + * worldY = shearY * col + scaleY * row + translateY + */ + double scaleX = affine.getScaleX(); + double shearX = affine.getShearX(); + double shearY = affine.getShearY(); + double scaleY = affine.getScaleY(); + double translateX = affine.getTranslateX(); + double translateY = affine.getTranslateY(); + + /* + * Apply grid-space vertical flip: + * + * row = (H - 1) - row' + */ + double newShearX = -shearX; + double newScaleY = -scaleY; + + double newTranslateX = translateX + shearX * (height - 1); + double newTranslateY = translateY + scaleY * (height - 1); + + AffineTransform2D flippedAffine = + new AffineTransform2D(scaleX, shearY, newShearX, newScaleY, newTranslateX, newTranslateY); + + GridGeometry2D newGridGeometry = + new GridGeometry2D( + gridGeometry.getGridRange(), + flippedAffine, + gridGeometry.getCoordinateReferenceSystem()); + + return RasterUtils.clone( + raster.getRenderedImage(), + newGridGeometry, + raster.getSampleDimensions(), + raster, + null, + true); + } + + /** + * PIXEL SPACE: Pixel space refers to pixel values anchored at their world coordinates. Pixel + * values are physically rearranged so that each value remains at the same CRS location. + * + *

EFFECT: - Raster rows are reversed in memory. - The affine transform is updated accordingly. + * - Pixel world coordinates are preserved. - The world envelope (bounding box) is preserved. + * + *

Use this when pixel values must remain fixed at their geographic locations. + */ + public static GridCoverage2D flipVerticallyPixelSpace(GridCoverage2D raster) { + GridGeometry2D gridGeometry = raster.getGridGeometry(); + + MathTransform mt = gridGeometry.getGridToCRS2D(); + if (!(mt instanceof AffineTransform2D)) { + throw new UnsupportedOperationException( + "flipVertically only supports AffineTransform2D grid-to-CRS transforms"); + } + + AffineTransform2D affine = (AffineTransform2D) mt; + + RenderedImage image = raster.getRenderedImage(); + int height = image.getHeight(); + + // Flip the affine transformation + double scaleX = affine.getScaleX(); + double shearX = affine.getShearX(); + double shearY = affine.getShearY(); + double scaleY = affine.getScaleY(); + double translateX = affine.getTranslateX(); + double translateY = affine.getTranslateY(); + + double newShearX = -shearX; + double newScaleY = -scaleY; + double newTranslateX = translateX + shearX * (height - 1); + double newTranslateY = translateY + scaleY * (height - 1); + + AffineTransform2D flippedAffine = + new AffineTransform2D(scaleX, shearY, newShearX, newScaleY, newTranslateX, newTranslateY); + + // Create a new WritableRaster to preserve pixel values + Raster originalRaster = RasterUtils.getRaster(image); + WritableRaster flippedRaster = originalRaster.createCompatibleWritableRaster(); + + // Copy pixel values to the flipped raster + for (int y = 0; y < height; y++) { + for (int x = 0; x < image.getWidth(); x++) { + flippedRaster.setPixel(x, height - y - 1, originalRaster.getPixel(x, y, (double[]) null)); + } + } + + GridGeometry2D newGridGeometry = + new GridGeometry2D( + gridGeometry.getGridRange(), + flippedAffine, + gridGeometry.getCoordinateReferenceSystem()); + + return RasterUtils.clone( + flippedRaster, newGridGeometry, raster.getSampleDimensions(), raster, null, true); + } } diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java index 0eaea403416..a3d2a474ea7 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterBandAccessorsTest.java @@ -18,6 +18,7 @@ */ package org.apache.sedona.common.raster; +import static org.apache.sedona.common.utils.RasterUtils.flipVerticallyPixelSpace; import static org.junit.Assert.*; import java.io.IOException; @@ -127,6 +128,7 @@ public void testZonalStatsForVerticalLineBug() throws Exception { public void testZonalStats() throws FactoryException, ParseException, IOException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); String polygon = "POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170 4197590, 236722 4204770))"; Geometry geom = Constructors.geomFromWKT(polygon, RasterAccessors.srid(raster)); @@ -180,6 +182,15 @@ public void testZonalStats() throws FactoryException, ParseException, IOExceptio () -> RasterBandAccessors.getZonalStats( raster, nonIntersectingGeom, 1, "sum", false, false, false)); + + // Test bottom-up raster case + geom = Constructors.geomFromWKT(polygon, RasterAccessors.srid(raster)); + actual = RasterBandAccessors.getZonalStats(raster_bottom_up, geom, 1, "sum", false, false); + assertEquals(1.0795427E7, actual, 0d); + + actual = RasterBandAccessors.getZonalStats(raster_bottom_up, geom, 2, "mean", false, false); + expected = 220.7711988936036; + assertEquals(expected, actual, FP_TOLERANCE); } @Test @@ -206,6 +217,7 @@ public void testZonalStatsWithMultiPolygonAndAllNoDataValues() public void testZonalStatsWithNoData() throws IOException, FactoryException, ParseException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster/raster_with_no_data/test5.tiff"); + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); String polygon = "POLYGON((-167.750000 87.750000, -155.250000 87.750000, -155.250000 40.250000, -180.250000 40.250000, -167.750000 87.750000))"; // Testing implicit CRS transformation @@ -238,6 +250,15 @@ public void testZonalStatsWithNoData() throws IOException, FactoryException, Par actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sd", false, true); expected = 74.81287592054916; assertEquals(expected, actual, FP_TOLERANCE); + + // Test bottom-up raster case + actual = RasterBandAccessors.getZonalStats(raster_bottom_up, geom, 1, "sum", false, true); + expected = 3229013.0; + assertEquals(expected, actual, 0d); + + actual = RasterBandAccessors.getZonalStats(raster_bottom_up, geom, 1, "mean", false, true); + expected = 226.61330619692416; + assertEquals(expected, actual, FP_TOLERANCE); } @Test diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java index 9b637c54a40..fbd49d456a5 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterBandEditorsTest.java @@ -19,6 +19,7 @@ package org.apache.sedona.common.raster; import static org.apache.sedona.common.raster.RasterBandEditors.rasterUnion; +import static org.apache.sedona.common.utils.RasterUtils.flipVerticallyPixelSpace; import static org.junit.Assert.*; import java.io.IOException; @@ -180,6 +181,9 @@ public void testClipWithGeometryTransform() throws FactoryException, IOException, ParseException, TransformException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); + + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); + String polygon = "POLYGON ((-8682522.873537656 4572703.890837922, -8673439.664183248 4572993.532747675, -8673155.57366801 4563873.2099182755, -8701890.325907696 4562931.7093397, -8682522.873537656 4572703.890837922))"; Geometry geom = Constructors.geomFromWKT(polygon, 3857); @@ -199,6 +203,27 @@ public void testClipWithGeometryTransform() Double[] actualValues = PixelFunctions.values(clippedRaster, points, 1).toArray(new Double[0]); Double[] expectedValues = new Double[] {null, null, 0.0, 0.0, null}; assertTrue(Arrays.equals(expectedValues, actualValues)); + + // Test for bottom-up raster + clippedRaster = RasterBandEditors.clip(raster_bottom_up, 1, geom, false, 200, false); + clippedMetadata = Arrays.stream(RasterAccessors.metadata(clippedRaster), 0, 9).toArray(); + originalMetadata = Arrays.stream(RasterAccessors.metadata(raster_bottom_up), 0, 9).toArray(); + assertArrayEquals(originalMetadata, clippedMetadata, 0.01d); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(clippedRaster)[5], + 0.0); + + points = new ArrayList<>(); + points.add(Constructors.geomFromWKT("POINT(223802 4.21769e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(224759 4.20453e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237201 4.20429e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237919 4.20357e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(254668 4.21769e+06)", 26918)); + actualValues = PixelFunctions.values(clippedRaster, points, 1).toArray(new Double[0]); + assertTrue(Arrays.equals(expectedValues, actualValues)); } @Test @@ -206,9 +231,13 @@ public void testClip_smallAOI() throws FactoryException, TransformException, Par // Test for AOI geometries smaller than a pixel GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, "F", 4, 4, 401805.039562261, 2095852.150947876, 30); + raster = RasterEditors.setSrid(raster, 5070); - double[] bandValues = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}; - raster = MapAlgebra.addBandFromArray(raster, bandValues, 1); + double[] bandValues_top_down = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}; + raster = MapAlgebra.addBandFromArray(raster, bandValues_top_down, 1); + + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); + String polygon = "POLYGON ((401866.1071465613 2095803.8989834636, 401850.9983055725 2095803.1375789386, 401850.039562261 2095822.150947876, 401865.14840359957 2095822.9124532426, 401880.2572475523 2095823.6738547424, 401881.2159852499 2095804.6604855175, 401866.1071465613 2095803.8989834636))"; Geometry geom = Constructors.geomFromWKT(polygon, 5070); @@ -222,6 +251,22 @@ public void testClip_smallAOI() throws FactoryException, TransformException, Par assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); assertTrue(Arrays.equals(expectedValues, actualValues)); + // Test for bottom-up raster + clippedRaster = RasterBandEditors.clip(raster_bottom_up, 1, geom, false, 0, true); + bandNoDataValue = RasterBandAccessors.getBandNoDataValue(clippedRaster); + actualValues = MapAlgebra.bandAsArray(clippedRaster, 1); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(clippedRaster)[5], + 0.0); + + expectedValues = new double[] {0.0, 7.0, 0.0, 0.0}; + assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); + assertTrue(Arrays.equals(expectedValues, actualValues)); + + // Test for geometry in different CRS Geometry geomTransformed = FunctionsGeoTools.transform(geom, "EPSG:5070", "EPSG:4326"); clippedRaster = RasterBandEditors.clip(raster, 1, geomTransformed, false, 0, true); @@ -232,6 +277,21 @@ public void testClip_smallAOI() throws FactoryException, TransformException, Par assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); assertTrue(Arrays.equals(expectedValues, actualValues)); + + // Test for bottom-up raster + clippedRaster = RasterBandEditors.clip(raster_bottom_up, 1, geom, false, 0, true); + bandNoDataValue = RasterBandAccessors.getBandNoDataValue(clippedRaster); + actualValues = MapAlgebra.bandAsArray(clippedRaster, 1); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(clippedRaster)[5], + 0.0); + + expectedValues = new double[] {0.0, 7.0, 0.0, 0.0}; + assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); + assertTrue(Arrays.equals(expectedValues, actualValues)); } @Test @@ -279,6 +339,24 @@ public void testClipWithClippedGeometryWithHole() assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); assertArrayEquals(expectedValues, actualValues, 0.0); + + // Test for bottom-up raster + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); + clippedRaster = RasterBandEditors.clip(raster_bottom_up, 1, geom, true, 0, true); + + bandNoDataValue = RasterBandAccessors.getBandNoDataValue(clippedRaster); + expectedBandNoDataValue = 0.0; + actualValues = MapAlgebra.bandAsArray(clippedRaster, 1); + expectedValues = + new double[] { + 91.0, 92.0, 93.0, 94.0, 0.0, 96.0, 97.0, 98.0, 99.0, 100.0, 81.0, 82.0, 83.0, 84.0, 0.0, + 86.0, 87.0, 88.0, 89.0, 90.0, 71.0, 72.0, 73.0, 74.0, 75.0, 76.0, 77.0, 78.0, 79.0, 80.0, + 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0, 68.0, 0.0, 0.0, 51.0, 52.0, 53.0, 54.0, 55.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 42.0, 43.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE); + assertArrayEquals(expectedValues, actualValues, 0.0); } @Test @@ -287,6 +365,8 @@ public void testClip() ClassNotFoundException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); + String polygon = "POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170 4197590, 236722 4204770))"; Geometry geom = Constructors.geomFromWKT(polygon, RasterAccessors.srid(raster)); @@ -350,6 +430,53 @@ public void testClip() actualValues = PixelFunctions.values(croppedRaster, points, 1).toArray(new Double[0]); expectedValues = new Double[] {85.0, 85.0, 127.0, 212.0, null}; assertArrayEquals(expectedValues, actualValues); + + // Test for bottom-up raster + // Testing red band without crop + clippedRaster = RasterBandEditors.clip(raster_bottom_up, 1, geom, false, 200, false); + + clippedMetadata = RasterAccessors.metadata(clippedRaster); + originalMetadata = RasterAccessors.metadata(raster_bottom_up); + assertArrayEquals( + Arrays.stream(originalMetadata, 0, 9).toArray(), + Arrays.stream(clippedMetadata, 0, 9).toArray(), + 0.01d); + assertEquals(1, clippedMetadata[9], FP_TOLERANCE); + actual = String.valueOf(clippedRaster.getSampleDimensions()[0]); + expected = + "RenderedSampleDimension(\"RED_BAND\":[200.0 ... 200.0])\n ‣ Category(\"No data\":[200...200])\n"; + assertEquals(expected, actual); + + points = new ArrayList<>(); + points.add(Constructors.geomFromWKT("POINT(223802 4.21769e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(224759 4.20453e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237201 4.20429e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237919 4.20357e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(254668 4.21769e+06)", 26918)); + + actualValues = PixelFunctions.values(clippedRaster, points, 1).toArray(new Double[0]); + expectedValues = new Double[] {null, null, 0.0, 0.0, null}; + assertTrue(Arrays.equals(expectedValues, actualValues)); + + // Testing green band with crop + croppedRaster = RasterBandEditors.clip(raster_bottom_up, 2, geom, false, 200, true); + assertEquals(0, croppedRaster.getRenderedImage().getMinX()); + assertEquals(0, croppedRaster.getRenderedImage().getMinY()); + croppedRaster2 = Serde.deserialize(Serde.serialize(croppedRaster)); + assertSameCoverage(croppedRaster, croppedRaster2); + + points = new ArrayList<>(); + points.add(Constructors.geomFromWKT("POINT(236842 4.20465e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(236961 4.20453e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237201 4.20429e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(237919 4.20357e+06)", 26918)); + points.add(Constructors.geomFromWKT("POINT(223802 4.20465e+06)", 26918)); + + actual = String.valueOf(croppedRaster.getSampleDimensions()[0]); + expected = + "RenderedSampleDimension(\"GREEN_BAND\":[200.0 ... 200.0])\n" + + " ‣ Category(\"No data\":[200...200])\n"; + assertEquals(expected, actual); } @Test @@ -357,6 +484,7 @@ public void testClipLenient() throws FactoryException, IOException, ParseException, TransformException { GridCoverage2D raster = rasterFromGeoTiff(resourceFolder + "raster_geotiff_color/FAA_UTM18N_NAD83.tif"); + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); // Construct a polygon that does not intersect with the raster Geometry nonIntersectingGeom = @@ -376,6 +504,17 @@ public void testClipLenient() RasterBandEditors.clip(raster, 1, nonIntersectingGeom, false, 200, false); assertNull(result); raster.dispose(true); + + // Test for bottom-up raster + assertThrows( + IllegalArgumentException.class, + () -> + RasterBandEditors.clip( + raster_bottom_up, 1, nonIntersectingGeom, false, 200, false, false)); + + result = RasterBandEditors.clip(raster_bottom_up, 1, nonIntersectingGeom, false, 200, false); + assertNull(result); + raster_bottom_up.dispose(true); } @Test diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java index 7e8b7643b39..7bfb37c30dc 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterConstructorsTest.java @@ -18,6 +18,7 @@ */ package org.apache.sedona.common.raster; +import static org.apache.sedona.common.utils.RasterUtils.flipVerticallyPixelSpace; import static org.junit.Assert.assertArrayEquals; import static org.junit.Assert.assertEquals; @@ -81,6 +82,7 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio // Polygon GridCoverage2D raster = RasterConstructors.makeEmptyRaster(2, 255, 255, 1, -1, 2, -2, 0, 0, 4326); + GridCoverage2D raster_bottom_up = flipVerticallyPixelSpace(raster); Geometry geom = Constructors.geomFromWKT("POLYGON((15 -15, 18 -20, 15 -24, 24 -25, 15 -15))", 0); @@ -108,6 +110,22 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + expected = + new double[] { + 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + actual = MapAlgebra.bandAsArray(rasterized, 1); + assertArrayEquals(expected, actual, 0.1d); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + // MultiPolygon geom = Constructors.geomFromWKT( @@ -134,6 +152,23 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio }; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, + 0.0, 0.0 + }; + assertArrayEquals(expected, actual, 0.1d); + // MultiLineString geom = Constructors.geomFromWKT("MULTILINESTRING ((5 -5, 10 -10), (10 -10, 15 -15, 20 -20))", 0); @@ -151,6 +186,28 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio }; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 3d); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0 + }; + + assertArrayEquals(expected, actual, 0.1d); + rasterized = RasterConstructors.asRaster(geom, raster, "d"); actual = MapAlgebra.bandAsArray(rasterized, 1); @@ -165,6 +222,27 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio }; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + assertArrayEquals(expected, actual, 0.1d); + // LinearRing geom = Constructors.geomFromWKT("LINEARRING (10 -10, 18 -20, 15 -24, 24 -25, 10 -10)", 0); rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 3093151, 3d); @@ -181,6 +259,28 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio }; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 3d); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 3093151.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, + 3093151.0, 3093151.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, + 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0 + }; + + assertArrayEquals(expected, actual, 0.1d); + rasterized = RasterConstructors.asRaster(geom, raster, "d"); actual = MapAlgebra.bandAsArray(rasterized, 1); @@ -194,6 +294,25 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + + // Making sure orientation is preserved + assertEquals( + RasterAccessors.metadata(raster_bottom_up)[5], + RasterAccessors.metadata(rasterized)[5], + 0.0); + + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + assertArrayEquals(expected, actual, 0.1d); + // MultiPoints geom = Constructors.geomFromWKT("MULTIPOINT ((5 -5), (10 -10), (15 -15))", 0); rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 3093151, 3d); @@ -208,6 +327,19 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio }; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 3d); + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, + 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0, + 3093151.0, 3093151.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + + assertArrayEquals(expected, actual, 0.1d); + rasterized = RasterConstructors.asRaster(geom, raster, "d"); actual = MapAlgebra.bandAsArray(rasterized, 1); @@ -220,6 +352,17 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + actual = MapAlgebra.bandAsArray(rasterized, 1); + expected = + new double[] { + 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, + 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, + 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 + }; + assertArrayEquals(expected, actual, 0.1d); + // Point geom = Constructors.geomFromWKT("POINT (5 -5)", 0); rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 3093151, 3d); @@ -229,10 +372,20 @@ public void testAsRasterWithEmptyRaster() throws FactoryException, ParseExceptio expected = new double[] {3093151.0, 3093151.0, 3093151.0, 3093151.0}; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 3d); + actual = MapAlgebra.bandAsArray(rasterized, 1); + assertArrayEquals(expected, actual, 0.1d); + rasterized = RasterConstructors.asRaster(geom, raster, "d"); actual = MapAlgebra.bandAsArray(rasterized, 1); expected = new double[] {1.0, 1.0, 1.0, 1.0}; assertArrayEquals(expected, actual, 0.1d); + + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d"); + actual = MapAlgebra.bandAsArray(rasterized, 1); + assertArrayEquals(expected, actual, 0.1d); } @Test @@ -240,18 +393,31 @@ public void testAsRasterLingString() throws FactoryException, ParseException { // Horizontal LineString GridCoverage2D raster = RasterConstructors.makeEmptyRaster(2, 255, 255, 1, -1, 2, -2, 0, 0, 4326); + GridCoverage2D raster_bottom_up = + RasterConstructors.makeEmptyRaster(2, 255, 255, 1, -511, 2, 2, 0, 0, 4326); + Geometry geom = Constructors.geomFromEWKT("LINESTRING(1 -1, 2 -1, 10 -1)"); GridCoverage2D rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 3093151, 0d); double[] actual = MapAlgebra.bandAsArray(rasterized, 1); double[] expected = new double[] {3093151.0, 3093151.0, 3093151.0, 3093151.0, 3093151.0}; assertArrayEquals(expected, actual, 0.1d); + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 0d); + actual = MapAlgebra.bandAsArray(rasterized, 1); + assertArrayEquals(expected, actual, 0.1d); + // Vertical LineString geom = Constructors.geomFromEWKT("LINESTRING(1 -1, 1 -2, 1 -10)"); rasterized = RasterConstructors.asRaster(geom, raster, "d", false, 3093151, 0d); actual = MapAlgebra.bandAsArray(rasterized, 1); expected = new double[] {3093151.0, 3093151.0, 3093151.0, 3093151.0, 3093151.0}; assertArrayEquals(expected, actual, 0.1d); + + // Test bottom-up raster case + rasterized = RasterConstructors.asRaster(geom, raster_bottom_up, "d", false, 3093151, 0d); + actual = MapAlgebra.bandAsArray(rasterized, 1); + assertArrayEquals(expected, actual, 0.1d); } @Test @@ -320,6 +486,9 @@ public void testAsRasterWithRasterExtent() throws IOException, ParseException, F public void testAsRasterWithRasterExtent2() throws FactoryException, ParseException { GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 5, 5, 0, 0.5, 0.1, -0.1, 0, 0, 4326); + GridCoverage2D raster_bottom_up = + RasterConstructors.makeEmptyRaster(1, 5, 5, 0, 0.0, 0.1, 0.1, 0, 0, 4326); + Geometry geom = Constructors.geomFromWKT("POLYGON((0.1 0.1, 0.1 0.4, 0.4 0.4, 0.4 0.1, 0.1 0.1))", 0); GridCoverage2D rasterized = @@ -332,6 +501,18 @@ public void testAsRasterWithRasterExtent2() throws FactoryException, ParseExcept assertEquals(5, RasterAccessors.getHeight(rasterized)); double sum = Arrays.stream(MapAlgebra.bandAsArray(rasterized, 1)).sum(); assertEquals(900, sum, 1e-6); // Covers 3x3 grid + + // Test bottom-up raster case + rasterized = + RasterConstructors.asRasterWithRasterExtent(geom, raster_bottom_up, "d", false, 100d, 0d); + assertEquals(0, rasterized.getEnvelope2D().getMinX(), 1e-6); + assertEquals(0, rasterized.getEnvelope2D().getMinY(), 1e-6); + assertEquals(0.5, rasterized.getEnvelope2D().getWidth(), 1e-6); + assertEquals(0.5, rasterized.getEnvelope2D().getHeight(), 1e-6); + assertEquals(5, RasterAccessors.getWidth(rasterized)); + assertEquals(5, RasterAccessors.getHeight(rasterized)); + sum = Arrays.stream(MapAlgebra.bandAsArray(rasterized, 1)).sum(); + assertEquals(900, sum, 1e-6); // Covers 3x3 grid } @Test diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterizationTests.java b/common/src/test/java/org/apache/sedona/common/raster/RasterizationTests.java index 9d03220c6e7..8050541ee19 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterizationTests.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterizationTests.java @@ -21,7 +21,6 @@ import static org.apache.sedona.common.raster.RasterAccessors.metadata; import org.geotools.api.referencing.FactoryException; -import org.geotools.api.referencing.operation.TransformException; import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.geometry.jts.ReferencedEnvelope; import org.junit.Assert; @@ -34,8 +33,7 @@ public class RasterizationTests extends RasterTestBase { private final WKTReader wktReader = new WKTReader(); @Test - public void testRasterizeGeomExtent() - throws ParseException, FactoryException, TransformException { + public void testRasterizeGeomExtent() throws ParseException, FactoryException { GridCoverage2D testRaster = RasterConstructors.makeEmptyRaster(1, "F", 4, 4, 4, 4, 1); double[] metadata = metadata(testRaster); @@ -114,11 +112,167 @@ public void testRasterizeGeomExtent() validateRasterizeGeomExtent( wktLine, new double[] {5.0, 6.0, 1.0, 3.0}, testRaster, metadata, false); + // horizontal line at the top-left edge + wktLine = "LINESTRING (4.0 4.0, 5.5 4.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 6.0, 3.0, 4.0}, testRaster, metadata, false); + + // horizontal line at the top-right edge + wktLine = "LINESTRING (6.5 4.0, 8.0 4.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {6.0, 8.0, 3.0, 4.0}, testRaster, metadata, false); + + // horizontal line at the bottom edge + wktLine = "LINESTRING (4.0 0.0, 8.0 0.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 8.0, 0.0, 1.0}, testRaster, metadata, false); + + // vertical line + wktLine = "LINESTRING (5.0 2.0, 5.0 3.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 6.0, 2.0, 3.0}, testRaster, metadata, false); + + // vertical line at the right edge + wktLine = "LINESTRING (8.0 0.0, 8.0 3.5)"; + validateRasterizeGeomExtent( + wktLine, new double[] {7.0, 8.0, 0.0, 4.0}, testRaster, metadata, false); + + // vertical line at the left edge + wktLine = "LINESTRING (4.0 0.0, 4.0 3.5)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 5.0, 0.0, 4.0}, testRaster, metadata, false); + + // diagonal line + wktLine = "LINESTRING (5.0 2.0, 6.0 3.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {5.0, 6.0, 2.0, 3.0}, testRaster, metadata, false); + + // point tests + String wktPoint = "POINT (5.0 2.0)"; + validateRasterizeGeomExtent( + wktPoint, new double[] {4.0, 6.0, 1.0, 3.0}, testRaster, metadata, false); + + // intersecting 2 pixels + wktPoint = "POINT (5.0 2.5)"; + validateRasterizeGeomExtent( + wktPoint, new double[] {4.0, 6.0, 2.0, 3.0}, testRaster, metadata, false); + + // within pixel + wktPoint = "POINT (5.25 2.25)"; + validateRasterizeGeomExtent( + wktPoint, new double[] {5.0, 6.0, 2.0, 3.0}, testRaster, metadata, false); + } + + @Test + public void testRasterizeGeomExtentWithBottomUpRaster() throws ParseException, FactoryException { + GridCoverage2D testRaster = RasterConstructors.makeEmptyRaster(1, 4, 4, 4, 0, 1, 1, 0, 0, 4326); + double[] metadata = metadata(testRaster); + + // Grid aligned polygon completely contained within raster extent + String wktPolygon = "POLYGON ((5 2, 6 1, 5 3, 7 3, 5 2))"; + validateRasterizeGeomExtent( + wktPolygon, new double[] {5.0, 7.0, 1.0, 3.0}, testRaster, metadata, false); + + // Grid aligned polygon in line with raster extent + wktPolygon = "POLYGON ((4 2, 6 0, 8 2, 6 4, 4 2))"; + validateRasterizeGeomExtent( + wktPolygon, new double[] {4.0, 8.0, 0.0, 4.0}, testRaster, metadata, false); + + // Polygon partially outside raster extent + wktPolygon = "POLYGON ((3 1, 5 -1, 8 2, 6 4, 3 1))"; + validateRasterizeGeomExtent( + wktPolygon, new double[] {4.0, 8.0, 0.0, 4.0}, testRaster, metadata, false); + + // Polygon completely outside raster extent + wktPolygon = "POLYGON ((-1 -1, 0 -2, 1 -1, 0 0, -1 -1))"; + validateRasterizeGeomExtent(wktPolygon, null, testRaster, metadata, false); + + // Partial pixel alignment + String wktPolygon5 = "POLYGON ((5.5 2.5, 4.5 0.5, 6.5 1.5, 5.5 2.5))"; + validateRasterizeGeomExtent( + wktPolygon5, new double[] {4.0, 7.0, 0.0, 3.0}, testRaster, metadata, false); + + GridCoverage2D testRaster_frac = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, 4.3333, 3.3334, 0.3333, 0.3333, 0, 0, 4326); + double[] metadata_frac = metadata(testRaster_frac); + + // Grid aligned polygon completely contained within raster extent + wktPolygon = "Polygon((4.6666 4.0, 4.9999 3.7777, 5.3332 4.0, 4.9999 4.3333, 4.6666 4.0))"; + validateRasterizeGeomExtent( + wktPolygon, + new double[] {4.6666, 5.3332, 3.6667, 4.3333}, + testRaster_frac, + metadata_frac, + false); + + // Partial pixel alignment + wktPolygon = "Polygon((4.7 3.9, 4.9999 3.9, 5.1 4.0, 4.9999 4.2, 4.7 3.9))"; + validateRasterizeGeomExtent( + wktPolygon, + new double[] {4.6666, 5.3332, 3.6667, 4.3333}, + testRaster_frac, + metadata_frac, + false); + + // polygon larger than raster extent + wktPolygon = "Polygon((4.0 3.0, 4.0 2.0, 8.0 2.0, 8.0 5.0, 4.0 5.0, 4.0 3.0))"; + validateRasterizeGeomExtent( + wktPolygon, + new double[] {4.3333, 5.6665, 3.3334, 4.6666}, + testRaster_frac, + metadata_frac, + false); + + GridCoverage2D testRaster_neg = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, -4.3333, -5.9998, 0.3333, 0.3333, 0, 0, 4326); + double[] metadata_neg = metadata(testRaster_neg); + + // polygon larger than raster extent + wktPolygon = "Polygon((-8.0 -2.0, -2.0 -2.0, -2.0 -6.0, -8.0 -6.0, -8.0 -2.0))"; + validateRasterizeGeomExtent( + wktPolygon, + new double[] {-4.3333, -3.0001, -5.9998, -4.6666}, + testRaster_neg, + metadata_neg, + false); + + // horizontal line + String wktLine = "LINESTRING (5.0 2.0, 6.0 2.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {5.0, 6.0, 1.0, 3.0}, testRaster, metadata, false); + + // horizontal line at the upper left edge + wktLine = "LINESTRING (4.0 4.0, 5.5 4.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 6.0, 3.0, 4.0}, testRaster, metadata, false); + + // horizontal line at the top-right edge + wktLine = "LINESTRING (6.5 4.0, 8.0 4.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {6.0, 8.0, 3.0, 4.0}, testRaster, metadata, false); + + // horizontal line at the bottom edge + wktLine = "LINESTRING (4.0 0.0, 8.0 0.0)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 8.0, 0.0, 1.0}, testRaster, metadata, false); + // vertical line wktLine = "LINESTRING (5.0 2.0, 5.0 3.0)"; validateRasterizeGeomExtent( wktLine, new double[] {4.0, 6.0, 2.0, 3.0}, testRaster, metadata, false); + // vertical line at the right edge + wktLine = "LINESTRING (8.0 0.0, 8.0 3.5)"; + validateRasterizeGeomExtent( + wktLine, new double[] {7.0, 8.0, 0.0, 4.0}, testRaster, metadata, false); + + // vertical line at the left edge + wktLine = "LINESTRING (4.0 0.0, 4.0 3.5)"; + validateRasterizeGeomExtent( + wktLine, new double[] {4.0, 5.0, 0.0, 4.0}, testRaster, metadata, false); + // diagonal line wktLine = "LINESTRING (5.0 2.0, 6.0 3.0)"; validateRasterizeGeomExtent( diff --git a/common/src/test/java/org/apache/sedona/common/utils/RasterUtilsTest.java b/common/src/test/java/org/apache/sedona/common/utils/RasterUtilsTest.java index 72b8ffc2886..e34bbb08037 100644 --- a/common/src/test/java/org/apache/sedona/common/utils/RasterUtilsTest.java +++ b/common/src/test/java/org/apache/sedona/common/utils/RasterUtilsTest.java @@ -18,12 +18,28 @@ */ package org.apache.sedona.common.utils; +import static org.apache.sedona.common.Functions.setSRID; +import static org.apache.sedona.common.raster.GeometryFunctions.envelope; +import static org.apache.sedona.common.raster.PixelFunctions.value; +import static org.apache.sedona.common.raster.PixelFunctions.values; +import static org.geotools.filter.function.StaticGeometry.geomFromWKT; +import static org.junit.Assert.assertEquals; + import java.awt.Color; +import java.util.Arrays; +import java.util.List; +import org.apache.sedona.common.raster.MapAlgebra; +import org.apache.sedona.common.raster.RasterAccessors; +import org.apache.sedona.common.raster.RasterConstructors; +import org.geotools.api.referencing.FactoryException; +import org.geotools.api.referencing.operation.TransformException; import org.geotools.coverage.Category; import org.geotools.coverage.GridSampleDimension; +import org.geotools.coverage.grid.GridCoverage2D; import org.geotools.util.NumberRange; import org.junit.Assert; import org.junit.Test; +import org.locationtech.jts.geom.Geometry; public class RasterUtilsTest { @Test @@ -133,4 +149,117 @@ public void testNoDataValueUsingFloatQuantitativeCategory() { Assert.assertEquals("GrayScale", band3.getCategory(100).getName().toString()); Assert.assertEquals("GrayScale", band3.getCategory(199.999).getName().toString()); } + + @Test + public void testFlipVerticallyPixelSpace() throws FactoryException, TransformException { + GridCoverage2D raster_top_down = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, 401805.039562261, 2095852.150947876, 30, -30, 0, 0, 5070); + GridCoverage2D raster_bottom_up = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, 401805.039562261, 2095732.150947876, 30, 30, 0, 0, 5070); + + double[] bandValues_top_down = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}; + double[] bandValues_bottom_up = {13, 14, 15, 16, 9, 10, 11, 12, 5, 6, 7, 8, 1, 2, 3, 4}; + + raster_top_down = MapAlgebra.addBandFromArray(raster_top_down, bandValues_top_down, 1); + raster_bottom_up = MapAlgebra.addBandFromArray(raster_bottom_up, bandValues_bottom_up, 1); + + GridCoverage2D raster_flipped = RasterUtils.flipVerticallyPixelSpace(raster_bottom_up); + + double[] raster_bottom_up_metadata = RasterAccessors.metadata(raster_bottom_up); + double[] raster_flipped_metadata = RasterAccessors.metadata(raster_flipped); + + assertEquals(-raster_bottom_up_metadata[5], raster_flipped_metadata[5], 0); + assertEquals(envelope(raster_bottom_up), envelope(raster_flipped)); + + List testGeometries = + Arrays.asList( + setSRID(geomFromWKT("POINT (401820.039562261 2095837.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095837.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401820.039562261 2095747.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095747.150947876)"), 5070)); + + List values_bottom_up = values(raster_bottom_up, testGeometries); + List values_flipped = values(raster_flipped, testGeometries); + Assert.assertArrayEquals(values_bottom_up.toArray(), values_flipped.toArray()); + + raster_flipped = RasterUtils.flipVerticallyPixelSpace(raster_top_down); + + double[] raster_top_down_metadata = RasterAccessors.metadata(raster_top_down); + raster_flipped_metadata = RasterAccessors.metadata(raster_flipped); + + assertEquals(-raster_top_down_metadata[5], raster_flipped_metadata[5], 0); + assertEquals(envelope(raster_top_down), envelope(raster_flipped)); + + List values_top_down = values(raster_top_down, testGeometries); + values_flipped = values(raster_flipped, testGeometries); + Assert.assertArrayEquals(values_top_down.toArray(), values_flipped.toArray()); + } + + @Test + public void testFlipVerticallyGridSpace() throws FactoryException, TransformException { + GridCoverage2D raster_top_down = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, 401805.039562261, 2095852.150947876, 30, -30, 0, 0, 5070); + GridCoverage2D raster_bottom_up = + RasterConstructors.makeEmptyRaster( + 1, "F", 4, 4, 401805.039562261, 2095732.150947876, 30, 30, 0, 0, 5070); + + double[] bandValues_top_down = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16}; + double[] bandValues_bottom_up = {13, 14, 15, 16, 9, 10, 11, 12, 5, 6, 7, 8, 1, 2, 3, 4}; + + raster_top_down = MapAlgebra.addBandFromArray(raster_top_down, bandValues_top_down, 1); + raster_bottom_up = MapAlgebra.addBandFromArray(raster_bottom_up, bandValues_bottom_up, 1); + + GridCoverage2D raster_flipped = RasterUtils.flipVerticallyGridSpace(raster_bottom_up); + + double[] raster_bottom_up_metadata = RasterAccessors.metadata(raster_bottom_up); + double[] raster_flipped_metadata = RasterAccessors.metadata(raster_flipped); + + assertEquals(-raster_bottom_up_metadata[5], raster_flipped_metadata[5], 0); + assertEquals(envelope(raster_bottom_up), envelope(raster_flipped)); + + List testGeometries = + Arrays.asList( + setSRID(geomFromWKT("POINT (401820.039562261 2095837.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095837.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401820.039562261 2095747.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095747.150947876)"), 5070)); + + List testGeometriesFlipped = + Arrays.asList( + setSRID(geomFromWKT("POINT (401820.039562261 2095747.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095747.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401820.039562261 2095837.150947876)"), 5070), + setSRID(geomFromWKT("POINT (401910.039562261 2095837.150947876)"), 5070)); + + for (int i = 0; i < testGeometries.size(); i++) { + Double valueBottomUp = value(raster_bottom_up, testGeometries.get(i)); + Double valueFlipped = value(raster_flipped, testGeometriesFlipped.get(i)); + Assert.assertEquals( + "Pixel values should match at corresponding positions", + valueBottomUp, + valueFlipped, + 1e-6f); + } + + raster_flipped = RasterUtils.flipVerticallyGridSpace(raster_top_down); + + double[] raster_top_down_metadata = RasterAccessors.metadata(raster_top_down); + raster_flipped_metadata = RasterAccessors.metadata(raster_flipped); + + assertEquals(-raster_top_down_metadata[5], raster_flipped_metadata[5], 0); + assertEquals(envelope(raster_top_down), envelope(raster_flipped)); + + for (int i = 0; i < testGeometries.size(); i++) { + Double valueTopDown = value(raster_top_down, testGeometries.get(i)); + Double valueFlipped = value(raster_flipped, testGeometriesFlipped.get(i)); + Assert.assertEquals( + "Pixel values should match at corresponding positions", + valueTopDown, + valueFlipped, + 1e-6f); + } + } } diff --git a/docs/api/sql/Raster-Operators/RS_AsRaster.md b/docs/api/sql/Raster-Operators/RS_AsRaster.md index 7e5a8e8afa7..a339ba6fd42 100644 --- a/docs/api/sql/Raster-Operators/RS_AsRaster.md +++ b/docs/api/sql/Raster-Operators/RS_AsRaster.md @@ -59,7 +59,6 @@ Since: `v1.5.0` The function doesn't support rasters that have any one of the following properties: ``` ScaleX < 0 - ScaleY > 0 SkewX != 0 SkewY != 0 ``` diff --git a/spark/common/src/test/resources/rasterization/test_rasterization.tiff b/spark/common/src/test/resources/rasterization/test_rasterization.tiff index 0904a95ff2c951b890aa263897408274b6a891ba..e35d883be15831d51a1ea6b0a69b89ba60842c9e 100644 GIT binary patch literal 1218 zcmZwDO-~b17{&23(^@Jt1OkEqOd2(Y5JCV!P(%h)u*iFPDT!af4Pj}_jbO&OV9CNI zA#8O}$Wp6aQ7*^i#MTSsu&AkmM_zaV$0txU!`&OhCeKaPauC~i~2qHzs+{rck$OuI{x*8Z>U~s|Gbdg z?1{2dTclz*KbF0&X-~4hH_~-SRK7}Ws#4ZZ9^EZpreWTeDSi8Rc>3V{Ae@(6Y&Yij zmg?2%iBBAV^ZOL$_9-VB>cDQij{{i6hxizu;2;j+Fpl6Tj^Q{y#R;6mXE=q^_#9v0 z48Fu!e1&tEGoJqjs62Lyy{Hjt)Jl!5kfW)PC|0 DHV|r> literal 524672 zcmeI)xr!7~6b9g1j*5uL1PY=CA|fIp?)!3_#eHKbQ0PaJPGE zJ~%ZWF1zaIl=7uJ@@xLuotk=Pq0QsRlx}fqUMf4=d%M0{&Uw%H%zD0e-0=O^vTuB5 ze*F3N;bXl&&fVLwx2t~WIM%%TV{-fV9NxUY@Zsj??Qg5sr$3&Hg?IaRua{So*Z($u zeQ#Iyn5$pZl^r|VcHQ_(jT_E?pF9`&1yNrpCANRR@rnbPhcb_59?LwQd7`D0009C7 z2oRWkfjHTAj8mE0e~pb-oXI?!c`oyO=7r3QnU`8R2@oJafB=En7l_Mk$GDPtHS=2L z^~@WYH#2W#-p;(!(n){-0RjXF%)UU}Z97Ixvadtj%YHBWLiUC1i`f^mFSYI@K!5-N z0t9AXAeQ>NWGuH`V>$c%?Dw-j$o?Sv!|V^UKWg1cfB*pk1PIK&Ks@T}lJU6h8jrI- z$^In!)9g>PuVi1z{;YK;0RjXF5Fjx70`aV`OUCoIYdp`snte6u5UZ2>L+fibN~}fp zzTI#VAV7cs0Rl5G5Nku-F2MhnAZ)N^C~ zUPbmkb7=XxhJyeB0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF z5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk z1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs z0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZ zfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&U zAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C7 z2oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N z0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+ z009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBly zK!5-N0t5&UAV7cs0RjXF5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAV7cs0RjXF l5FkK+009C72oNAZfB*pk1PBlyK!5-N0t5&UAn?xxegbxdU3CBe