[SEDONA-321] Implement RS_Intersects(rast, geom)

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);
+    }
+}

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);
+    }
+}

common/src/test/java/org/apache/sedona/common/raster/RasterTestBase.java

@@ -17,13 +17,18 @@
 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 @@ GridCoverage2D createEmptyRaster(int numBands)
         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);

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")
 
-```
+```

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]()
   )

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)
+  }
+}

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))
+    }
   }
 }