You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ji...@apache.org on 2023/08/16 07:53:03 UTC
[sedona] branch master updated: [SEDONA-361] RS_MapAlgebra for performing map algebra operations using simple expressions (#970)
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 0794529eb [SEDONA-361] RS_MapAlgebra for performing map algebra operations using simple expressions (#970)
0794529eb is described below
commit 0794529eb052777fb15a0c81b5d5de2100aef6fa
Author: Kristin Cowalcijk <bo...@wherobots.com>
AuthorDate: Wed Aug 16 15:52:58 2023 +0800
[SEDONA-361] RS_MapAlgebra for performing map algebra operations using simple expressions (#970)
---
common/pom.xml | 17 ++++
.../common/raster/DeepCopiedRenderedImage.java | 46 +++++++++--
.../apache/sedona/common/raster/MapAlgebra.java | 62 +++++++++++++++
.../sedona/common/raster/MapAlgebraTest.java | 64 +++++++++++++++
docs/api/sql/Raster-map-algebra.md | 90 ++++++++++++++++++++++
docs/api/sql/Raster-operators.md | 35 +++++++++
mkdocs.yml | 1 +
pom.xml | 19 +++++
.../scala/org/apache/sedona/sql/UDF/Catalog.scala | 1 +
.../sedona_sql/expressions/raster/MapAlgebra.scala | 12 ++-
.../org/apache/sedona/sql/rasteralgebraTest.scala | 68 ++++++++++++++++
11 files changed, 407 insertions(+), 8 deletions(-)
diff --git a/common/pom.xml b/common/pom.xml
index 14b7a08f5..0df2c5566 100644
--- a/common/pom.xml
+++ b/common/pom.xml
@@ -81,6 +81,23 @@
<groupId>net.sf.geographiclib</groupId>
<artifactId>GeographicLib-Java</artifactId>
</dependency>
+ <dependency>
+ <groupId>it.geosolutions.jaiext.jiffle</groupId>
+ <artifactId>jt-jiffle-language</artifactId>
+ </dependency>
+ <!-- These test dependencies are for running map algebra tests -->
+ <dependency>
+ <groupId>org.antlr</groupId>
+ <artifactId>antlr4-runtime</artifactId>
+ <version>${antlr-runtime.version}</version>
+ <scope>test</scope>
+ </dependency>
+ <dependency>
+ <groupId>org.codehaus.janino</groupId>
+ <artifactId>janino</artifactId>
+ <version>${janino-version}</version>
+ <scope>test</scope>
+ </dependency>
</dependencies>
<build>
<sourceDirectory>src/main/java</sourceDirectory>
diff --git a/common/src/main/java/org/apache/sedona/common/raster/DeepCopiedRenderedImage.java b/common/src/main/java/org/apache/sedona/common/raster/DeepCopiedRenderedImage.java
index e4ba21412..2774afd48 100644
--- a/common/src/main/java/org/apache/sedona/common/raster/DeepCopiedRenderedImage.java
+++ b/common/src/main/java/org/apache/sedona/common/raster/DeepCopiedRenderedImage.java
@@ -14,6 +14,7 @@
package org.apache.sedona.common.raster;
import com.sun.media.jai.util.ImageUtil;
+import it.geosolutions.jaiext.range.NoDataContainer;
import javax.media.jai.JAI;
import javax.media.jai.PlanarImage;
@@ -334,17 +335,23 @@ public final class DeepCopiedRenderedImage implements RenderedImage, Serializabl
boolean propertiesCloned = false;
Enumeration<String> keys = propertyTable.keys();
while (keys.hasMoreElements()) {
- Object key = keys.nextElement();
- if (!(this.properties.get(key) instanceof Serializable)) {
+ String key = keys.nextElement();
+ Object value = this.properties.get(key);
+ if (!(value instanceof Serializable)) {
if (!propertiesCloned) {
propertyTable = (Hashtable<String, Object>) this.properties.clone();
propertiesCloned = true;
}
- propertyTable.remove(key);
+ // GC_NODATA is a special property used by GeoTools. We need to serialize it.
+ if (value instanceof NoDataContainer) {
+ NoDataContainer noDataContainer = (NoDataContainer) value;
+ propertyTable.put(key, new SingleValueNoDataContainer(noDataContainer.getAsSingleValue()));
+ } else {
+ propertyTable.remove(key);
+ }
}
}
- out.writeObject(SerializerFactory.getState(this.sampleModel, null));
out.writeObject(SerializerFactory.getState(this.colorModel, null));
out.writeObject(propertyTable);
if (this.source != null) {
@@ -369,15 +376,40 @@ public final class DeepCopiedRenderedImage implements RenderedImage, Serializabl
@SuppressWarnings("unchecked")
private void readObject(ObjectInputStream in) throws IOException, ClassNotFoundException {
this.source = null;
- this.colorModel = null;
in.defaultReadObject();
- SerializableState smState = (SerializableState)in.readObject();
- this.sampleModel = (SampleModel)smState.getObject();
SerializableState cmState = (SerializableState)in.readObject();
this.colorModel = (ColorModel)cmState.getObject();
this.properties = (Hashtable<String, Object>)in.readObject();
+ for (String key : this.properties.keySet()) {
+ Object value = this.properties.get(key);
+ // Restore the value of GC_NODATA property as a NoDataContainer object.
+ if (value instanceof SingleValueNoDataContainer) {
+ SingleValueNoDataContainer noDataContainer = (SingleValueNoDataContainer) value;
+ this.properties.put(key, new NoDataContainer(noDataContainer.singleValue));
+ }
+ }
SerializableState rasState = (SerializableState)in.readObject();
this.imageRaster = (Raster)rasState.getObject();
+
+ // The deserialized rendered image contains only one tile (imageRaster). We need to update
+ // the sample model and tile properties to reflect this.
+ this.sampleModel = this.imageRaster.getSampleModel();
+ this.tileWidth = this.width;
+ this.tileHeight = this.height;
+ this.numXTiles = 1;
+ this.numYTiles = 1;
+ }
+
+ /**
+ * This class is for serializing NoDataContainer objects. It is mainly used to serialize the value of
+ * GC_NODATA property. We only considered the case where the NoDataContainer object contains only one
+ * value. However, it will cover most of the real world cases.
+ */
+ private static class SingleValueNoDataContainer implements Serializable {
+ private final double singleValue;
+ SingleValueNoDataContainer(double value) {
+ singleValue = value;
+ }
}
}
diff --git a/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java b/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java
index b1c06b152..0721df7a8 100644
--- a/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java
+++ b/common/src/main/java/org/apache/sedona/common/raster/MapAlgebra.java
@@ -18,11 +18,15 @@
*/
package org.apache.sedona.common.raster;
+import it.geosolutions.jaiext.jiffle.JiffleBuilder;
+import it.geosolutions.jaiext.jiffle.runtime.JiffleDirectRuntime;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.GridSampleDimension;
import org.geotools.coverage.grid.GridCoverage2D;
+import org.jaitools.imageutils.ImageUtils;
import javax.media.jai.RasterFactory;
+import javax.media.jai.TiledImage;
import java.awt.*;
import java.awt.image.Raster;
import java.awt.image.RenderedImage;
@@ -172,4 +176,62 @@ public class MapAlgebra
private static GridCoverage2D copyRasterAndReplaceBand(GridCoverage2D gridCoverage2D, int bandIndex, double[] bandValues) {
return copyRasterAndReplaceBand(gridCoverage2D, bandIndex, bandValues, null, false);
}
+
+ private static final ThreadLocal<String> previousScript = new ThreadLocal<>();
+ private static final ThreadLocal<JiffleDirectRuntime> previousRuntime = new ThreadLocal<>();
+
+ /**
+ * Applies a map algebra script to the given raster.
+ * @param gridCoverage2D The raster to apply the script to
+ * @param pixelType The pixel type of the output raster. If null, the pixel type of the input raster is used.
+ * @param script The script to apply
+ * @param noDataValue The no data value of the output raster.
+ * @return The result of the map algebra script
+ */
+ public static GridCoverage2D mapAlgebra(GridCoverage2D gridCoverage2D, String pixelType, String script, Double noDataValue) {
+ if (gridCoverage2D == null || script == null) {
+ return null;
+ }
+ RenderedImage renderedImage = gridCoverage2D.getRenderedImage();
+ int rasterDataType = pixelType != null? RasterUtils.getDataTypeCode(pixelType) : renderedImage.getSampleModel().getDataType();
+ int width = renderedImage.getWidth();
+ int height = renderedImage.getHeight();
+ TiledImage resultImage = ImageUtils.createConstantImage(width, height, 0.0);
+ try {
+ String prevScript = previousScript.get();
+ JiffleDirectRuntime prevRuntime = previousRuntime.get();
+ JiffleDirectRuntime runtime;
+ if (prevRuntime != null && script.equals(prevScript)) {
+ // Reuse the runtime to avoid recompiling the script
+ runtime = prevRuntime;
+ runtime.setSourceImage("rast", renderedImage);
+ runtime.setDestinationImage("out", resultImage);
+ runtime.setDefaultBounds();
+ } else {
+ JiffleBuilder builder = new JiffleBuilder();
+ runtime = builder.script(script).source("rast", renderedImage).dest("out", resultImage).getRuntime();
+ previousScript.set(script);
+ previousRuntime.set(runtime);
+ }
+
+ runtime.evaluateAll(null);
+
+ // If pixelType does not match with the data type of the result image (should be double since Jiffle only supports
+ // double destination image), we need to convert the resultImage to the specified pixel type.
+ if (rasterDataType != resultImage.getSampleModel().getDataType()) {
+ // Copy the resultImage to a new raster with the specified pixel type
+ WritableRaster convertedRaster = RasterFactory.createBandedRaster(rasterDataType, width, height, 1, null);
+ double[] samples = resultImage.getData().getSamples(0, 0, width, height, 0, (double[]) null);
+ convertedRaster.setSamples(0, 0, width, height, 0, samples);
+ resultImage.dispose();
+ return RasterUtils.create(convertedRaster, gridCoverage2D.getGridGeometry(), null, noDataValue);
+ } else {
+ // build a new GridCoverage2D from the resultImage
+ return RasterUtils.create(resultImage, gridCoverage2D.getGridGeometry(), null, noDataValue);
+ }
+ } catch (Exception e) {
+ resultImage.dispose();
+ throw new RuntimeException("Failed to run map algebra", e);
+ }
+ }
}
diff --git a/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java b/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java
index 4e30f4602..dc6c2c55a 100644
--- a/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java
+++ b/common/src/test/java/org/apache/sedona/common/raster/MapAlgebraTest.java
@@ -20,9 +20,13 @@ package org.apache.sedona.common.raster;
import org.apache.sedona.common.utils.RasterUtils;
import org.geotools.coverage.grid.GridCoverage2D;
+import org.junit.Assert;
import org.junit.Test;
import org.opengis.referencing.FactoryException;
+import java.awt.image.DataBuffer;
+import java.util.Random;
+
import static org.junit.Assert.assertEquals;
import static org.junit.Assert.assertNull;
import static org.junit.Assert.assertTrue;
@@ -190,4 +194,64 @@ public class MapAlgebraTest extends RasterTestBase
assertEquals(band[i], bandNew[i], 1e-9);
}
}
+
+ @Test
+ public void testMapAlgebra() throws FactoryException {
+ Random random = new Random();
+ String[] pixelTypes = {null, "b", "i", "s", "us", "f", "d"};
+ for (String pixelType : pixelTypes) {
+ int width = random.nextInt(100) + 10;
+ int height = random.nextInt(100) + 10;
+ testMapAlgebra(width, height, pixelType, null);
+ testMapAlgebra(width, height, pixelType, 100.0);
+ }
+ }
+
+ private void testMapAlgebra(int width, int height, String pixelType, Double noDataValue) throws FactoryException {
+ GridCoverage2D raster = RasterConstructors.makeEmptyRaster(2, "b", width, height, 10, 20, 1);
+ double[] band1 = new double[width * height];
+ double[] band2 = new double[width * height];
+ for (int i = 0; i < band1.length; i++) {
+ band1[i] = Math.random() * 10;
+ band2[i] = Math.random() * 10;
+ }
+ raster = MapAlgebra.addBandFromArray(raster, band1, 1);
+ raster = MapAlgebra.addBandFromArray(raster, band2, 2);
+ GridCoverage2D result = MapAlgebra.mapAlgebra(raster, pixelType, "out = (rast[0] + rast[1]) * 0.4;", noDataValue);
+ double actualNoDataValue = RasterUtils.getNoDataValue(result.getSampleDimension(0));
+ if (noDataValue != null) {
+ Assert.assertEquals(noDataValue, actualNoDataValue, 1e-9);
+ } else {
+ Assert.assertTrue(Double.isNaN(actualNoDataValue));
+ }
+
+ int resultDataType = result.getRenderedImage().getSampleModel().getDataType();
+ int expectedDataType;
+ if (pixelType != null) {
+ expectedDataType = RasterUtils.getDataTypeCode(pixelType);
+ } else {
+ expectedDataType = raster.getRenderedImage().getSampleModel().getDataType();
+ }
+ Assert.assertEquals(expectedDataType, resultDataType);
+
+ Assert.assertEquals(raster.getGridGeometry().getGridToCRS2D(), result.getGridGeometry().getGridToCRS2D());
+ band1 = MapAlgebra.bandAsArray(raster, 1);
+ band2 = MapAlgebra.bandAsArray(raster, 2);
+ double[] bandResult = MapAlgebra.bandAsArray(result, 1);
+ Assert.assertEquals(band1.length, bandResult.length);
+ for (int i = 0; i < band1.length; i++) {
+ double expected = (band1[i] + band2[i]) * 0.4;
+ double actual = bandResult[i];
+ switch (resultDataType) {
+ case DataBuffer.TYPE_BYTE:
+ case DataBuffer.TYPE_SHORT:
+ case DataBuffer.TYPE_USHORT:
+ case DataBuffer.TYPE_INT:
+ Assert.assertEquals((int) expected, (int) actual);
+ break;
+ default:
+ Assert.assertEquals(expected, actual, 1e-3);
+ }
+ }
+ }
}
diff --git a/docs/api/sql/Raster-map-algebra.md b/docs/api/sql/Raster-map-algebra.md
new file mode 100644
index 000000000..4e0aa9ffa
--- /dev/null
+++ b/docs/api/sql/Raster-map-algebra.md
@@ -0,0 +1,90 @@
+## Map Algebra
+
+Map algebra is a way to perform raster calculations using mathematical expressions. The expression can be a simple arithmetic operation or a complex combination of multiple operations. The expression can be applied to a single raster band or multiple raster bands. The result of the expression is a new raster.
+
+Apache Sedona provides two ways to perform map algebra operations:
+
+1. Using the `RS_MapAlgebra` function.
+2. Using `RS_BandAsArray` and array based map algebra functions, such as `RS_Add`, `RS_Multiply`, etc.
+
+Generally, the `RS_MapAlgebra` function is more flexible and can be used to perform more complex operations. The `RS_MapAlgebra(rast, pixelType, script, [noDataValue])` function takes three to four arguments:
+
+* `rast`: The raster to apply the map algebra expression to.
+* `pixelType`: The data type of the output raster. This can be one of `D` (double), `F` (float), `I` (integer), `S` (short), `US` (unsigned short) or `B` (byte). If specified `NULL`, the output raster will have the same data type as the input raster.
+* `script`: The map algebra script.
+* `noDataValue`: (Optional) The nodata value of the output raster.
+
+`RS_MapAlgebra` also has good performance, since it is backed by [Jiffle](https://github.com/geosolutions-it/jai-ext/wiki/Jiffle) and can be compiled to Java bytecode for
+execution. We'll demonstrate both approaches to implementing commonly used map algebra operations.
+
+### NDVI
+
+The Normalized Difference Vegetation Index (NDVI) is a simple graphical indicator that can be used to analyze remote sensing measurements, typically, but not necessarily, from a space platform, and assess whether the target being observed contains live green vegetation or not. NDVI has become a de facto standard index used to determine whether a given area contains live green vegetation or not. The NDVI is calculated from these individual measurements as follows:
+
+```
+NDVI = (NIR - Red) / (NIR + Red)
+```
+
+where NIR is the near-infrared band and Red is the red band.
+
+Assume that we have a bunch of rasters with 4 bands: red, green, blue, and near-infrared. We want to calculate the NDVI for each raster. We can use the `RS_MapAlgebra` function to do this:
+
+```sql
+SELECT RS_MapAlgebra(rast, 'D', 'out = (rast[3] - rast[0]) / (rast[3] + rast[0]);') as ndvi FROM raster_table
+```
+
+The Jiffle script is `out = (rast[3] - rast[0]) / (rast[3] + rast[0]);`. The `rast` variable is always bound to the input raster, and
+the `out` variable is bound to the output raster. Jiffle iterates over all the pixels in the input raster and executes the script for each pixel. the `rast[3]` and `rast[0]`
+refers to the current pixel values of the near-infrared and red bands, respectively. The `out` variable is the current output pixel value.
+
+The result of the `RS_MapAlgebra` function is a raster with a single band. The band is of type double, since we specified `D` as the `pixelType` argument.
+
+We can implement the same NDVI calculation using the array based map algebra functions:
+
+```sql
+SELECT RS_Divide(
+ RS_Subtract(RS_BandAsArray(rast, 1), RS_BandAsArray(rast, 4)),
+ RS_Add(RS_BandAsArray(rast, 1), RS_BandAsArray(rast, 4))) as ndvi FROM raster_table
+```
+
+The `RS_BandAsArray` function extracts the specified band of the input raster to an array of double, and the `RS_Add`, `RS_Subtract`, and `RS_Divide` functions perform the arithmetic operations on the arrays. The code using the array based map algebra functions is more verbose. However, there is a `RS_NormalizedDifference` function that can be used to calculate the NDVI more concisely:
+
+```sql
+SELECT RS_NormalizedDifference(RS_BandAsArray(rast, 1), RS_BandAsArray(rast, 4)) as ndvi FROM raster_table
+```
+
+The result of array based map algebra functions is an array of double. User can use `RS_AddBandFromArray` to add the array to a raster as a new band.
+
+### AWEI
+
+The Automated Water Extraction Index (AWEI) is a spectral index that can be used to extract water bodies from remote sensing imagery. The AWEI is calculated from these individual measurements as follows:
+
+```
+AWEI = 4 * (Green - SWIR2) - (0.25 * NIR + 2.75 * SWIR1)
+```
+
+AWEI can be implemented easily using `RS_MapAlgebra`:
+
+```sql
+-- Assume that the raster includes all 13 Sentinel-2 bands
+SELECT RS_MapAlgebra(rast, 'D', 'out = 4 * (rast[2] - rast[11]) - (0.25 * rast[7] + 2.75 * rast[12]);') as awei FROM raster_table
+```
+
+We can also implement the same AWEI calculation using array based map algebra functions. The code looks more verbose:
+
+```sql
+SELECT RS_Subtract(
+ RS_Add(RS_MultiplyFactor(band_nir, 0.25), RS_MultiplyFactor(band_swir1, 2.75)),
+ RS_MultiplyFactor(RS_Subtract(band_swir2, band_green), 4)) as awei
+FROM (
+SELECT RS_BandAsArray(rast, 3) AS band_green,
+ RS_BandAsArray(rast, 12) AS band_swir2,
+ RS_BandAsArray(rast, 13) AS band_swir1,
+ RS_BandAsArray(rast, 8) AS band_nir
+FROM raster_table) t
+```
+
+### Further Reading
+
+* [Jiffle language summary](https://github.com/geosolutions-it/jai-ext/wiki/Jiffle---language-summary)
+* [Raster operators](../Raster-operators/)
diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md
index b4e54df5c..05653732b 100644
--- a/docs/api/sql/Raster-operators.md
+++ b/docs/api/sql/Raster-operators.md
@@ -791,6 +791,41 @@ Output:
+--------------------+
```
+### RS_MapAlgebra
+
+Introduction: Apply a map algebra script on a raster.
+
+Format: `RS_MapAlgebra (raster: Raster, pixelType: String, script: String)`
+
+Format: `RS_MapAlgebra (raster: Raster, pixelType: String, script: String, noDataValue: Double)`
+
+Since: `v1.5.0`
+
+`RS_MapAlgebra` runs a script on a raster. The script is written in a map algebra language called [Jiffle](https://github.com/geosolutions-it/jai-ext/wiki/Jiffle). The script takes a raster
+as input and returns a raster of the same size as output. The script can be used to apply a map algebra expression on a raster. The input raster is named `rast` in the Jiffle script, and the output raster is named `out`.
+
+SQL example:
+
+Calculate the NDVI of a raster with 4 bands (R, G, B, NIR):
+
+```sql
+-- Assume that the input raster has 4 bands: R, G, B, NIR
+-- rast[0] refers to the R band, rast[3] refers to the NIR band.
+SELECT RS_MapAlgebra(rast, 'D', 'out = (rast[3] - rast[0]) / (rast[3] + rast[0]);') AS ndvi FROM raster_table
+```
+
+Output:
+```
++--------------------+
+| raster|
++--------------------+
+|GridCoverage2D["g...|
++--------------------+
+```
+
+For more details and examples about `RS_MapAlgebra`, please refer to the [Map Algebra documentation](../Raster-map-algebra/).
+To learn how to write map algebra script, please refer to [Jiffle language summary](https://github.com/geosolutions-it/jai-ext/wiki/Jiffle---language-summary).
+
## Map Algebra operators
Map algebra operators work on a single band of a raster. Each band is represented as an array of doubles. The operators return an array of doubles.
diff --git a/mkdocs.yml b/mkdocs.yml
index 44e04ef70..15c78e5d3 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -62,6 +62,7 @@ nav:
- Raster loader: api/sql/Raster-loader.md
- Raster writer: api/sql/Raster-writer.md
- Raster operators: api/sql/Raster-operators.md
+ - Raster map algebra: api/sql/Raster-map-algebra.md
- Parameter: api/sql/Parameter.md
- RDD (core):
- Scala/Java doc: api/java-api.md
diff --git a/pom.xml b/pom.xml
index 922bde107..3361565e8 100644
--- a/pom.xml
+++ b/pom.xml
@@ -74,6 +74,10 @@
<jts.version>1.19.0</jts.version>
<jts2geojson.version>0.16.1</jts2geojson.version>
+ <jt-jiffle.version>1.1.24</jt-jiffle.version>
+ <antlr-runtime.version>4.9.3</antlr-runtime.version>
+ <janino-version>3.1.9</janino-version>
+
<!-- Actual scala, spark and log4j version will be set by activated profiles.
Setting a default value helps IDE:s that can't make sense of profiles. -->
<scala.compat.version>2.12</scala.compat.version>
@@ -234,6 +238,21 @@
<version>${geotools.version}</version>
<scope>${geotools.scope}</scope>
</dependency>
+ <dependency>
+ <groupId>it.geosolutions.jaiext.jiffle</groupId>
+ <artifactId>jt-jiffle-language</artifactId>
+ <version>${jt-jiffle.version}</version>
+ <exclusions>
+ <exclusion>
+ <groupId>org.antlr</groupId>
+ <artifactId>*</artifactId>
+ </exclusion>
+ <exclusion>
+ <groupId>org.codehaus.janino</groupId>
+ <artifactId>*</artifactId>
+ </exclusion>
+ </exclusions>
+ </dependency>
<dependency>
<groupId>org.apache.hadoop</groupId>
<artifactId>hadoop-client</artifactId>
diff --git a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
index 2d2c22dda..337b1ceee 100644
--- a/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
+++ b/sql/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala
@@ -196,6 +196,7 @@ object Catalog {
function[RS_Append](),
function[RS_AddBandFromArray](),
function[RS_BandAsArray](),
+ function[RS_MapAlgebra](null),
function[RS_FromArcInfoAsciiGrid](),
function[RS_FromGeoTiff](),
function[RS_MakeEmptyRaster](),
diff --git a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala
index 438df4289..ae0e50e9b 100644
--- a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala
+++ b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/MapAlgebra.scala
@@ -20,13 +20,15 @@ package org.apache.spark.sql.sedona_sql.expressions.raster
import org.apache.sedona.common.raster.MapAlgebra
import org.apache.spark.sql.catalyst.InternalRow
-import org.apache.spark.sql.catalyst.expressions.{Expression, ImplicitCastInputTypes}
import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback
+import org.apache.spark.sql.catalyst.expressions.{Expression, ImplicitCastInputTypes}
import org.apache.spark.sql.catalyst.util.{ArrayData, GenericArrayData}
import org.apache.spark.sql.sedona_sql.expressions.InferrableFunctionConverter._
import org.apache.spark.sql.sedona_sql.expressions.{InferredExpression, UserDataGeneratator}
import org.apache.spark.sql.types._
+import javax.media.jai.RasterFactory
+
/// Calculate Normalized Difference between two bands
case class RS_NormalizedDifference(inputExpressions: Seq[Expression])
extends Expression with CodegenFallback with UserDataGeneratator {
@@ -816,3 +818,11 @@ case class RS_BandAsArray(inputExpressions: Seq[Expression]) extends InferredExp
copy(inputExpressions = newChildren)
}
}
+
+case class RS_MapAlgebra(inputExpressions: Seq[Expression])
+ extends InferredExpression(nullTolerantInferrableFunction4(MapAlgebra.mapAlgebra)) {
+
+ protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+ copy(inputExpressions = newChildren)
+ }
+}
diff --git a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala b/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala
index 08e83f6f8..65f8b99f7 100644
--- a/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala
+++ b/sql/common/src/test/scala/org/apache/sedona/sql/rasteralgebraTest.scala
@@ -18,6 +18,7 @@
*/
package org.apache.sedona.sql
+import org.apache.sedona.common.raster.MapAlgebra
import org.apache.sedona.common.utils.RasterUtils
import org.apache.spark.sql.functions.{collect_list, expr}
import org.geotools.coverage.grid.GridCoverage2D
@@ -25,6 +26,7 @@ import org.junit.Assert.{assertEquals, assertNull}
import org.locationtech.jts.geom.{Coordinate, Geometry}
import org.scalatest.{BeforeAndAfter, GivenWhenThen}
+import java.awt.image.DataBuffer
import scala.collection.mutable
@@ -731,5 +733,71 @@ class rasteralgebraTest extends TestBaseScala with BeforeAndAfter with GivenWhen
val result = df.selectExpr("RS_BandPixelType(raster)").first().getString(0)
assertEquals("UNSIGNED_8BITS", result)
}
+
+ it("Passed RS_MapAlgebra") {
+ val df = sparkSession.read.format("binaryFile")
+ .option("recursiveFileLookup", "true")
+ .option("pathGlobFilter", "*.tif*")
+ .load(resourceFolder + "raster")
+ .selectExpr("path", "RS_FromGeoTiff(content) as rast")
+ Seq(null, "b", "s", "i", "f", "d").foreach { pixelType =>
+ val pixelTypeExpr = if (pixelType == null) null else s"'$pixelType'"
+ val dfResult = df.withColumn("rast_2", expr(s"RS_MapAlgebra(rast, $pixelTypeExpr, 'out[0] = rast[0] * 0.5;')"))
+ .select("path", "rast", "rast_2")
+ dfResult.collect().foreach { row =>
+ val rast = row.getAs[GridCoverage2D]("rast")
+ val rast2 = row.getAs[GridCoverage2D]("rast_2")
+ assert(rast.getGridGeometry.getGridToCRS2D == rast2.getGridGeometry.getGridToCRS2D)
+ val dataType = pixelType match {
+ case null => rast.getRenderedImage.getSampleModel.getDataType
+ case _ => RasterUtils.getDataTypeCode(pixelType)
+ }
+ assert(rast2.getRenderedImage.getSampleModel.getDataType == dataType)
+ val noDataValue = RasterUtils.getNoDataValue(rast2.getSampleDimension(0))
+ assert(noDataValue.isNaN)
+ val band = MapAlgebra.bandAsArray(rast, 1)
+ val band2 = MapAlgebra.bandAsArray(rast2, 1)
+ assert(band.size == band2.size)
+ for (i <- band.indices) {
+ pixelType match {
+ case "b" => assert(((band(i) * 0.5).toInt & 0xFF) == band2(i).toInt)
+ case "s" => assert((band(i) * 0.5).toShort == band2(i))
+ case "i" | null => assert((band(i) * 0.5).toInt == band2(i))
+ case "f" | "d" =>
+ if (band(i) != 0) {
+ assert((band(i) * 0.5) == band2(i))
+ } else {
+ // If the source image has NoDataContainer.GC_NODATA property, Jiffle may convert nodata values in
+ // source raster to NaN.
+ assert(band2(i) == 0 || band2(i).isNaN)
+ }
+ }
+ }
+ }
+ }
+ }
+
+ it("Passed RS_MapAlgebra with nodata value") {
+ val df = sparkSession.read.format("binaryFile")
+ .option("recursiveFileLookup", "true")
+ .option("pathGlobFilter", "*.tif*")
+ .load(resourceFolder + "raster")
+ .selectExpr("RS_FromGeoTiff(content) as rast")
+ val dfResult = df.withColumn("rast_2", expr("RS_MapAlgebra(rast, 'us', 'out[0] = rast[0] * 0.2;', 10)"))
+ .select("rast", "rast_2")
+ dfResult.collect().foreach { row =>
+ val rast = row.getAs[GridCoverage2D]("rast")
+ val rast2 = row.getAs[GridCoverage2D]("rast_2")
+ assert(rast2.getRenderedImage.getSampleModel.getDataType == DataBuffer.TYPE_USHORT)
+ val noDataValue = RasterUtils.getNoDataValue(rast2.getSampleDimension(0))
+ assert(noDataValue == 10)
+ val band = MapAlgebra.bandAsArray(rast, 1)
+ val band2 = MapAlgebra.bandAsArray(rast2, 1)
+ assert(band.size == band2.size)
+ for (i <- band.indices) {
+ assert((band(i) * 0.2).toShort == band2(i))
+ }
+ }
+ }
}
}