Skip to content

Commit

Permalink
[SEDONA-515] Add handling of noDataValues in RS_Resample (#1279)
Browse files Browse the repository at this point in the history 
* test

* feat: Add handling for noDataValues

* Update imports

* Update imports

* Add tests; Update mean to median
  • Loading branch information
@prantogg
prantogg committed Mar 15, 2024
1 parent 25f607a commit 98da1d7
Show file tree
Hide file tree
Showing 3 changed files with 238 additions and 2 deletions.
25 changes: 24 additions & 1 deletion common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -228,18 +228,41 @@ 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);
}else if (algorithm.equalsIgnoreCase("bilinear")) {
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);

Expand Down
133 changes: 133 additions & 0 deletions common/src/main/java/org/apache/sedona/common/utils/RasterUtils.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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<Double> 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<Double> getNeighboringPixels(int x, int y, int band, Raster raster, Double noDataValue) {
List<Double> 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<Double> 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;
}

}
82 changes: 81 additions & 1 deletion common/src/test/java/org/apache/sedona/common/raster/RasterEditorsTest.java
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);
Expand All @@ -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};
Expand Down

0 comments on commit 98da1d7

Please sign in to comment.