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 01414111 [SEDONA-321] Implement RS_Intersects(rast, geom) (#893)
01414111 is described below
commit 01414111b0e5ab8715428a72a3f03377c77dc456
Author: Kristin Cowalcijk <[email protected]>
AuthorDate: Fri Jul 7 12:26:08 2023 +0800
[SEDONA-321] Implement RS_Intersects(rast, geom) (#893)
---
.../sedona/common/raster/RasterPredicates.java | 65 +++++++++++++++
.../sedona/common/raster/RasterPredicatesTest.java | 96 ++++++++++++++++++++++
.../sedona/common/raster/RasterTestBase.java | 37 +++++++++
docs/api/sql/Raster-operators.md | 21 ++++-
.../scala/org/apache/sedona/sql/UDF/Catalog.scala | 1 +
.../expressions/raster/RasterPredicates.scala | 50 +++++++++++
.../org/apache/sedona/sql/rasteralgebraTest.scala | 17 ++++
7 files changed, 286 insertions(+), 1 deletion(-)
diff --git
a/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java
b/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java
new file mode 100644
index 00000000..ec492df9
--- /dev/null
+++ b/common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java
@@ -0,0 +1,65 @@
+/*
+ * 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 org.geotools.coverage.grid.GridCoverage2D;
+import org.geotools.geometry.Envelope2D;
+import org.geotools.geometry.jts.JTS;
+import org.geotools.referencing.CRS;
+import org.geotools.referencing.crs.DefaultEngineeringCRS;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+import org.opengis.referencing.FactoryException;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.TransformException;
+
+public class RasterPredicates {
+ private static final GeometryFactory GEOMETRY_FACTORY = new
GeometryFactory();
+
+ /**
+ * Test if a raster intersects a query window. If both the raster and the
query window have a
+ * CRS, the query window will be transformed to the CRS of the raster
before testing for intersection.
+ * Please note that the CRS transformation will be lenient, which means
that the transformation may
+ * not be accurate.
+ * @param raster the raster
+ * @param queryWindow the query window
+ * @return true if the raster intersects the query window
+ */
+ public static boolean rsIntersects(GridCoverage2D raster, Geometry
queryWindow) {
+ Envelope2D rasterEnvelope2D = raster.getEnvelope2D();
+ CoordinateReferenceSystem rasterCRS =
rasterEnvelope2D.getCoordinateReferenceSystem();
+ int queryWindowSRID = queryWindow.getSRID();
+ if (rasterCRS != null && !(rasterCRS instanceof DefaultEngineeringCRS)
&& queryWindowSRID > 0) {
+ try {
+ CoordinateReferenceSystem queryWindowCRS = CRS.decode("EPSG:"
+ queryWindowSRID);
+ if (!CRS.equalsIgnoreMetadata(rasterCRS, queryWindowCRS)) {
+ MathTransform transform =
CRS.findMathTransform(queryWindowCRS, rasterCRS, true);
+ queryWindow = JTS.transform(queryWindow, transform);
+ }
+ } catch (FactoryException | TransformException e) {
+ throw new RuntimeException("Cannot transform CRS of query
window", e);
+ }
+ }
+ Envelope rasterEnvelope = JTS.toEnvelope(rasterEnvelope2D);
+ Geometry rasterGeometry = GEOMETRY_FACTORY.toGeometry(rasterEnvelope);
+ return rasterGeometry.intersects(queryWindow);
+ }
+}
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java
b/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java
new file mode 100644
index 00000000..2c18059e
--- /dev/null
+++
b/common/src/test/java/org/apache/sedona/common/raster/RasterPredicatesTest.java
@@ -0,0 +1,96 @@
+/*
+ * 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 org.geotools.coverage.grid.GridCoverage2D;
+import org.junit.Assert;
+import org.junit.Test;
+import org.locationtech.jts.geom.Envelope;
+import org.locationtech.jts.geom.Geometry;
+import org.locationtech.jts.geom.GeometryFactory;
+
+import java.awt.image.DataBuffer;
+
+public class RasterPredicatesTest extends RasterTestBase {
+ private static final GeometryFactory GEOMETRY_FACTORY = new
GeometryFactory();
+ @Test
+ public void testIntersectsNoCrs() {
+ Geometry queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(0, 10,
0, 10));
+ GridCoverage2D raster = createRandomRaster(DataBuffer.TYPE_BYTE, 100,
100, 0, 100, 1, 1, null);
+ boolean result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(1000, 1010,
1000, 1010));
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertFalse(result);
+ }
+
+ @Test
+ public void testIntersectsQueryWindowNoCrs() {
+ Geometry queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(0, 10,
0, 10));
+ GridCoverage2D raster = createRandomRaster(DataBuffer.TYPE_BYTE, 100,
100, 0, 100, 1, 1, "EPSG:3857");
+ boolean result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(1000, 1010,
1000, 1010));
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertFalse(result);
+ }
+
+ @Test
+ public void testIntersectsRasterNoCrs() {
+ Geometry queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(0, 10,
0, 10));
+ queryWindow.setSRID(3857);
+ GridCoverage2D raster = createRandomRaster(DataBuffer.TYPE_BYTE, 100,
100, 0, 100, 1, 1, null);
+ boolean result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(1000, 1010,
1000, 1010));
+ queryWindow.setSRID(3857);
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertFalse(result);
+ }
+
+ @Test
+ public void testIntersectsSameCrs() {
+ Geometry queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(0, 10,
0, 10));
+ queryWindow.setSRID(3857);
+ GridCoverage2D raster = createRandomRaster(DataBuffer.TYPE_BYTE, 100,
100, 0, 100, 1, 1, "EPSG:3857");
+ boolean result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(10, 20, 10,
20));
+ queryWindow.setSRID(3857);
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(1000, 1010,
1000, 1010));
+ queryWindow.setSRID(3857);
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertFalse(result);
+ }
+
+ @Test
+ public void testIntersectsWithTransformations() {
+ Geometry queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(0, 10,
0, 10));
+ queryWindow.setSRID(4326);
+ GridCoverage2D raster = createRandomRaster(DataBuffer.TYPE_BYTE, 100,
100, 0, 100, 1, 1, "EPSG:3857");
+ boolean result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertTrue(result);
+ queryWindow = GEOMETRY_FACTORY.toGeometry(new Envelope(10, 20, 10,
20));
+ queryWindow.setSRID(4326);
+ result = RasterPredicates.rsIntersects(raster, queryWindow);
+ Assert.assertFalse(result);
+ }
+}
diff --git
a/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java
b/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java
index bccb4c7c..d75a34b9 100644
--- a/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java
+++ b/common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java
@@ -17,13 +17,18 @@ import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.GridCoverageFactory;
import org.geotools.gce.geotiff.GeoTiffWriter;
import org.geotools.geometry.Envelope2D;
+import org.geotools.geometry.jts.ReferencedEnvelope;
+import org.geotools.referencing.CRS;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.junit.Before;
import org.opengis.parameter.GeneralParameterValue;
import org.opengis.referencing.FactoryException;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import javax.media.jai.RasterFactory;
import java.awt.Color;
import java.awt.image.BufferedImage;
+import java.awt.image.WritableRaster;
import java.io.ByteArrayOutputStream;
import java.io.IOException;
import java.nio.charset.StandardCharsets;
@@ -53,6 +58,38 @@ public class RasterTestBase {
double cellSize = 2;
return RasterConstructors.makeEmptyRaster(numBands, widthInPixel,
heightInPixel, upperLeftX, upperLeftY, cellSize);
}
+
+ GridCoverage2D createRandomRaster(int dataBufferType, int widthInPixel,
int heightInPixel,
+ double upperLeftX, double upperLeftY,
double pixelSize,
+ int numBand, String crsCode) {
+ WritableRaster raster =
+ RasterFactory.createBandedRaster(
+ dataBufferType, widthInPixel, heightInPixel, numBand,
null);
+ for (int b = 0; b < numBand; b++) {
+ for (int y = 0; y < heightInPixel; y++) {
+ for (int x = 0; x < widthInPixel; x++) {
+ double value = Math.random() * 255;
+ raster.setSample(x, y, b, value);
+ }
+ }
+ }
+ CoordinateReferenceSystem crs = null;
+ if (crsCode != null && !crsCode.isEmpty()) {
+ try {
+ crs = CRS.decode(crsCode);
+ } catch (FactoryException e) {
+ throw new RuntimeException(e);
+ }
+ }
+ double x1 = upperLeftX;
+ double x2 = upperLeftX + widthInPixel * pixelSize;
+ double y1 = upperLeftY - heightInPixel * pixelSize;
+ double y2 = upperLeftY;
+ ReferencedEnvelope referencedEnvelope = new ReferencedEnvelope(x1, x2,
y1, y2, crs);
+ GridCoverageFactory factory = new GridCoverageFactory();
+ return factory.create("random-raster", raster, referencedEnvelope);
+ }
+
GridCoverage2D createMultibandRaster() throws IOException {
GridCoverageFactory factory = new GridCoverageFactory();
BufferedImage image = new BufferedImage(10, 10,
BufferedImage.TYPE_INT_ARGB);
diff --git a/docs/api/sql/Raster-operators.md b/docs/api/sql/Raster-operators.md
index fcaddaa3..6cbff57d 100644
--- a/docs/api/sql/Raster-operators.md
+++ b/docs/api/sql/Raster-operators.md
@@ -17,6 +17,25 @@ Output:
POLYGON((0 0,20 0,20 60,0 60,0 0))
```
+### RS_Intersects
+
+Introduction: Returns true if the envelope of the raster intersects the given
geometry. If the geometry does not have a
+defined SRID, it is considered to be in the same CRS with the raster. If the
geometry has a defined SRID, the geometry
+will be transformed to the CRS of the raster before the intersection test.
+
+Format: `RS_Intersects (raster: Raster, geom: Geometry)`
+
+Since: `v1.5.0`
+
+Spark SQL example:
+```sql
+SELECT RS_Intersects(raster, ST_SetSRID(ST_PolygonFromEnvelope(0, 0, 10, 10),
4326)) FROM raster_table
+```
+Output:
+```
+true
+```
+
### RS_MetaData
Introduction: Returns the metadata of the raster as an array of double. The
array contains the following values:
@@ -563,4 +582,4 @@ Spark SQL example:
val subtractDF = spark.sql("select RS_Subtract(band1, band2) as
differenceOfOfBands from dataframe")
-```
\ No newline at end of file
+```
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 e472fde6..a66d4cae 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
@@ -200,6 +200,7 @@ object Catalog {
function[RS_SRID](),
function[RS_Value](1),
function[RS_Values](1),
+ function[RS_Intersects](),
function[RS_AsGeoTiff](),
function[RS_AsArcGrid]()
)
diff --git
a/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterPredicates.scala
b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterPredicates.scala
new file mode 100644
index 00000000..d35a0256
--- /dev/null
+++
b/sql/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterPredicates.scala
@@ -0,0 +1,50 @@
+/*
+ * 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.spark.sql.sedona_sql.expressions.raster
+
+import org.apache.sedona.common.raster.RasterPredicates
+import org.apache.spark.sql.catalyst.InternalRow
+import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes,
Expression}
+import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback
+import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT}
+import
org.apache.spark.sql.sedona_sql.expressions.implicits.InputExpressionEnhancer
+import
org.apache.spark.sql.sedona_sql.expressions.raster.implicits.RasterInputExpressionEnhancer
+import org.apache.spark.sql.types.{AbstractDataType, BooleanType, DataType}
+
+case class RS_Intersects(inputExpressions: Seq[Expression]) extends Expression
with CodegenFallback with ExpectsInputTypes {
+
+ override def eval(input: InternalRow): Any = {
+ val raster = inputExpressions.head.toRaster(input)
+ val geom = inputExpressions(1).toGeometry(input)
+ if (raster == null || geom == null) {
+ null
+ } else {
+ RasterPredicates.rsIntersects(raster, geom)
+ }
+ }
+
+ override def nullable: Boolean = true
+ override def dataType: DataType = BooleanType
+ override def inputTypes: Seq[AbstractDataType] = Seq(RasterUDT, GeometryUDT)
+ override def children: Seq[Expression] = inputExpressions
+
+ 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 0c35e680..8ef24399 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
@@ -435,5 +435,22 @@ class rasteralgebraTest extends TestBaseScala with
BeforeAndAfter with GivenWhen
df.selectExpr("RS_BandAsArray(RS_AddBandFromArray(raster, band, 100),
1)").collect()
}
}
+
+ it("Passed RS_Intersects") {
+ val df = sparkSession.read.format("binaryFile").load(resourceFolder +
"raster/test1.tiff")
+ .selectExpr("path", "RS_FromGeoTiff(content) as raster")
+
+ // query window without SRID
+ assert(df.selectExpr("RS_Intersects(raster,
ST_Point(-13076178,4003651))").first().getBoolean(0))
+ assert(!df.selectExpr("RS_Intersects(raster,
ST_Point(-13055247,3979620))").first().getBoolean(0))
+
+ // query window and raster are in the same CRS
+ assert(df.selectExpr("RS_Intersects(raster,
ST_SetSRID(ST_Point(-13067806,4009116), 3857))").first().getBoolean(0))
+ assert(!df.selectExpr("RS_Intersects(raster,
ST_SetSRID(ST_Point(-13057457,4027023), 3857))").first().getBoolean(0))
+
+ // query window and raster not in the same CRS
+ assert(df.selectExpr("RS_Intersects(raster,
ST_SetSRID(ST_Point(33.81798,-117.47993), 4326))").first().getBoolean(0))
+ assert(!df.selectExpr("RS_Intersects(raster,
ST_SetSRID(ST_Point(33.97896,-117.27868), 4326))").first().getBoolean(0))
+ }
}
}
