Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[SEDONA-515] Add handling of noDataValues in RS_Resample #1279

Merged
merged 5 commits into from
Mar 15, 2024
Merged

[SEDONA-515] Add handling of noDataValues in RS_Resample #1279

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
9d1d8df
test
prantogg Mar 14, 2024
2897e29
feat: Add handling for noDataValues
prantogg Mar 14, 2024
84c8529
Update imports
prantogg Mar 14, 2024
f78f9eb
Update imports
prantogg Mar 14, 2024
0dce197
Add tests; Update mean to median
prantogg Mar 15, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
25 changes: 24 additions & 1 deletion common/src/main/java/org/apache/sedona/common/raster/RasterEditors.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -186,18 +186,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
Open in desktop
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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,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 @@ -134,6 +134,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
Loading