This is an automated email from the ASF dual-hosted git repository.
jiayu pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sedona.git
The following commit(s) were added to refs/heads/master by this push:
new 5a15d06553 [SEDONA-746] Fix RS_Clip behavior (#2400)
5a15d06553 is described below
commit 5a15d065538cce8f59674cf2d463164727dce2a0
Author: Pranav Toggi <[email protected]>
AuthorDate: Tue Oct 14 23:53:27 2025 -0700
[SEDONA-746] Fix RS_Clip behavior (#2400)
Co-authored-by: Copilot <[email protected]>
---
.../sedona/common/raster/RasterBandEditors.java | 231 +++++++++++----------
.../apache/sedona/common/raster/RasterOutputs.java | 86 ++++++++
.../apache/sedona/common/raster/Rasterization.java | 155 +++++++++-----
.../sedona/common/raster/FunctionEditorsTest.java | 4 +-
.../common/raster/RasterBandAccessorsTest.java | 20 ++
.../common/raster/RasterBandEditorsTest.java | 104 +++++++++-
.../common/raster/RasterConstructorsTest.java | 2 +-
.../sedona/common/raster/RasterizationTests.java | 162 +++++++++++++++
spark/common/src/test/resources/allNoData.tiff | Bin 0 -> 131468 bytes
.../rasterization/test_rasterization.tiff | Bin 0 -> 524672 bytes
10 files changed, 591 insertions(+), 173 deletions(-)
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java
b/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java
index 51c906429f..82f57610d6 100644
---
a/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java
+++
b/common/src/main/java/org/apache/sedona/common/raster/RasterBandEditors.java
@@ -27,7 +27,6 @@ import javax.media.jai.RasterFactory;
import org.apache.commons.lang3.ArrayUtils;
import org.apache.commons.lang3.tuple.Pair;
import org.apache.sedona.common.utils.RasterUtils;
-import org.geotools.api.geometry.BoundingBox;
import org.geotools.api.parameter.ParameterValueGroup;
import org.geotools.api.referencing.FactoryException;
import org.geotools.api.referencing.operation.TransformException;
@@ -35,8 +34,6 @@ import org.geotools.coverage.GridSampleDimension;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.processing.CannotCropException;
import org.geotools.coverage.processing.operation.Crop;
-import org.geotools.geometry.jts.JTS;
-import org.geotools.geometry.jts.ReferencedEnvelope;
import org.locationtech.jts.geom.Geometry;
public class RasterBandEditors {
@@ -269,6 +266,46 @@ public class RasterBandEditors {
}
}
+ public static GridCoverage2D crop(
+ GridCoverage2D raster, double noDataValue, Geometry geomExtent, boolean
lenient) {
+ // Crop the raster
+ // this will shrink the extent of the raster to the geometry
+ Crop cropObject = new Crop();
+ ParameterValueGroup parameters = cropObject.getParameters();
+ parameters.parameter("Source").setValue(raster);
+ parameters.parameter(Crop.PARAMNAME_DEST_NODATA).setValue(new double[]
{noDataValue});
+ parameters.parameter(Crop.PARAMNAME_ROI).setValue(geomExtent);
+
+ // crop the raster to the geometry extent
+ try {
+ raster = (GridCoverage2D) cropObject.doOperation(parameters, null);
+ } catch (CannotCropException e) {
+ if (lenient) {
+ return null;
+ } else {
+ throw e;
+ }
+ }
+
+ RenderedImage image = raster.getRenderedImage();
+ int minX = image.getMinX();
+ int minY = image.getMinY();
+ if (minX != 0 || minY != 0) {
+ raster = RasterUtils.shiftRasterToZeroOrigin(raster, noDataValue);
+ } else {
+ raster =
+ RasterUtils.clone(
+ raster.getRenderedImage(),
+ raster.getGridGeometry(),
+ raster.getSampleDimensions(),
+ raster,
+ noDataValue,
+ true);
+ }
+
+ return raster;
+ }
+
/**
* Return a clipped raster with the specified ROI by the geometry
*
@@ -292,121 +329,97 @@ public class RasterBandEditors {
boolean lenient)
throws FactoryException, TransformException {
- // Selecting the band from original raster
- RasterUtils.ensureBand(raster, band);
- GridCoverage2D singleBandRaster = RasterBandAccessors.getBand(raster, new
int[] {band});
-
- Pair<GridCoverage2D, Geometry> pair =
- RasterUtils.setDefaultCRSAndTransform(singleBandRaster, geometry);
- singleBandRaster = pair.getLeft();
- geometry = pair.getRight();
-
- // Use rasterizeGeomExtent for AOI geometries smaller than a pixel
- double[] metadata = RasterAccessors.metadata(singleBandRaster);
- ReferencedEnvelope geomEnvelope =
- Rasterization.rasterizeGeomExtent(geometry, singleBandRaster,
metadata, allTouched);
- Geometry geomExtent = JTS.toGeometry((BoundingBox) geomEnvelope);
-
- // Crop the raster
- // this will shrink the extent of the raster to the geometry
- Crop cropObject = new Crop();
- ParameterValueGroup parameters = cropObject.getParameters();
- parameters.parameter("Source").setValue(singleBandRaster);
- parameters.parameter(Crop.PARAMNAME_DEST_NODATA).setValue(new double[]
{noDataValue});
- parameters.parameter(Crop.PARAMNAME_ROI).setValue(geomExtent);
-
- GridCoverage2D newRaster;
- try {
- newRaster = (GridCoverage2D) cropObject.doOperation(parameters, null);
- } catch (CannotCropException e) {
+ if (!RasterPredicates.rsIntersects(raster, geometry)) {
if (lenient) {
return null;
- } else {
- throw e;
}
+ throw new IllegalArgumentException("Geometry does not intersect
Raster.");
}
- if (!crop) {
- double[] metadataOriginal = RasterAccessors.metadata(raster);
- int widthOriginalRaster = (int) metadataOriginal[2],
- heightOriginalRaster = (int) metadataOriginal[3];
- Raster rasterData = RasterUtils.getRaster(raster.getRenderedImage());
+ // Selecting the band from original raster
+ RasterUtils.ensureBand(raster, band);
+ Pair<GridCoverage2D, Geometry> pair =
RasterUtils.setDefaultCRSAndTransform(raster, geometry);
+ geometry = pair.getRight();
- // create a new raster and set a default value that's the no-data value
- String bandType = RasterBandAccessors.getBandType(raster, 1);
- int dataTypeCode =
RasterUtils.getDataTypeCode(RasterBandAccessors.getBandType(raster, 1));
- boolean isDataTypeIntegral =
RasterUtils.isDataTypeIntegral(dataTypeCode);
- WritableRaster resultRaster =
- RasterFactory.createBandedRaster(
- dataTypeCode, widthOriginalRaster, heightOriginalRaster, 1,
null);
- int sizeOfArray = widthOriginalRaster * heightOriginalRaster;
- if (isDataTypeIntegral) {
- int[] array =
- ArrayUtils.toPrimitive(
- Collections.nCopies(sizeOfArray, (int) noDataValue)
- .toArray(new Integer[sizeOfArray]));
- resultRaster.setSamples(0, 0, widthOriginalRaster,
heightOriginalRaster, 0, array);
- } else {
- double[] array =
- ArrayUtils.toPrimitive(
- Collections.nCopies(sizeOfArray, noDataValue).toArray(new
Double[sizeOfArray]));
- resultRaster.setSamples(0, 0, widthOriginalRaster,
heightOriginalRaster, 0, array);
- }
+ double[] rasterMetadata = RasterAccessors.metadata(raster);
+ int rasterWidth = (int) rasterMetadata[2], rasterHeight = (int)
rasterMetadata[3];
+ Raster rasterData = RasterUtils.getRaster(raster.getRenderedImage());
+
+ // create a new raster and set a default value that's the no-data value
+ String bandType = RasterBandAccessors.getBandType(raster, band);
+ int dataTypeCode =
RasterUtils.getDataTypeCode(RasterBandAccessors.getBandType(raster, band));
+ boolean isDataTypeIntegral = RasterUtils.isDataTypeIntegral(dataTypeCode);
+ WritableRaster writableRaster =
+ RasterFactory.createBandedRaster(dataTypeCode, rasterWidth,
rasterHeight, 1, null);
+ int sizeOfArray = rasterWidth * rasterHeight;
+ if (isDataTypeIntegral) {
+ int[] array =
+ ArrayUtils.toPrimitive(
+ Collections.nCopies(sizeOfArray, (int) noDataValue)
+ .toArray(new Integer[sizeOfArray]));
+ writableRaster.setSamples(0, 0, rasterWidth, rasterHeight, 0, array);
+ } else {
+ double[] array =
+ ArrayUtils.toPrimitive(
+ Collections.nCopies(sizeOfArray, noDataValue).toArray(new
Double[sizeOfArray]));
+ writableRaster.setSamples(0, 0, rasterWidth, rasterHeight, 0, array);
+ }
- // rasterize the geometry to iterate over the clipped raster
- GridCoverage2D rasterized =
- RasterConstructors.asRaster(geometry, raster, bandType, allTouched,
150);
- Raster rasterizedData =
RasterUtils.getRaster(rasterized.getRenderedImage());
- double[] metadataRasterized = RasterAccessors.metadata(rasterized);
- int widthRasterized = (int) metadataRasterized[2],
- heightRasterized = (int) metadataRasterized[3];
-
- for (int j = 0; j < heightRasterized; j++) {
- for (int i = 0; i < widthRasterized; i++) {
- Point2D point = RasterUtils.getWorldCornerCoordinates(rasterized, i,
j);
- int[] rasterCoord =
- RasterUtils.getGridCoordinatesFromWorld(raster, point.getX(),
point.getY());
- int x = Math.abs(rasterCoord[0]), y = Math.abs(rasterCoord[1]);
-
- if (rasterizedData.getPixel(i, j, (int[]) null)[0] == 0) {
- continue;
- }
-
- if (isDataTypeIntegral) {
- int[] pixelValue = rasterData.getPixel(x, y, (int[]) null);
-
- resultRaster.setPixel(x, y, new int[] {pixelValue[band - 1]});
- } else {
- double[] pixelValue = rasterData.getPixel(x, y, (double[]) null);
-
- resultRaster.setPixel(x, y, new double[] {pixelValue[band - 1]});
- }
+ // rasterize the geometry to iterate over the clipped raster
+ GridCoverage2D maskRaster =
+ RasterConstructors.asRaster(geometry, raster, bandType, allTouched,
150);
+ Raster maskData = RasterUtils.getRaster(maskRaster.getRenderedImage());
+ double[] maskMetadata = RasterAccessors.metadata(maskRaster);
+ int maskWidth = (int) maskMetadata[2], maskHeight = (int) maskMetadata[3];
+
+ // Calculate offset
+ Point2D point = RasterUtils.getWorldCornerCoordinates(maskRaster, 1, 1);
+
+ // Add half the pixel length to account for floating-point precision
issues.
+ // This adjustment increases the margin of error because the upper-left
corner of a pixel
+ // is very close to neighboring pixels. Without this adjustment,
floating-point inaccuracies
+ // can cause the calculated world coordinates to incorrectly fall into an
adjacent pixel.
+ // For example, the upperLeftY of a maskRaster pixel might be
243924.000000000001,
+ // while the upperLeftY of the original raster pixel is
243924.000000000000.
+ int[] rasterCoord =
+ RasterUtils.getGridCoordinatesFromWorld(
+ raster, point.getX() + (rasterMetadata[4] / 2), point.getY() +
(rasterMetadata[5] / 2));
+ double offsetX = rasterCoord[0];
+ double offsetY = rasterCoord[1];
+
+ for (int j = 0; j < maskHeight; j++) {
+ for (int i = 0; i < maskWidth; i++) {
+ // Calculate mapped raster index
+ int x = (int) (i + offsetX), y = (int) (j + offsetY);
+
+ // Check if the pixel in the maskRaster data is valid
+ if (maskData.getPixel(i, j, (int[]) null)[0] == 0) {
+ continue;
+ }
+
+ if (isDataTypeIntegral) {
+ int[] pixelValue = rasterData.getPixel(x, y, (int[]) null);
+ writableRaster.setPixel(x, y, new int[] {pixelValue[band - 1]});
+ } else {
+ double[] pixelValue = rasterData.getPixel(x, y, (double[]) null);
+ writableRaster.setPixel(x, y, new double[] {pixelValue[band - 1]});
}
}
- newRaster =
- RasterUtils.clone(
- resultRaster,
- raster.getGridGeometry(),
- newRaster.getSampleDimensions(),
- newRaster,
- noDataValue,
- true);
- } else {
- RenderedImage image = newRaster.getRenderedImage();
- int minX = image.getMinX();
- int minY = image.getMinY();
- if (minX != 0 || minY != 0) {
- newRaster = RasterUtils.shiftRasterToZeroOrigin(newRaster,
noDataValue);
- } else {
- newRaster =
- RasterUtils.clone(
- newRaster.getRenderedImage(),
- newRaster.getGridGeometry(),
- newRaster.getSampleDimensions(),
- newRaster,
- noDataValue,
- true);
- }
+ }
+
+ // Not cropped but clipped
+ GridCoverage2D newRaster =
+ RasterUtils.clone(
+ writableRaster,
+ raster.getGridGeometry(),
+ new GridSampleDimension[] {raster.getSampleDimension(band - 1)},
+ raster,
+ noDataValue,
+ true);
+
+ if (crop) {
+ Geometry maskRasterExtent = GeometryFunctions.envelope(maskRaster);
+ return crop(newRaster, noDataValue, maskRasterExtent, lenient);
}
return newRaster;
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java
b/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java
index 77d115352f..c2dc22782f 100644
--- a/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java
+++ b/common/src/main/java/org/apache/sedona/common/raster/RasterOutputs.java
@@ -29,6 +29,7 @@ import java.io.File;
import java.io.FileOutputStream;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
+import java.text.DecimalFormat;
import java.util.Base64;
import javax.imageio.ImageIO;
import javax.imageio.ImageWriteParam;
@@ -37,6 +38,7 @@ import javax.media.jai.JAI;
import javax.media.jai.RenderedOp;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.api.coverage.grid.GridCoverageWriter;
+import org.geotools.api.metadata.spatial.PixelOrientation;
import org.geotools.api.parameter.GeneralParameterValue;
import org.geotools.api.parameter.ParameterValueGroup;
import org.geotools.coverage.grid.GridCoverage2D;
@@ -45,6 +47,7 @@ import org.geotools.gce.arcgrid.ArcGridWriteParams;
import org.geotools.gce.arcgrid.ArcGridWriter;
import org.geotools.gce.geotiff.GeoTiffWriteParams;
import org.geotools.gce.geotiff.GeoTiffWriter;
+import org.geotools.referencing.operation.transform.AffineTransform2D;
public class RasterOutputs {
public static byte[] asGeoTiff(
@@ -294,6 +297,89 @@ public class RasterOutputs {
return Base64.getEncoder().encodeToString(png);
}
+ public static String asMatrixPretty(GridCoverage2D raster, int band, int
precision) {
+ try {
+ RasterUtils.ensureBand(raster, band);
+ RenderedImage img = raster.getRenderedImage();
+ Raster r = RasterUtils.getRaster(img);
+
+ final int minX = img.getMinX();
+ final int minY = img.getMinY();
+ final int width = r.getWidth();
+ final int height = r.getHeight();
+ final int bandIdx = band - 1;
+
+ // Affine transform (upper-left corners)
+ AffineTransform2D at = RasterUtils.getAffineTransform(raster,
PixelOrientation.UPPER_LEFT);
+ final double a = at.getScaleX(), b = at.getShearX(), d = at.getShearY(),
e = at.getScaleY();
+ final double tx = at.getTranslateX(), ty = at.getTranslateY();
+
+ final boolean integral =
RasterUtils.isDataTypeIntegral(r.getDataBuffer().getDataType());
+ final String numPat = precision <= 0 ? "0" : "0." +
"#".repeat(precision);
+ final DecimalFormat df = new DecimalFormat(numPat);
+
+ String[][] cells = new String[height + 1][width + 1];
+ cells[0][0] = "y\\x";
+
+ // Header: use upper-left corner (no +0.5)
+ for (int x = 0; x < width; x++) {
+ double xc = minX + x;
+ double yc = minY;
+ double worldX = a * xc + b * yc + tx;
+ cells[0][x + 1] = "[" + (minX + x) + ", " + df.format(worldX) + "]";
+ }
+
+ // Rows
+ for (int y = 0; y < height; y++) {
+ double yc = minY + y;
+ double xc = minX;
+ double worldY = d * xc + e * yc + ty;
+ cells[y + 1][0] = "[" + (minY + y) + ", " + df.format(worldY) + "]";
+
+ if (integral) {
+ int[] vals = r.getSamples(0, y, width, 1, bandIdx, (int[]) null);
+ for (int x = 0; x < width; x++) cells[y + 1][x + 1] =
String.valueOf(vals[x]);
+ } else {
+ double[] vals = r.getSamples(0, y, width, 1, bandIdx, (double[])
null);
+ for (int x = 0; x < width; x++) cells[y + 1][x + 1] =
df.format(vals[x]);
+ }
+ }
+
+ // --- Pretty print ---
+ int rows = height + 1, cols = width + 1;
+ int[] colWidths = new int[cols];
+ for (int c = 0; c < cols; c++) {
+ int max = 0;
+ for (int rIdx = 0; rIdx < rows; rIdx++) max = Math.max(max,
cells[rIdx][c].length());
+ colWidths[c] = max;
+ }
+
+ StringBuilder out = new StringBuilder();
+ String sep = "+";
+ for (int c = 0; c < cols; c++) sep += "-".repeat(colWidths[c] + 2) + "+";
+ sep += "\n";
+
+ out.append(sep);
+ for (int rIdx = 0; rIdx < rows; rIdx++) {
+ out.append("| ");
+ for (int c = 0; c < cols; c++) {
+ String val = cells[rIdx][c];
+ boolean leftAlign = (c == 0 || rIdx == 0);
+ if (leftAlign) out.append(String.format("%-" + colWidths[c] + "s",
val));
+ else out.append(String.format("%" + colWidths[c] + "s", val));
+ out.append(" | ");
+ }
+ out.append("\n");
+ if (rIdx == 0) out.append(sep);
+ }
+ out.append(sep);
+ return out.toString();
+
+ } catch (Exception e) {
+ throw new RuntimeException("Failed to format raster matrix", e);
+ }
+ }
+
public static String asMatrix(GridCoverage2D raster, int band, int
postDecimalPrecision) {
RasterUtils.ensureBand(raster, band);
Raster rasterData = RasterUtils.getRaster(raster.getRenderedImage());
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 730272c6e0..f5f32aecff 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
@@ -87,6 +87,11 @@ public class Rasterization {
boolean allTouched)
throws FactoryException {
+ // For instances where sub geometry is completely outside raster
+ if (geomExtent == null) {
+ return;
+ }
+
switch (geom.getGeometryType()) {
case "GeometryCollection":
case "MultiPolygon":
@@ -303,9 +308,24 @@ public class Rasterization {
return x >= minX && x <= maxX && y >= minY && y <= maxY;
}
- protected static ReferencedEnvelope rasterizeGeomExtent(
+ static ReferencedEnvelope rasterizeGeomExtent(
Geometry geom, GridCoverage2D raster, double[] metadata, boolean
allTouched) {
+ validateRasterMetadata(metadata);
+
+ // Always enable allTouched for MultiLineString and MultiPoint.
+ //
+ // Rationale:
+ // - Points and lines cannot be rasterized using "pixel-center inside
geometry" logic
+ // because they have no area. For these geometry types, we must use "any
touch"
+ // semantics to mark pixels they intersect/touch.
+ // - Single Point/LineString cases are already safeguarded by the
zero-envelope
+ // expansion logic at lines 367:373 (ensuring a non-degenerate search
window).
+ // - For Polygons, whether we expand the snapped extent depends on the
caller’s
+ // allTouched setting (center-based vs any-intersection semantics).
+ //
+ // Enabling allTouched here guarantees correct handling for multi-part
point/line
+ // geometries whose parts may sit exactly on pixel boundaries.
if (Objects.equals(geom.getGeometryType(), "MultiLineString")) {
allTouched = true;
}
@@ -313,57 +333,45 @@ public class Rasterization {
allTouched = true;
}
- ReferencedEnvelope rasterExtent =
+ ReferencedEnvelope geomExtent =
new ReferencedEnvelope(geom.getEnvelopeInternal(),
raster.getCoordinateReferenceSystem2D());
// Using BigDecimal to avoid floating point errors
- BigDecimal upperLeftX = BigDecimal.valueOf(metadata[0]);
- BigDecimal upperLeftY = BigDecimal.valueOf(metadata[1]);
- BigDecimal scaleX = BigDecimal.valueOf(metadata[4]);
- BigDecimal scaleY = BigDecimal.valueOf(metadata[5]);
+ double upperLeftX = metadata[0];
+ double upperLeftY = metadata[1];
+ double scaleX = metadata[4];
+ double scaleY = metadata[5];
// Compute the aligned min/max values
double alignedMinX =
- (scaleX
- .multiply(
- BigDecimal.valueOf(
- Math.floor((rasterExtent.getMinX() + metadata[0]) /
metadata[4])))
- .subtract(upperLeftX))
- .doubleValue();
-
+ toWorldCoordinate(
+ toPixelIndex(geomExtent.getMinX(), scaleX, upperLeftX, true),
scaleX, upperLeftX);
double alignedMinY =
- (scaleY
- .multiply(
- BigDecimal.valueOf(
- Math.ceil((rasterExtent.getMinY() + metadata[1]) /
metadata[5])))
- .subtract(upperLeftY))
- .doubleValue();
-
+ toWorldCoordinate(
+ toPixelIndex(geomExtent.getMinY(), scaleY, upperLeftY, true),
scaleY, upperLeftY);
double alignedMaxX =
- (scaleX
- .multiply(
- BigDecimal.valueOf(
- Math.ceil((rasterExtent.getMaxX() + metadata[0]) /
metadata[4])))
- .subtract(upperLeftX))
- .doubleValue();
-
+ toWorldCoordinate(
+ toPixelIndex(geomExtent.getMaxX(), scaleX, upperLeftX, false),
scaleX, upperLeftX);
double alignedMaxY =
- (scaleY
- .multiply(
- BigDecimal.valueOf(
- Math.floor((rasterExtent.getMaxY() + metadata[1]) /
metadata[5])))
- .subtract(upperLeftY))
- .doubleValue();
-
- // For points and LineStrings at intersection of 2 or more pixels,
- // extend search grid by 1 pixel in each direction
+ toWorldCoordinate(
+ toPixelIndex(geomExtent.getMaxY(), scaleY, upperLeftY, false),
scaleY, upperLeftY);
+
+ // Ensure the snapped AOI window is at least one pixel wide/tall.
+ // After snapping the continuous bbox to pixel edges, a point or thin line
can
+ // collapse to min == max along an axis (i.e., a degenerate
0-width/0-height window).
+ // That would produce an empty search region and skip rasterization
entirely.
+ //
+ // We grow the aligned window by exactly one pixel on each side, using the
+ // pixel size (scaleX > 0, scaleY < 0 for north-up). This yields at least a
+ // 1-pixel-wide/1-pixel-tall search window and guarantees we visit the
pixel(s)
+ // the geometry lies in.
if (alignedMaxX == alignedMinX) {
- alignedMinX -= metadata[4];
- alignedMaxX += metadata[4];
+ alignedMinX -= scaleX;
+ alignedMaxX += scaleX;
}
if (alignedMaxY == alignedMinY) {
- alignedMinY += metadata[5];
- alignedMaxY -= metadata[5];
+ alignedMinY += scaleY;
+ alignedMaxY -= scaleY;
}
// Get the extent of the original raster
@@ -372,15 +380,37 @@ public class Rasterization {
double originalMaxX = raster.getEnvelope().getMaximum(0);
double originalMaxY = raster.getEnvelope().getMaximum(1);
- // Extend the aligned extent by 1 pixel if allTouched is true and if any
geometry extent meets
- // the aligned extent
+ // Quick bbox intersection test for when rasterizeGeomExtent gets called
independently
+ // return null if they do not intersect
+ if (alignedMinX >= originalMaxX
+ || alignedMaxX <= originalMinX
+ || alignedMinY >= originalMaxY
+ || alignedMaxY <= originalMinY) {
+ return null;
+ }
+
+ // Handle "allTouched" behavior when geometry edges line up exactly with
pixel boundaries.
+ //
+ // Normally, each pixel is considered only if its *center* falls inside
the geometry.
+ // But if an edge of the geometry sits exactly on a pixel boundary, the
geometry
+ // might "touch" neighboring pixels without covering their centers — and
those pixels
+ // would be skipped.
+ //
+ // When allTouched = true, we expand the aligned bounding box outward by
one pixel
+ // on any side where the geometry’s edge exactly matches a pixel boundary.
+ // This guarantees we include those neighboring pixels that the geometry
merely touches.
+ //
+ // We only expand sides that line up perfectly with grid lines (equal
coordinates).
+ // The scaleX / scaleY values (which already encode pixel size and
direction) ensure
+ // this expansion moves exactly one pixel outward in each direction.
if (allTouched) {
- alignedMinX -= (rasterExtent.getMinX() == alignedMinX) ? metadata[4] : 0;
- alignedMinY += (rasterExtent.getMinY() == alignedMinY) ? metadata[5] : 0;
- alignedMaxX += (rasterExtent.getMaxX() == alignedMaxX) ? metadata[4] : 0;
- alignedMaxY -= (rasterExtent.getMaxY() == alignedMaxY) ? metadata[5] : 0;
+ alignedMinX -= (geomExtent.getMinX() == alignedMinX) ? scaleX : 0;
+ alignedMinY += (geomExtent.getMinY() == alignedMinY) ? scaleY : 0;
+ alignedMaxX += (geomExtent.getMaxX() == alignedMaxX) ? scaleX : 0;
+ alignedMaxY -= (geomExtent.getMaxY() == alignedMaxY) ? scaleY : 0;
}
+ // Clamp the aligned extent to the original raster extent
alignedMinX = Math.max(alignedMinX, originalMinX);
alignedMinY = Math.max(alignedMinY, originalMinY);
alignedMaxX = Math.min(alignedMaxX, originalMaxX);
@@ -388,16 +418,38 @@ public class Rasterization {
// Create the aligned raster extent
ReferencedEnvelope alignedRasterExtent =
- ReferencedEnvelope.rect(
+ new ReferencedEnvelope(
alignedMinX,
+ alignedMaxX,
alignedMinY,
- alignedMaxX - alignedMinX,
- alignedMaxY - alignedMinY,
- rasterExtent.getCoordinateReferenceSystem());
+ alignedMaxY,
+ geomExtent.getCoordinateReferenceSystem());
return alignedRasterExtent;
}
+ private static double toPixelIndex(double coord, double scale, double
upperLeft, boolean isMin) {
+ BigDecimal rel =
BigDecimal.valueOf(coord).subtract(BigDecimal.valueOf(upperLeft));
+
+ BigDecimal px = rel.divide(BigDecimal.valueOf(scale), 16,
RoundingMode.FLOOR);
+ if (scale > 0) {
+ return isMin
+ ? px.setScale(0, RoundingMode.FLOOR).doubleValue()
+ : px.setScale(0, RoundingMode.CEILING).doubleValue();
+ } else {
+ return isMin
+ ? px.setScale(0, RoundingMode.CEILING).doubleValue()
+ : px.setScale(0, RoundingMode.FLOOR).doubleValue();
+ }
+ }
+
+ private static double toWorldCoordinate(double pixelIndex, double scale,
double upperLeft) {
+ return BigDecimal.valueOf(pixelIndex)
+ .multiply(BigDecimal.valueOf(scale))
+ .add(BigDecimal.valueOf(upperLeft))
+ .doubleValue();
+ }
+
private static RasterizationParams calculateRasterizationParams(
GridCoverage2D raster,
boolean useGeometryExtent,
@@ -574,7 +626,7 @@ public class Rasterization {
// calculating slope
for (double y = yStart; y >= yEnd; y--) {
double xIntercept = p1X; // Vertical line, xIntercept is constant
- if (xIntercept < 0 || xIntercept >=
params.writableRaster.getWidth()) {
+ if (xIntercept < 0 || xIntercept >
params.writableRaster.getWidth()) {
continue; // Skip xIntercepts outside geomExtent
}
scanlineIntersections.computeIfAbsent(y, k -> new
TreeSet<>()).add(xIntercept);
@@ -583,6 +635,7 @@ public class Rasterization {
double slope = (worldP2.y - worldP1.y) / (worldP2.x - worldP1.x);
double xMin = (geomExtent.getMinX() - params.upperLeftX) /
params.scaleX;
double xMax = (geomExtent.getMaxX() - params.upperLeftX) /
params.scaleX;
+
for (double y = yStart; y >= yEnd; y--) {
double xIntercept = p1X + ((p1Y - y) / slope);
if ((xIntercept < xMin) || (xIntercept >= xMax)) {
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java
b/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java
index 37c6d4ed73..c6310e8371 100644
---
a/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java
+++
b/common/src/test/java/org/apache/sedona/common/raster/FunctionEditorsTest.java
@@ -68,8 +68,8 @@ public class FunctionEditorsTest extends RasterTestBase {
expected =
new double[] {
0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0, 10.0,
10.0, 10.0, 10.0, 10.0,
- 10.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 10.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 10.0,
- 10.0, 0.0, 10.0, 0.0, 0.0, 0.0, 0.0, 10.0, 0.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
+ 10.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 0.0, 10.0,
+ 10.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 10.0, 0.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0, 10.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 10.0,
10.0, 10.0, 10.0, 0.0,
0.0, 0.0, 0.0, 10.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
10.0, 10.0, 10.0, 10.0,
10.0, 10.0, 10.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
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 d953d91ece..0eaea40341 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
@@ -182,6 +182,26 @@ public class RasterBandAccessorsTest extends
RasterTestBase {
raster, nonIntersectingGeom, 1, "sum", false, false, false));
}
+ @Test
+ public void testZonalStatsWithMultiPolygonAndAllNoDataValues()
+ throws IOException, ParseException, FactoryException {
+ GridCoverage2D raster = rasterFromGeoTiff(resourceFolder +
"allNoData.tiff");
+ String polygon =
+ "SRID=2263;MULTIPOLYGON (((994691.8668410989 222732.55147949007,
994838.9964562772 222657.33114599282, 994835.8129632992 222651.32007158513,
994866.5497303953 222635.3891649635, 994894.9626119346 222620.6621292782,
994897.9281735104 222627.20212832885, 994912.1027229384 222619.95519748307,
994956.740641462 222597.13404523468, 994906.6541878665 222501.29768674026,
994887.6889911637 222465.0093761085, 994837.7340744453 222492.35408672754,
994853.625535307 222529.50849369404, 994840 [...]
+ Geometry geom = Constructors.geomFromEWKT(polygon);
+
+ Double actual = RasterBandAccessors.getZonalStats(raster, geom, 1,
"count", false, false);
+ double expected = 20470.0;
+ assertEquals(expected, actual, 0d);
+
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "sum", false,
false);
+ expected = 0;
+ assertEquals(expected, actual, 0d);
+
+ actual = RasterBandAccessors.getZonalStats(raster, geom, 1, "mean", false,
true);
+ assertNull(actual);
+ }
+
@Test
public void testZonalStatsWithNoData() throws IOException, FactoryException,
ParseException {
GridCoverage2D raster =
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 c46dd601b0..9b637c54a4 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
@@ -30,11 +30,11 @@ import java.util.Map;
import java.util.function.Function;
import java.util.stream.Collectors;
import org.apache.sedona.common.Constructors;
+import org.apache.sedona.common.FunctionsGeoTools;
import org.apache.sedona.common.raster.serde.Serde;
import org.geotools.api.referencing.FactoryException;
import org.geotools.api.referencing.operation.TransformException;
import org.geotools.coverage.grid.GridCoverage2D;
-import org.geotools.coverage.processing.CannotCropException;
import org.junit.Test;
import org.locationtech.jts.geom.Geometry;
import org.locationtech.jts.io.ParseException;
@@ -217,12 +217,70 @@ public class RasterBandEditorsTest extends RasterTestBase
{
double bandNoDataValue =
RasterBandAccessors.getBandNoDataValue(clippedRaster);
double expectedBandNoDataValue = 0.0;
double[] actualValues = MapAlgebra.bandAsArray(clippedRaster, 1);
- double[] expectedValues = {2.0, 3.0, 6.0, 7.0};
+ double[] expectedValues = {0.0, 0.0, 0.0, 7.0};
+
+ assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE);
+ assertTrue(Arrays.equals(expectedValues, actualValues));
+
+ Geometry geomTransformed = FunctionsGeoTools.transform(geom, "EPSG:5070",
"EPSG:4326");
+
+ clippedRaster = RasterBandEditors.clip(raster, 1, geomTransformed, false,
0, true);
+ bandNoDataValue = RasterBandAccessors.getBandNoDataValue(clippedRaster);
+ expectedBandNoDataValue = 0.0;
+ actualValues = MapAlgebra.bandAsArray(clippedRaster, 1);
+ expectedValues = new double[] {0.0, 0.0, 0.0, 7.0};
assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE);
assertTrue(Arrays.equals(expectedValues, actualValues));
}
+ @Test
+ public void testClipEdgeCase()
+ throws FactoryException, TransformException, ParseException, IOException
{
+ GridCoverage2D raster =
+ rasterFromGeoTiff(resourceFolder +
"rasterization/test_rasterization.tiff");
+ String polygon =
+ "POLYGON ((-112.33284156960647 35.2847218363903, -111.83787460663726
32.58335532290123, -110.70619885787484 32.647569129204015, -110.29612639967809
30.671467842387234, -108.19258739561617 31.285215692431052, -107.96086869190255
34.366572252095835, -108.22326988847186 35.68843747308277, -110.52624509609558
35.31045505826484, -112.33284156960647 35.2847218363903))";
+ Geometry geom = Constructors.geomFromWKT(polygon,
RasterAccessors.srid(raster));
+ GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
true, 0, true);
+
+ double bandNoDataValue =
RasterBandAccessors.getBandNoDataValue(clippedRaster);
+ double expectedBandNoDataValue = 0.0;
+ double[] actualValues = MapAlgebra.bandAsArray(clippedRaster, 1);
+ double[] expectedValues = {
+ 12.0, 13.0, 14.0, 15.0, 16.0, 0.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0,
32.0, 33.0, 34.0, 35.0,
+ 36.0, 37.0, 0.0, 43.0, 44.0, 45.0, 46.0, 47.0, 0.0, 0.0, 54.0, 55.0,
56.0, 0.0, 0.0, 0.0,
+ 64.0, 65.0, 0.0, 0.0
+ };
+
+ assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE);
+ assertArrayEquals(expectedValues, actualValues, 0.0);
+ }
+
+ @Test
+ public void testClipWithClippedGeometryWithHole()
+ throws FactoryException, TransformException, ParseException, IOException
{
+ GridCoverage2D raster =
+ rasterFromGeoTiff(resourceFolder +
"rasterization/test_rasterization.tiff");
+ String polygon =
+ "POLYGON ((-116.974798 27.410563, -113.10645 25.244469, -110.249147
27.566499, -110.205189 28.420171, -111.392068 29.036741, -109.28206 29.420242,
-107.831429 28.49744, -108.446848 27.566499, -107.391844 25.720509, -102.995994
24.08636, -99.479314 27.722214, -112.183321 32.063743, -116.974798 27.410563))";
+ Geometry geom = Constructors.geomFromWKT(polygon, 4326);
+ GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
true, 0, true);
+
+ double bandNoDataValue =
RasterBandAccessors.getBandNoDataValue(clippedRaster);
+ double expectedBandNoDataValue = 0.0;
+ double[] actualValues = MapAlgebra.bandAsArray(clippedRaster, 1);
+ double[] expectedValues = {
+ 0.0, 42.0, 43.0, 0.0, 0.0, 0.0, 0.0, 0.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, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0, 68.0, 0.0, 0.0,
71.0, 72.0, 73.0,
+ 74.0, 75.0, 76.0, 77.0, 78.0, 79.0, 80.0, 81.0, 82.0, 83.0, 84.0, 0.0,
86.0, 87.0, 88.0, 89.0,
+ 90.0, 91.0, 92.0, 93.0, 94.0, 0.0, 96.0, 97.0, 98.0, 99.0, 100.0
+ };
+
+ assertEquals(expectedBandNoDataValue, bandNoDataValue, FP_TOLERANCE);
+ assertArrayEquals(expectedValues, actualValues, 0.0);
+ }
+
@Test
public void testClip()
throws IOException, FactoryException, TransformException, ParseException,
@@ -233,11 +291,17 @@ public class RasterBandEditorsTest extends RasterTestBase
{
"POLYGON ((236722 4204770, 243900 4204770, 243900 4197590, 221170
4197590, 236722 4204770))";
Geometry geom = Constructors.geomFromWKT(polygon,
RasterAccessors.srid(raster));
+ // Testing red band without crop
GridCoverage2D clippedRaster = RasterBandEditors.clip(raster, 1, geom,
false, 200, false);
- double[] clippedMetadata =
- Arrays.stream(RasterAccessors.metadata(clippedRaster), 0, 9).toArray();
- double[] originalMetadata =
Arrays.stream(RasterAccessors.metadata(raster), 0, 9).toArray();
- assertArrayEquals(originalMetadata, clippedMetadata, 0.01d);
+
+ double[] clippedMetadata = RasterAccessors.metadata(clippedRaster);
+ double[] originalMetadata = RasterAccessors.metadata(raster);
+ assertArrayEquals(
+ Arrays.stream(originalMetadata, 0, 9).toArray(),
+ Arrays.stream(clippedMetadata, 0, 9).toArray(),
+ 0.01d);
+
+ assertEquals(1, clippedMetadata[9], FP_TOLERANCE);
String actual = String.valueOf(clippedRaster.getSampleDimensions()[0]);
String expected =
@@ -254,7 +318,9 @@ public class RasterBandEditorsTest extends RasterTestBase {
Double[] expectedValues = new Double[] {null, null, 0.0, 0.0, null};
assertTrue(Arrays.equals(expectedValues, actualValues));
- GridCoverage2D croppedRaster = RasterBandEditors.clip(raster, 1, geom,
false, 200, true);
+ // Testing green band with crop
+ GridCoverage2D croppedRaster = RasterBandEditors.clip(raster, 2, geom,
false, 200, true);
+
assertEquals(0, croppedRaster.getRenderedImage().getMinX());
assertEquals(0, croppedRaster.getRenderedImage().getMinY());
GridCoverage2D croppedRaster2 =
Serde.deserialize(Serde.serialize(croppedRaster));
@@ -265,9 +331,25 @@ public class RasterBandEditorsTest extends RasterTestBase {
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);
+
+ double[] croppedMetadata = RasterAccessors.metadata(clippedRaster);
+ originalMetadata = RasterAccessors.metadata(raster);
+ assertArrayEquals(
+ Arrays.stream(originalMetadata, 0, 9).toArray(),
+ Arrays.stream(croppedMetadata, 0, 9).toArray(),
+ 0.01d);
+
+ assertEquals(1, croppedMetadata[9], FP_TOLERANCE);
+
actualValues = PixelFunctions.values(croppedRaster, points, 1).toArray(new
Double[0]);
- expectedValues = new Double[] {0.0, 0.0, 0.0, 0.0, 255.0};
- assertTrue(Arrays.equals(expectedValues, actualValues));
+ expectedValues = new Double[] {85.0, 85.0, 127.0, 212.0, null};
+ assertArrayEquals(expectedValues, actualValues);
}
@Test
@@ -283,8 +365,10 @@ public class RasterBandEditorsTest extends RasterTestBase {
0);
// Throws an exception in non-lenient mode
+ // Throws "Geometry does not intersect Raster." instead - emitted by
Rasterization.rasterize()
+ // Should Rasterization.rasterize() have a lenient mode as well?
assertThrows(
- CannotCropException.class,
+ IllegalArgumentException.class,
() -> RasterBandEditors.clip(raster, 1, nonIntersectingGeom, false,
200, false, false));
// Returns null in lenient mode
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 f42d4acc5e..a5cc02aada 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
@@ -612,7 +612,7 @@ public class RasterConstructorsTest extends RasterTestBase {
if (offsetX + x < width && offsetY + y < height) {
float[] expectedValues = gridCoverage2D.evaluate(gridCoord,
(float[]) null);
if (bandIndices == null) {
- Assert.assertArrayEquals(expectedValues, values, 1e-6f);
+ assertArrayEquals(expectedValues, values, 1e-6f);
} else {
Assert.assertEquals(bandIndices.length, values.length);
for (int j = 0; j < bandIndices.length; j++) {
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
new file mode 100644
index 0000000000..9d03220c6e
--- /dev/null
+++
b/common/src/test/java/org/apache/sedona/common/raster/RasterizationTests.java
@@ -0,0 +1,162 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one
+ * or more contributor license agreements. See the NOTICE file
+ * distributed with this work for additional information
+ * regarding copyright ownership. The ASF licenses this file
+ * to you under the Apache License, Version 2.0 (the
+ * "License"); you may not use this file except in compliance
+ * with the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing,
+ * software distributed under the License is distributed on an
+ * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+ * KIND, either express or implied. See the License for the
+ * specific language governing permissions and limitations
+ * under the License.
+ */
+package org.apache.sedona.common.raster;
+
+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;
+import org.junit.Test;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.io.ParseException;
+import org.locationtech.jts.io.WKTReader;
+
+public class RasterizationTests extends RasterTestBase {
+ private final WKTReader wktReader = new WKTReader();
+
+ @Test
+ public void testRasterizeGeomExtent()
+ throws ParseException, FactoryException, TransformException {
+ GridCoverage2D testRaster = RasterConstructors.makeEmptyRaster(1, "F", 4,
4, 4, 4, 1);
+ 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, 4.6666, 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, -4.6666, 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);
+
+ // 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);
+
+ // 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);
+ }
+
+ private void validateRasterizeGeomExtent(
+ String wkt,
+ double[] expectedEnvelope,
+ GridCoverage2D raster,
+ double[] metadata,
+ boolean allTouched)
+ throws ParseException {
+ Geometry geom = wktReader.read(wkt);
+ ReferencedEnvelope envelope =
+ Rasterization.rasterizeGeomExtent(geom, raster, metadata, allTouched);
+ if (expectedEnvelope == null) {
+ Assert.assertNull(envelope);
+ } else {
+ Assert.assertEquals(expectedEnvelope[0], envelope.getMinX(),
FP_TOLERANCE);
+ Assert.assertEquals(expectedEnvelope[1], envelope.getMaxX(),
FP_TOLERANCE);
+ Assert.assertEquals(expectedEnvelope[2], envelope.getMinY(),
FP_TOLERANCE);
+ Assert.assertEquals(expectedEnvelope[3], envelope.getMaxY(),
FP_TOLERANCE);
+ }
+ }
+}
diff --git a/spark/common/src/test/resources/allNoData.tiff
b/spark/common/src/test/resources/allNoData.tiff
new file mode 100644
index 0000000000..2f5623148f
Binary files /dev/null and b/spark/common/src/test/resources/allNoData.tiff
differ
diff --git
a/spark/common/src/test/resources/rasterization/test_rasterization.tiff
b/spark/common/src/test/resources/rasterization/test_rasterization.tiff
new file mode 100644
index 0000000000..0904a95ff2
Binary files /dev/null and
b/spark/common/src/test/resources/rasterization/test_rasterization.tiff differ
