From 98da1d73b650a50c17724b66ffb915821117ac09 Mon Sep 17 00:00:00 2001 From: Pranav Toggi Date: Fri, 15 Mar 2024 15:27:00 -0400 Subject: [PATCH] [SEDONA-515] Add handling of noDataValues in RS_Resample (#1279) * test * feat: Add handling for noDataValues * Update imports * Update imports * Add tests; Update mean to median --- .../sedona/common/raster/RasterEditors.java | 25 +++- .../sedona/common/utils/RasterUtils.java | 133 ++++++++++++++++++ .../common/raster/RasterEditorsTest.java | 82 ++++++++++- 3 files changed, 238 insertions(+), 2 deletions(-) diff --git a/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java index eb0db1eab4..c7a9e76f0d 100644 --- a/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java +++ b/common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java @@ -228,6 +228,10 @@ public static GridCoverage2D resample(GridCoverage2D raster, double widthOrScale PixelInCell.CELL_CORNER, transform, crs, null); Interpolation resamplingAlgorithm = Interpolation.getInstance(0); + GridCoverage2D newRaster; + GridCoverage2D noDataValueMask; + GridCoverage2D resampledNoDataValueMask; + if (!Objects.isNull(algorithm) && !algorithm.isEmpty()) { if (algorithm.equalsIgnoreCase("nearestneighbor")) { resamplingAlgorithm = Interpolation.getInstance(0); @@ -235,11 +239,30 @@ public static GridCoverage2D resample(GridCoverage2D raster, double widthOrScale resamplingAlgorithm = Interpolation.getInstance(1); }else if (algorithm.equalsIgnoreCase("bicubic")) { resamplingAlgorithm = Interpolation.getInstance(2); + } else { + throw new IllegalArgumentException("Invalid 'algorithm': '" + algorithm + "'. Expected one of: 'NearestNeighbor', 'Bilinear', 'Bicubic'."); } } - GridCoverage2D newRaster = (GridCoverage2D) Operations.DEFAULT.resample(raster, null, gridGeometry, resamplingAlgorithm); + + + if ((!Objects.isNull(algorithm) && !algorithm.isEmpty()) && (algorithm.equalsIgnoreCase("Bilinear") || algorithm.equalsIgnoreCase("Bicubic"))) { + // Create and resample noDataValue mask + noDataValueMask = RasterUtils.extractNoDataValueMask(raster); + resampledNoDataValueMask = resample(noDataValueMask, widthOrScale, heightOrScale, gridX, -gridY, useScale,"NearestNeighbor"); + + // Replace noDataValues with mean of neighbors and resample + raster = RasterUtils.replaceNoDataValues(raster); + newRaster = (GridCoverage2D) Operations.DEFAULT.resample(raster, null, gridGeometry, resamplingAlgorithm); + + // Apply resampled noDataValue mask to resampled raster + newRaster = RasterUtils.applyRasterMask(newRaster, resampledNoDataValueMask); + } else { + newRaster = (GridCoverage2D) Operations.DEFAULT.resample(raster, null, gridGeometry, resamplingAlgorithm); + } + return newRaster; } + public static GridCoverage2D resample(GridCoverage2D raster, double widthOrScale, double heightOrScale, boolean useScale, String algorithm) throws TransformException { return resample(raster, widthOrScale, heightOrScale, RasterAccessors.getUpperLeftX(raster), RasterAccessors.getUpperLeftY(raster), useScale, algorithm); 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 81a6f179cf..5f991d824c 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 @@ -20,6 +20,7 @@ import com.sun.media.imageioimpl.common.BogusColorSpace; import org.apache.commons.lang3.tuple.Pair; +import org.apache.commons.math3.stat.descriptive.rank.Median; import org.apache.sedona.common.Functions; import org.apache.sedona.common.FunctionsGeoTools; import org.apache.sedona.common.raster.RasterAccessors; @@ -676,4 +677,136 @@ public static GridCoverage2D shiftRasterToZeroOrigin(GridCoverage2D raster, Doub return RasterUtils.clone(raster.getRenderedImage(), raster.getGridGeometry(), raster.getSampleDimensions(), raster, noDataValue, true); } } + + /** + * Retrieves a List of neighboring pixel values excluding noDataValue neighbors + * + * @param x grid coordinate + * @param y grid coordinate + * @param band Band index of raster + * @param raster + * @param noDataValue + * @return + */ + public static List getNeighboringPixels(int x, int y, int band, Raster raster, Double noDataValue) { + List neighbors = new ArrayList<>(); + int width = raster.getWidth(); + int height = raster.getHeight(); + + for (int dx = -1; dx <= 1; dx++) { + for (int dy = -1; dy <= 1; dy++) { + int nx = x + dx; + int ny = y + dy; + + if (nx >= 0 && nx < width && ny >= 0 && ny < height && !(dx == 0 && dy == 0)) { + double value = raster.getSampleDouble(nx, ny, band); // Now checks the specified band + if (value != noDataValue) { + neighbors.add(value); + } + } + } + } + return neighbors; + } + + + /** + * Replaces noDataValue pixels in each band with mean of neighboring pixel values + * + * @param raster + * @return A new grid coverage with noDataValues pixels replaced + */ + public static GridCoverage2D replaceNoDataValues(GridCoverage2D raster) { + Raster rasterData = raster.getRenderedImage().getData(); + WritableRaster writableRaster = rasterData.createCompatibleWritableRaster(); + + Median medianCalculator = new Median(); + + // Iterate over each band + for (int band = 0; band < raster.getNumSampleDimensions(); band++) { + GridSampleDimension sampleDimension = raster.getSampleDimension(band); + double noDataValue = RasterUtils.getNoDataValue(sampleDimension); + + // Replace no data values with the median of neighboring pixels for each band + for (int y = 0; y < rasterData.getHeight(); y++) { + for (int x = 0; x < rasterData.getWidth(); x++) { + double originalValue = rasterData.getSampleDouble(x, y, band); + if (originalValue == noDataValue) { + List neighbors = RasterUtils.getNeighboringPixels(x, y, band, rasterData, noDataValue); + double[] neighborArray = neighbors.stream().mapToDouble(Double::doubleValue).toArray(); + double medianValue = neighborArray.length > 0 ? medianCalculator.evaluate(neighborArray) : Double.NaN; + writableRaster.setSample(x, y, band, !Double.isNaN(medianValue) ? medianValue : originalValue); + } else { + writableRaster.setSample(x, y, band, originalValue); + } + } + } + } + + GridCoverage2D modifiedRaster = RasterUtils.clone(writableRaster, raster.getGridGeometry(), raster.getSampleDimensions(), raster, null, true); + return modifiedRaster; + } + + /** + * Filters out the noDataValue pixels from each band as a grid coverage + * + * @param raster + * @return Returns a grid coverage with noDataValues and valid data values as Double.Nan + */ + public static GridCoverage2D extractNoDataValueMask(GridCoverage2D raster) { + Raster rasterData = raster.getRenderedImage().getData(); + WritableRaster writableRaster = rasterData.createCompatibleWritableRaster(RasterAccessors.getWidth(raster), RasterAccessors.getHeight(raster)); + + // Iterate over each band + for (int band = 0; band < raster.getNumSampleDimensions(); band++) { + GridSampleDimension sampleDimension = raster.getSampleDimension(band); + Double noDataValue = RasterUtils.getNoDataValue(sampleDimension); + + for (int y = 0; y < rasterData.getHeight(); y++) { + for (int x = 0; x < rasterData.getWidth(); x++) { + double originalValue = rasterData.getSampleDouble(x, y, band); + if (originalValue == noDataValue) { + writableRaster.setSample(x, y, band, originalValue); + } else { + writableRaster.setSample(x, y, band, Double.NaN); + } + } + } + } + + GridCoverage2D modifiedRaster = RasterUtils.clone(writableRaster, raster.getGridGeometry(), raster.getSampleDimensions(), raster, null, true); + return modifiedRaster; + } + + /** + * Superimposes the mask values onto the original raster, maintaining the original values where the mask is NaN. + * + * @param raster The original raster to which the mask will be applied. + * @param mask Grid coverage mask to be applied, containing the values to overlay. This mask should have the same dimensions and number of bands as the original raster. + * @return A new GridCoverage2D object with the mask applied. + */ + public static GridCoverage2D applyRasterMask(GridCoverage2D raster, GridCoverage2D mask) { + Raster rasterData = raster.getRenderedImage().getData(); + Raster maskData = mask.getRenderedImage().getData(); + WritableRaster writableRaster = rasterData.createCompatibleWritableRaster(RasterAccessors.getWidth(raster), RasterAccessors.getHeight(raster)); + + // Iterate over each band + for (int band = 0; band < raster.getNumSampleDimensions(); band++) { + for (int y = 0; y < rasterData.getHeight(); y++) { + for (int x = 0; x < rasterData.getWidth(); x++) { + double originalValue = rasterData.getSampleDouble(x, y, band); + double maskValue = maskData.getSampleDouble(x, y, band); + if (!Double.isNaN(maskValue)) { + writableRaster.setSample(x, y, band, maskValue); + } else { + writableRaster.setSample(x, y, band, originalValue); + } + } + } + } + + GridCoverage2D modifiedRaster = RasterUtils.clone(writableRaster, raster.getGridGeometry(), raster.getSampleDimensions(), raster, null, false); + return modifiedRaster; + } + } diff --git a/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java b/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java index 895127525f..5a59a24c04 100644 --- a/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java +++ b/common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java @@ -180,7 +180,7 @@ public void testResample() throws FactoryException, TransformException { //verify correct interpolation assertEquals(expectedRes, res); double[] metadata = RasterAccessors.metadata(newRaster); - double[] expectedMetadata = {-0.33333333333333326,0.19999999999999996,6,5,1.388888888888889,-1.24,0,0,0,1}; + double[] expectedMetadata = {-0.33333333333333326, 0.19999999999999996, 6, 5, 1.388888888888889, -1.24, 0, 0, 0, 1}; //verify correct raster geometry for (int i = 0; i < metadata.length; i++) { assertEquals(expectedMetadata[i], metadata[i], 1e-6); @@ -204,6 +204,86 @@ public void testResample() throws FactoryException, TransformException { } } + @Test + public void testResampleNoDataValueHandled() throws FactoryException, TransformException { + double[] values1 = {1,99,3,4,99,6,7,99,9,10,99,12}; + double[] values2 = {10,11,-9999,13,14,15,-9999,17,18,19,-9999,21}; + GridCoverage2D raster1 = RasterConstructors.makeEmptyRaster(1, "d", 4, 3, 0, 0, 2, -2, 0, 0, 0); + raster1 = MapAlgebra.addBandFromArray(raster1, values1, 1, 99.0); + raster1 = MapAlgebra.addBandFromArray(raster1, values2, 2, -9999.0); + + values1 = new double[] {1,2,3,4,5,6,7,8,9,10,99,12,13,14,15,16,17,18,19,20,21,22,23,24,25}; + values2 = new double[] {10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,-9999}; + GridCoverage2D raster2 = RasterConstructors.makeEmptyRaster(1, "d", 5, 5, 0, 0, 2, -2, 0, 0, 0); + raster2 = MapAlgebra.addBandFromArray(raster2, values1, 1, 99.0); + raster2 = MapAlgebra.addBandFromArray(raster2, values2, 2, -9999.0); + + GridCoverage2D newRaster = RasterEditors.resample(raster1, 6, 5, 1, -1, false, "bilinear"); + + String res1 = RasterOutputs.asMatrix(newRaster, 1); + String res2 = RasterOutputs.asMatrix(newRaster, 2); + String expectedRes1 = "|99.000000 99.000000 99.000000 99.000000 99.000000 99.000000|\n" + + "|99.000000 3.838750 99.000000 4.479375 4.400208 4.495000|\n" + + "|99.000000 99.000000 5.985764 6.593403 6.169791 99.000000|\n" + + "|99.000000 99.000000 8.250487 7.955348 8.473751 99.000000|\n" + + "|99.000000 9.375000 9.895833 99.000000 99.000000 12.000000|\n"; + String expectedRes2 = "|-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000|\n" + + "|-9999.000000 11.695000 12.482500 -9999.000000 -9999.000000 14.320000|\n" + + "|-9999.000000 14.175000 14.876389 -9999.000000 -9999.000000 16.800000|\n" + + "|-9999.000000 16.655001 17.270278 -9999.000000 -9999.000000 19.280001|\n" + + "|-9999.000000 18.375000 18.930555 -9999.000000 -9999.000000 21.000000|\n"; + + double[] metadata = RasterAccessors.metadata(newRaster); + double[] expectedMetadata = {-0.33333298563957214, 0.20000000298023224, 6.0, 5.0, 1.3888888309399288, -1.2400000005960465, 0.0, 0.0, 0.0, 2.0}; + assertEquals(expectedRes1, res1); + assertEquals(expectedRes2, res2); + for (int i = 0; i < metadata.length; i++) { + assertEquals(expectedMetadata[i], metadata[i], 1e-6); + } + + newRaster = RasterEditors.resample(raster1, 6, 5, 1, -1, false, "bicubic"); + + res1 = RasterOutputs.asMatrix(newRaster, 1); + res2 = RasterOutputs.asMatrix(newRaster, 2); + expectedRes1 = "|99.000000 99.000000 99.000000 99.000000 99.000000 99.000000|\n" + + "|99.000000 3.539351 99.000000 4.273053 4.102367 4.205373|\n" + + "|99.000000 99.000000 5.866029 6.737416 6.227383 99.000000|\n" + + "|99.000000 99.000000 8.386340 8.053590 8.636799 99.000000|\n" + + "|99.000000 9.564744 10.207731 99.000000 99.000000 12.610778|\n"; + expectedRes2 = "|-9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000 -9999.000000|\n" + + "|-9999.000000 11.252455 12.144701 -9999.000000 -9999.000000 13.989338|\n" + + "|-9999.000000 14.089166 14.862710 -9999.000000 -9999.000000 16.841017|\n" + + "|-9999.000000 16.901485 17.557348 -9999.000000 -9999.000000 19.668177|\n" + + "|-9999.000000 18.642647 19.225651 -9999.000000 -9999.000000 21.418527|\n"; + + metadata = RasterAccessors.metadata(newRaster); + expectedMetadata = new double[] {-0.33333298563957214, 0.20000000298023224, 6.0, 5.0, 1.3888888309399288, -1.2400000005960465, 0.0, 0.0, 0.0, 2.0}; + assertEquals(expectedRes1, res1); + assertEquals(expectedRes2, res2); + for (int i = 0; i < metadata.length; i++) { + assertEquals(expectedMetadata[i], metadata[i], 1e-6); + } + + newRaster = RasterEditors.resample(raster2, 3, 3, 1, -1, false, "bilinear"); + + res1 = RasterOutputs.asMatrix(newRaster, 1); + res2 = RasterOutputs.asMatrix(newRaster, 2); + expectedRes1 = "|99.000000 99.000000 99.000000|\n" + + "|99.000000 9.500002 11.555557|\n" + + "|99.000000 19.777779 21.833334|\n"; + expectedRes2 = "|-9999.000000 -9999.000000 -9999.000000|\n" + + "|-9999.000000 18.500002 20.555557|\n" + + "|-9999.000000 28.777779 29.718364|\n"; + + metadata = RasterAccessors.metadata(newRaster); + expectedMetadata = new double[] {-2.3333330154418945, 2.3333330154418945, 3.0, 3.0, 4.111111005147298, -4.111111005147298, 0.0, 0.0, 0.0, 2.0}; + assertEquals(expectedRes1, res1); + assertEquals(expectedRes2, res2); + for (int i = 0; i < metadata.length; i++) { + assertEquals(expectedMetadata[i], metadata[i], 1e-6); + } + } + @Test public void testResampleResizeFlavor() throws FactoryException, TransformException { double[] values = {1, 2, 3, 5, 4, 5, 6, 9, 7, 8, 9, 10};