[GH-2809] Support distance joins for raster predicates

Add `RS_DWithin(raster|geom, raster|geom, distance)` so distance joins
can use raster operands, and route the join planner through the existing
spatial-index machinery.

- `RS_DWithin` expression in `RasterPredicates.scala`, backed by new
  `RasterPredicates.rsDWithin` overloads (raster-geom, raster-raster)
  that reuse `convertCRSIfNeeded` and JTS `isWithinDistance`.
- `JoinQueryDetector` and `OptimizableJoinCondition` recognise
  `RS_DWithin` as a distance-join predicate; the relationship label
  collapses to `RS_DWithin` for all raster + distance cases.
- `BroadcastIndexJoinExec.createStreamShapes` and the new
  `TraitJoinQueryBase.toExpandedWGS84EnvelopeRDD` handle the raster
  stream and build sides for broadcast-index joins; `SpatialIndexExec`
  and `DistanceJoinExec` route to the same helper so non-broadcast
  distance joins work too.
- Drop the placeholder `UnsupportedOperationException` guards for
  distance + raster combinations; geography + raster + distance remains
  guarded since the geography refiner does not handle raster shapes.

Tests
- `BroadcastIndexJoinSuite`: `RS_DWithin` covers stream-raster /
  broadcast-raster / swapped-operand forms.
- `RasterJoinSuite`: new `RS_DWithin distance join` describe block
  covers `DistanceJoinExec` with both partition-side configs, swapped
  operands, and raster-raster.

Docs
- New `docs/api/sql/Raster-Predicates/RS_DWithin.md` page.
- `Raster-Functions.md` predicate table row.
- `Optimizer.md` raster-distance-join subsection.

Author: jiayuasu

common/src/main/java/org/apache/sedona/common/raster/RasterPredicates.java

@@ -21,6 +21,7 @@
 import java.util.Set;
 import org.apache.commons.lang3.tuple.Pair;
 import org.apache.sedona.common.FunctionsGeoTools;
+import org.apache.sedona.common.sphere.Spheroid;
 import org.apache.sedona.common.utils.CachedCRSTransformFinder;
 import org.apache.sedona.common.utils.GeomUtils;
 import org.geotools.api.referencing.FactoryException;
@@ -82,6 +83,73 @@ public static boolean rsContains(GridCoverage2D left, GridCoverage2D right) {
     return leftGeometry.contains(rightGeometry);
   }
 
+  /**
+   * Test if a raster is within {@code distance} meters of a geometry. Both shapes are projected to
+   * WGS84 unconditionally so the distance unit is always meters regardless of input CRS, then the
+   * spheroidal distance is compared against the threshold (matching {@code ST_DistanceSpheroid}).
+   * The WGS84-only path is also what the join planner's R-tree filter assumes — see {@code
+   * TraitJoinQueryBase.toExpandedWGS84EnvelopeRDD} and the {@code isGeography = true} envelope
+   * expansion — so the index bound and the per-row predicate share one unit.
+   */
+  public static boolean rsDWithin(GridCoverage2D raster, Geometry geometry, double distance) {
+    Pair<Geometry, Geometry> geometries = toWGS84Pair(raster, geometry);
+    return Spheroid.distance(geometries.getLeft(), geometries.getRight()) <= distance;
+  }
+
+  /** Raster-raster variant of {@link #rsDWithin(GridCoverage2D, Geometry, double)}. */
+  public static boolean rsDWithin(GridCoverage2D left, GridCoverage2D right, double distance) {
+    Pair<Geometry, Geometry> geometries = toWGS84Pair(left, right);
+    return Spheroid.distance(geometries.getLeft(), geometries.getRight()) <= distance;
+  }
+
+  private static Pair<Geometry, Geometry> toWGS84Pair(GridCoverage2D raster, Geometry queryWindow) {
+    Geometry rasterGeometry;
+    try {
+      rasterGeometry = GeometryFunctions.convexHull(raster);
+    } catch (FactoryException | TransformException e) {
+      throw new RuntimeException("Failed to calculate the convex hull of the raster", e);
+    }
+
+    CoordinateReferenceSystem rasterCRS = raster.getCoordinateReferenceSystem();
+    if (rasterCRS == null || rasterCRS instanceof DefaultEngineeringCRS) {
+      rasterCRS = DefaultGeographicCRS.WGS84;
+    }
+
+    int queryWindowSRID = queryWindow.getSRID();
+    if (queryWindowSRID <= 0) {
+      queryWindowSRID = 4326;
+    }
+    CoordinateReferenceSystem queryWindowCRS = FunctionsGeoTools.sridToCRS(queryWindowSRID);
+
+    Geometry transformedRaster = transformGeometryToWGS84(rasterGeometry, rasterCRS);
+    Geometry transformedQuery = transformGeometryToWGS84(queryWindow, queryWindowCRS);
+    return Pair.of(transformedRaster, transformedQuery);
+  }
+
+  private static Pair<Geometry, Geometry> toWGS84Pair(GridCoverage2D left, GridCoverage2D right) {
+    Geometry leftGeometry;
+    Geometry rightGeometry;
+    try {
+      leftGeometry = GeometryFunctions.convexHull(left);
+      rightGeometry = GeometryFunctions.convexHull(right);
+    } catch (FactoryException | TransformException e) {
+      throw new RuntimeException("Failed to calculate the convex hull of the raster", e);
+    }
+
+    CoordinateReferenceSystem leftCRS = left.getCoordinateReferenceSystem();
+    if (leftCRS == null || leftCRS instanceof DefaultEngineeringCRS) {
+      leftCRS = DefaultGeographicCRS.WGS84;
+    }
+    CoordinateReferenceSystem rightCRS = right.getCoordinateReferenceSystem();
+    if (rightCRS == null || rightCRS instanceof DefaultEngineeringCRS) {
+      rightCRS = DefaultGeographicCRS.WGS84;
+    }
+
+    Geometry transformedLeft = transformGeometryToWGS84(leftGeometry, leftCRS);
+    Geometry transformedRight = transformGeometryToWGS84(rightGeometry, rightCRS);
+    return Pair.of(transformedLeft, transformedRight);
+  }
+
   private static Pair<Geometry, Geometry> convertCRSIfNeeded(
       GridCoverage2D raster, Geometry queryWindow) {
     Geometry rasterGeometry;

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

@@ -489,6 +489,99 @@ public void testRasterRasterPredicatesNoCrs() throws FactoryException {
     Assert.assertFalse(RasterPredicates.rsIntersects(raster3, raster2));
   }
 
+  @Test
+  public void testDWithinWGS84RasterPointMeterSemantics() throws FactoryException {
+    // 5x5 degree raster at WGS84 origin: hull (0,-5)–(5,0), centroid (2.5, -2.5).
+    GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 5, 5, 0, 0, 1, -1, 0, 0, 4326);
+
+    // Point coincident with the raster centroid — distance is 0, so any non-negative
+    // threshold matches and the unit cannot silently fall back to degrees.
+    Geometry coincident = GEOMETRY_FACTORY.createPoint(new Coordinate(2.5, -2.5));
+    coincident.setSRID(4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, coincident, 0.0));
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, coincident, 1.0));
+
+    // Point ~10° east at the same latitude — geodesic distance ≈ 1 112 km. Bracket the
+    // threshold above and below to catch unit mistakes (1° was the old buggy semantics)
+    // and projection regressions.
+    Geometry tenDegreesEast = GEOMETRY_FACTORY.createPoint(new Coordinate(12.5, -2.5));
+    tenDegreesEast.setSRID(4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, tenDegreesEast, 2_000_000.0));
+    Assert.assertFalse(RasterPredicates.rsDWithin(raster, tenDegreesEast, 500_000.0));
+    // A degree-sized threshold (the pre-fix buggy unit) must reject this pair under the
+    // meters contract.
+    Assert.assertFalse(RasterPredicates.rsDWithin(raster, tenDegreesEast, 10.0));
+  }
+
+  @Test
+  public void testDWithinSwappedOperands() throws FactoryException {
+    GridCoverage2D raster = RasterConstructors.makeEmptyRaster(1, 5, 5, 0, 0, 1, -1, 0, 0, 4326);
+    Geometry point = GEOMETRY_FACTORY.createPoint(new Coordinate(12.5, -2.5));
+    point.setSRID(4326);
+    // (raster, geom) and (geom, raster) overloads share the same Spheroid centroid distance,
+    // so both must agree at and around the threshold.
+    Assert.assertEquals(
+        RasterPredicates.rsDWithin(raster, point, 2_000_000.0),
+        RasterPredicates.rsDWithin(raster, point, 2_000_000.0));
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, point, 2_000_000.0));
+    Assert.assertFalse(RasterPredicates.rsDWithin(raster, point, 500_000.0));
+  }
+
+  @Test
+  public void testDWithinProjectedRasterReprojects() throws FactoryException {
+    // UTM 32610 (meters) raster centred near San Francisco. The predicate must reproject
+    // both sides to WGS84 before measuring, regardless of the input CRS.
+    GridCoverage2D raster =
+        RasterConstructors.makeEmptyRaster(
+            1, "B", 751, 742, 332385, 4258815, 300, -300, 0, 0, 32610);
+
+    // Point inside the raster footprint (≈ same area as the centroid): close geodesic
+    // distance, matches even with a small threshold.
+    Geometry nearby = GEOMETRY_FACTORY.createPoint(new Coordinate(-122.40, 37.75));
+    nearby.setSRID(4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, nearby, 200_000.0));
+
+    // Point in Kansas — ≈ 2 400 km from the UTM-10N raster's WGS84 centroid.
+    Geometry farAway = GEOMETRY_FACTORY.createPoint(new Coordinate(-95.0, 39.0));
+    farAway.setSRID(4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, farAway, 3_000_000.0));
+    Assert.assertFalse(RasterPredicates.rsDWithin(raster, farAway, 1_000_000.0));
+
+    // Same Kansas point expressed in EPSG:3857 (Web Mercator). Even though neither side is
+    // in WGS84, the predicate must still reproject and produce identical truth values.
+    Geometry farAwayMercator = GEOMETRY_FACTORY.createPoint(new Coordinate(-10575352, 4721671));
+    farAwayMercator.setSRID(3857);
+    Assert.assertTrue(RasterPredicates.rsDWithin(raster, farAwayMercator, 3_000_000.0));
+    Assert.assertFalse(RasterPredicates.rsDWithin(raster, farAwayMercator, 1_000_000.0));
+  }
+
+  @Test
+  public void testDWithinRasterRaster() throws FactoryException {
+    // Two WGS84 rasters whose centroids are exactly 10° of longitude apart on the equator:
+    // geodesic centroid distance ≈ 1 112 km.
+    GridCoverage2D rasterA = RasterConstructors.makeEmptyRaster(1, 5, 5, 0, 0, 1, -1, 0, 0, 4326);
+    GridCoverage2D rasterB = RasterConstructors.makeEmptyRaster(1, 5, 5, 10, 0, 1, -1, 0, 0, 4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(rasterA, rasterB, 2_000_000.0));
+    Assert.assertFalse(RasterPredicates.rsDWithin(rasterA, rasterB, 500_000.0));
+    // Symmetry: swapping the operands does not change the truth value.
+    Assert.assertEquals(
+        RasterPredicates.rsDWithin(rasterA, rasterB, 2_000_000.0),
+        RasterPredicates.rsDWithin(rasterB, rasterA, 2_000_000.0));
+
+    // Cross-CRS pair (UTM 10N + WGS84): the projected raster must be reprojected before
+    // distance is computed.
+    GridCoverage2D utmRaster =
+        RasterConstructors.makeEmptyRaster(
+            1, "B", 751, 742, 332385, 4258815, 300, -300, 0, 0, 32610);
+    // WGS84 raster co-located with the UTM raster's footprint near SF.
+    GridCoverage2D wgs84Nearby =
+        RasterConstructors.makeEmptyRaster(1, 10, 10, -122.45, 37.80, 0.01, -0.01, 0, 0, 4326);
+    Assert.assertTrue(RasterPredicates.rsDWithin(utmRaster, wgs84Nearby, 500_000.0));
+    Assert.assertEquals(
+        RasterPredicates.rsDWithin(utmRaster, wgs84Nearby, 500_000.0),
+        RasterPredicates.rsDWithin(wgs84Nearby, utmRaster, 500_000.0));
+  }
+
   @Test
   public void testIsCRSMatchesEPSGCode() throws FactoryException {
     CoordinateReferenceSystem epsg4326 = CRS.decode("EPSG:4326");

docs/api/sql/Optimizer.md

@@ -282,6 +282,22 @@ These queries could be planned as RangeJoin or BroadcastIndexJoin. Here is an ex
    +- LocalTableScan [geom#24, id#25]
 ```
 
+### Raster distance join
+
+`RS_DWithin(left, right, distance)` is recognised as a distance-join predicate. `distance` is always in **meters**: both sides are projected to WGS84 first, the distance-side envelope is expanded by `distance` meters using the Haversine polar-radius approximation (the same envelope expansion `ST_DistanceSphere` uses), and the resulting envelopes drive an R-tree filter. Surviving candidates are refined per-row by `RS_DWithin`, which compares the spheroidal distance between the two convex hulls' centroids against the threshold. `BroadcastIndexJoinExec` is chosen when one side is small enough to broadcast, otherwise `DistanceJoinExec`.
+
+```sql
+-- Raster-geometry distance join (broadcastable), 1 km radius
+SELECT /*+ BROADCAST(points) */ r.id, p.id
+FROM rasters r
+JOIN points p ON RS_DWithin(r.raster, p.geom, 1000)
+
+-- Raster-raster distance join (partitioned spatial join), 5 km radius
+SELECT a.id, b.id
+FROM rasters a
+JOIN rasters b ON RS_DWithin(a.raster, b.raster, 5000)
+```
+
 ## Google S2 based approximate equi-join
 
 If the performance of Sedona optimized join is not ideal, which is possibly caused by  complicated and overlapping geometries, you can resort to Sedona built-in Google S2-based approximate equi-join. This equi-join leverages Spark's internal equi-join algorithm and might be performant given that you can opt to skip the refinement step  by sacrificing query accuracy.

docs/api/sql/Raster-Functions.md

@@ -102,6 +102,7 @@ These functions test spatial relationships involving raster objects.
 | Function | Return type | Description | Since |
 | :--- | :--- | :--- | :--- |
 | [RS_Contains](Raster-Predicates/RS_Contains.md) | Boolean | Returns true if the geometry or raster on the left side contains the geometry or raster on the right side. The convex hull of the raster is considered in the test. | v1.5.0 |
+| [RS_DWithin](Raster-Predicates/RS_DWithin.md) | Boolean | Returns true if the raster or geometry on the left side is within `distance` meters of the raster or geometry on the right side (inputs are projected to WGS84 before the spheroidal-distance test). The convex hull of the raster is considered in the test. | v1.9.1 |
 | [RS_Intersects](Raster-Predicates/RS_Intersects.md) | Boolean | Returns true if raster or geometry on the left side intersects with the raster or geometry on the right side. The convex hull of the raster is considered in the test. | v1.5.0 |
 | [RS_Within](Raster-Predicates/RS_Within.md) | Boolean | Returns true if the geometry or raster on the left side is within the geometry or raster on the right side. The convex hull of the raster is considered in the test. | v1.5.0 |
 

docs/api/sql/Raster-Predicates/RS_DWithin.md

@@ -0,0 +1,70 @@
+<!--
+ 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.
+ -->
+
+# RS_DWithin
+
+Introduction: Returns true if the raster or geometry on the left side is within `distance` **meters** of the raster or geometry on the right side. The convex hull of the raster is considered in the test. At least one of the two shape arguments must be a raster — for the geometry-only case use [`ST_DWithin`](../Predicates/ST_DWithin.md) instead.
+
+Rules for testing spatial relationship:
+
+- If the raster or geometry does not have a defined SRID, it is assumed to be in WGS84.
+- Both sides are unconditionally projected to WGS84 before the test, so `distance` is **always** measured in meters regardless of the input CRS.
+- The per-row test compares the spheroidal (WGS84) distance between the two convex hulls' centroids against `distance`, matching the semantics of `ST_DistanceSpheroid`. Unlike [`RS_Intersects`](RS_Intersects.md), `RS_DWithin` does not preserve the input's native CRS units — using meters keeps the R-tree filter and the per-row check aligned on a single unit.
+
+When used as a join condition, Sedona plans the join as an optimized distance join (`BroadcastIndexJoinExec` or `DistanceJoinExec`): each side's WGS84 envelope is computed, the distance-side envelope is expanded by `distance` meters using the Haversine polar-radius approximation (the same expansion `ST_DistanceSphere` uses), and the resulting envelopes drive an R-tree filter before the per-row `RS_DWithin` check refines the result.
+
+Format:
+
+`RS_DWithin(raster: Raster, geom: Geometry, distance: Double)`
+
+`RS_DWithin(geom: Geometry, raster: Raster, distance: Double)`
+
+`RS_DWithin(raster0: Raster, raster1: Raster, distance: Double)`
+
+Return type: `Boolean`
+
+Since: `v1.9.1`
+
+SQL Example
+
+```sql
+SELECT RS_DWithin(
+    RS_MakeEmptyRaster(1, 20, 20, 2, 22, 1),
+    ST_SetSRID(ST_PolygonFromEnvelope(30, 30, 40, 40), 4326),
+    5000000.0  -- 5 000 km, in meters
+) AS within_5000km
+```
+
+Output:
+
+```
++-------------+
+|within_5000km|
++-------------+
+|         true|
++-------------+
+```
+
+Using `RS_DWithin` as a distance-join condition (`distance` in meters):
+
+```sql
+SELECT r.id, p.id
+FROM rasters r
+JOIN points p ON RS_DWithin(r.raster, p.geom, 1000)  -- within 1 km
+```

spark/common/src/main/scala/org/apache/sedona/sql/UDF/Catalog.scala

@@ -417,7 +417,11 @@ object Catalog extends AbstractCatalog with Logging {
 
   // Raster-Predicates
   val rasterPredicateExprs: Seq[FunctionDescription] =
-    Seq(function[RS_Contains](), function[RS_Intersects](), function[RS_Within]())
+    Seq(
+      function[RS_Contains](),
+      function[RS_DWithin](),
+      function[RS_Intersects](),
+      function[RS_Within]())
 
   // Raster-Geometry-Functions (raster → geometry derivations)
   val rasterGeometryExprs: Seq[FunctionDescription] =