[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

@@ -82,6 +82,22 @@ public static boolean rsContains(GridCoverage2D left, GridCoverage2D right) {
     return leftGeometry.contains(rightGeometry);
   }
 
+  /**
+   * Test if a raster is within {@code distance} of a geometry. Both shapes are converted to a
+   * common CRS (matching the rest of the RS_* predicate family); the distance is then measured in
+   * that CRS's units (degrees for WGS84, the projected unit otherwise).
+   */
+  public static boolean rsDWithin(GridCoverage2D raster, Geometry geometry, double distance) {
+    Pair<Geometry, Geometry> geometries = convertCRSIfNeeded(raster, geometry);
+    return geometries.getLeft().isWithinDistance(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 = convertCRSIfNeeded(left, right);
+    return geometries.getLeft().isWithinDistance(geometries.getRight(), distance);
+  }
+
   private static Pair<Geometry, Geometry> convertCRSIfNeeded(
       GridCoverage2D raster, Geometry queryWindow) {
     Geometry rasterGeometry;

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. Sedona converts each side to a WGS84 envelope, expands the envelope on the distance side by `distance`, and runs an R-tree filter before applying the per-row `RS_DWithin` check. The same plan is chosen as for `ST_DWithin` — `BroadcastIndexJoinExec` when one side is small enough to broadcast, otherwise `DistanceJoinExec`.
+
+```sql
+-- Raster-geometry distance join (broadcastable)
+SELECT /*+ BROADCAST(points) */ r.id, p.id
+FROM rasters r
+JOIN points p ON RS_DWithin(r.raster, p.geom, 0.01)
+
+-- Raster-raster distance join (partitioned spatial join)
+SELECT a.id, b.id
+FROM rasters a
+JOIN rasters b ON RS_DWithin(a.raster, b.raster, 0.01)
+```
+
 ## 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` of the raster or geometry on the right side. 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` 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`](../Predicate/#st_dwithin) instead.
+
+Rules for testing spatial relationship:
+
+- If the raster or geometry does not have a defined SRID, it is assumed to be in WGS84.
+- If both sides are in the same CRS, then the distance is measured in the units of that CRS.
+- Otherwise, both sides will be transformed to WGS84 before the distance test and the distance is measured in degrees.
+
+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 the radius, 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),
+    15.0
+) AS within_15_degrees
+```
+
+Output:
+
+```
++-----------------+
+|within_15_degrees|
++-----------------+
+|             true|
++-----------------+
+```
+
+Using `RS_DWithin` as a distance-join condition:
+
+```sql
+SELECT r.id, p.id
+FROM rasters r
+JOIN points p ON RS_DWithin(r.raster, p.geom, 0.01)
+```

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] =

spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/expressions/raster/RasterPredicates.scala

@@ -22,10 +22,10 @@ import org.apache.sedona.common.raster.RasterPredicates
 import org.apache.sedona.sql.utils.{GeometrySerializer, RasterSerializer}
 import org.apache.spark.sql.catalyst.InternalRow
 import org.apache.spark.sql.catalyst.expressions.codegen.CodegenFallback
-import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression}
+import org.apache.spark.sql.catalyst.expressions.{ExpectsInputTypes, Expression, ImplicitCastInputTypes}
 import org.apache.spark.sql.sedona_sql.UDT.{GeometryUDT, RasterUDT}
 import org.apache.spark.sql.sedona_sql.expressions.{FoldableExpression, NullIntolerantShim}
-import org.apache.spark.sql.types.{AbstractDataType, BooleanType, DataType}
+import org.apache.spark.sql.types.{AbstractDataType, BooleanType, DataType, DoubleType}
 import org.geotools.coverage.grid.GridCoverage2D
 import org.locationtech.jts.geom.Geometry
 
@@ -155,3 +155,83 @@ private[apache] case class RS_Within(inputExpressions: Seq[Expression])
     copy(inputExpressions = newChildren)
   }
 }
+
+/**
+ * Distance predicate for rasters: `RS_DWithin(left, right, distance)`. `left` and `right` can each
+ * be a raster or a geometry (at least one must be a raster). Returns true when the shapes are
+ * within `distance` of each other, with both sides projected to a common CRS prior to the JTS
+ * distance check (mirroring [[RS_Intersects]]). This expression is recognised by
+ * [[org.apache.spark.sql.sedona_sql.strategy.join.JoinQueryDetector]] as a distance-join key, so a
+ * coarse R-tree filter over expanded WGS84 envelopes prunes candidates before the per-row check.
+ */
+private[apache] case class RS_DWithin(inputExpressions: Seq[Expression])
+    extends Expression
+    with FoldableExpression
+    with ImplicitCastInputTypes
+    with NullIntolerantShim
+    with CodegenFallback {
+
+  override def children: Seq[Expression] = inputExpressions
+
+  override def nullable: Boolean = children.exists(_.nullable)
+
+  override def dataType: DataType = BooleanType
+
+  override def toString: String = s" **${getClass.getName}**  "
+
+  override def inputTypes: Seq[AbstractDataType] = {
+    if (inputExpressions.length != 3) {
+      throw new IllegalArgumentException(
+        s"RS_DWithin requires exactly 3 inputs, but got ${inputExpressions.length}")
+    }
+    val leftType = inputExpressions.head.dataType
+    val rightType = inputExpressions(1).dataType
+    (leftType, rightType) match {
+      case (_: RasterUDT, _: GeometryUDT) => Seq(RasterUDT(), GeometryUDT(), DoubleType)
+      case (_: GeometryUDT, _: RasterUDT) => Seq(GeometryUDT(), RasterUDT(), DoubleType)
+      case (_: RasterUDT, _: RasterUDT) => Seq(RasterUDT(), RasterUDT(), DoubleType)
+      case _ =>
+        throw new IllegalArgumentException(
+          s"RS_DWithin requires at least one raster input; got: $leftType, $rightType")
+    }
+  }
+
+  override def eval(inputRow: InternalRow): Any = {
+    val leftValue = inputExpressions.head.eval(inputRow)
+    if (leftValue == null) return null
+    val rightValue = inputExpressions(1).eval(inputRow)
+    if (rightValue == null) return null
+    val distanceValue = inputExpressions(2).eval(inputRow)
+    if (distanceValue == null) return null
+
+    val leftBytes = leftValue.asInstanceOf[Array[Byte]]
+    val rightBytes = rightValue.asInstanceOf[Array[Byte]]
+    val distance = distanceValue.asInstanceOf[Double]
+    val leftType = inputExpressions.head.dataType
+    val rightType = inputExpressions(1).dataType
+    (leftType, rightType) match {
+      case (_: RasterUDT, _: GeometryUDT) =>
+        RasterPredicates.rsDWithin(
+          RasterSerializer.deserialize(leftBytes),
+          GeometrySerializer.deserialize(rightBytes),
+          distance)
+      case (_: GeometryUDT, _: RasterUDT) =>
+        RasterPredicates.rsDWithin(
+          RasterSerializer.deserialize(rightBytes),
+          GeometrySerializer.deserialize(leftBytes),
+          distance)
+      case (_: RasterUDT, _: RasterUDT) =>
+        RasterPredicates.rsDWithin(
+          RasterSerializer.deserialize(leftBytes),
+          RasterSerializer.deserialize(rightBytes),
+          distance)
+      case _ =>
+        throw new IllegalArgumentException(
+          s"RS_DWithin requires at least one raster input; got: $leftType, $rightType")
+    }
+  }
+
+  protected def withNewChildrenInternal(newChildren: IndexedSeq[Expression]) = {
+    copy(inputExpressions = newChildren)
+  }
+}

spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/BroadcastIndexJoinExec.scala

@@ -122,8 +122,7 @@ case class BroadcastIndexJoinExec(
     case (None, _, false) => s"ST_$spatialPredicate($windowExpression, $objectExpression)"
     case (None, _, true) => s"RS_$spatialPredicate($windowExpression, $objectExpression)"
     case (Some(r), _, true) =>
-      throw new UnsupportedOperationException(
-        "Distance joins are not supported for raster predicates")
+      s"RS_Distance($windowExpression, $objectExpression) < $r"
   }
 
   override def simpleString(maxFields: Int): String =
@@ -321,6 +320,27 @@ case class BroadcastIndexJoinExec(
             (shape, row)
           }
         })
+      case Some(distanceExpression) if isRasterPredicate =>
+        val boundDistance =
+          BindReferences.bindReference(distanceExpression, streamed.output)
+        streamResultsRaw.map(row => {
+          val serialized = boundStreamShape.eval(row).asInstanceOf[Array[Byte]]
+          if (serialized == null) {
+            (null, row)
+          } else {
+            val baseShape = if (boundStreamShape.dataType.isInstanceOf[RasterUDT]) {
+              val raster = RasterSerializer.deserialize(serialized)
+              JoinedGeometryRaster.rasterToWGS84Envelope(raster)
+            } else {
+              val geom = GeometrySerializer.deserialize(serialized)
+              JoinedGeometryRaster.geometryToWGS84Envelope(geom)
+            }
+            val radius = boundDistance.eval(row).asInstanceOf[Double]
+            val envelope = baseShape.getEnvelopeInternal
+            envelope.expandBy(radius)
+            (baseShape.getFactory.toGeometry(envelope), row)
+          }
+        })
       case Some(distanceExpression) =>
         streamResultsRaw.map(row => {
           val geometry = TraitJoinQueryBase.shapeToGeometry(boundStreamShape, row)

spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/DistanceJoinExec.scala

@@ -24,6 +24,7 @@ import org.apache.spark.internal.Logging
 import org.apache.spark.rdd.RDD
 import org.apache.spark.sql.catalyst.expressions.{BindReferences, Expression, UnsafeRow}
 import org.apache.spark.sql.execution.SparkPlan
+import org.apache.spark.sql.sedona_sql.UDT.RasterUDT
 import org.apache.spark.sql.sedona_sql.execution.SedonaBinaryExecNode
 import org.locationtech.jts.geom.Geometry
 
@@ -83,7 +84,23 @@ case class DistanceJoinExec(
       leftShapeExpr: Expression,
       rightRdd: RDD[UnsafeRow],
       rightShapeExpr: Expression): (SpatialRDD[Geometry], SpatialRDD[Geometry]) = {
-    if (distanceBoundToLeft) {
+    // Raster predicates project both sides to a WGS84 envelope before partitioning so the coarse
+    // R-tree filter applies regardless of the input CRS — mirroring the non-distance raster join
+    // (toWGS84EnvelopeRDD) and the broadcast-index raster distance path.
+    val isRasterPredicate =
+      leftShapeExpr.dataType.isInstanceOf[RasterUDT] ||
+        rightShapeExpr.dataType.isInstanceOf[RasterUDT]
+    if (isRasterPredicate) {
+      if (distanceBoundToLeft) {
+        (
+          toExpandedWGS84EnvelopeRDD(leftRdd, leftShapeExpr, boundRadius),
+          toWGS84EnvelopeRDD(rightRdd, rightShapeExpr))
+      } else {
+        (
+          toWGS84EnvelopeRDD(leftRdd, leftShapeExpr),
+          toExpandedWGS84EnvelopeRDD(rightRdd, rightShapeExpr, boundRadius))
+      }
+    } else if (distanceBoundToLeft) {
       (
         toExpandedEnvelopeRDD(leftRdd, leftShapeExpr, boundRadius, isGeography),
         toSpatialRDD(rightRdd, rightShapeExpr))

spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/JoinQueryDetector.scala

@@ -174,6 +174,32 @@ class JoinQueryDetector(sparkSession: SparkSession) extends SparkStrategy {
         Some(condition)))
   }
 
+  /**
+   * Build a [[JoinQueryDetection]] for [[RS_DWithin]]. The coarse spatial join uses an expanded
+   * WGS84 envelope (driven by `distance`) and falls back to the per-row predicate for refinement,
+   * matching the pattern used for [[org.apache.spark.sql.sedona_sql.expressions.ST_DWithin]].
+   */
+  private def getRasterDistanceJoinDetection(
+      left: LogicalPlan,
+      right: LogicalPlan,
+      predicate: RS_DWithin,
+      extraCondition: Option[Expression] = None): Option[JoinQueryDetection] = {
+    val leftShape = predicate.inputExpressions.head
+    val rightShape = predicate.inputExpressions(1)
+    val distance = predicate.inputExpressions(2)
+    val condition = extraCondition.map(And(_, predicate)).getOrElse(predicate)
+    Some(
+      JoinQueryDetection(
+        left,
+        right,
+        leftShape,
+        rightShape,
+        SpatialPredicate.INTERSECTS,
+        isGeography = false,
+        extraCondition = Some(condition),
+        distance = Some(distance)))
+  }
+
   def apply(plan: LogicalPlan): Seq[SparkPlan] = plan match {
     case Join(left, right, joinType, condition, JoinHint(leftHint, rightHint))
         if optimizationEnabled(left, right, condition) =>
@@ -320,6 +346,8 @@ class JoinQueryDetector(sparkSession: SparkSession) extends SparkStrategy {
               getJoinDetection(left, right, pred, extraCondition)
             case pred: RS_Predicate =>
               getRasterJoinDetection(left, right, pred, extraCondition)
+            case pred: RS_DWithin =>
+              getRasterDistanceJoinDetection(left, right, pred, extraCondition)
             case ST_DWithin(Seq(leftShape, rightShape, distance)) =>
               val geographyShape =
                 isGeographyInput(leftShape) || isGeographyInput(rightShape)
@@ -903,9 +931,10 @@ class JoinQueryDetector(sparkSession: SparkSession) extends SparkStrategy {
         case (None, _, true, _, true) =>
           throw new UnsupportedOperationException(
             "Geography joins are not yet supported for raster predicates")
-        case (Some(_), _, _, _, true) =>
+        case (Some(_), _, true, _, true) =>
           throw new UnsupportedOperationException(
-            "Distance joins are not supported for raster predicates")
+            "Geography distance joins are not yet supported for raster predicates")
+        case (Some(_), _, false, _, true) => "RS_DWithin"
       }
     val (distanceOnIndexSide, distanceOnStreamSide) = distance
       .map { distanceExpr =>

spark/common/src/main/scala/org/apache/spark/sql/sedona_sql/strategy/join/OptimizableJoinCondition.scala

@@ -21,7 +21,7 @@ package org.apache.spark.sql.sedona_sql.strategy.join
 import org.apache.spark.sql.catalyst.expressions.{And, Expression, LessThan, LessThanOrEqual, Literal}
 import org.apache.spark.sql.catalyst.plans.logical.LogicalPlan
 import org.apache.spark.sql.sedona_sql.expressions._
-import org.apache.spark.sql.sedona_sql.expressions.raster.RS_Predicate
+import org.apache.spark.sql.sedona_sql.expressions.raster.{RS_DWithin, RS_Predicate}
 import org.apache.spark.sql.sedona_sql.optimization.ExpressionUtils
 
 case class OptimizableJoinCondition(left: LogicalPlan, right: LogicalPlan) {
@@ -73,6 +73,8 @@ case class OptimizableJoinCondition(left: LogicalPlan, right: LogicalPlan) {
       case ST_DWithin(Seq(leftShape, rightShape, distance, useSpheroid)) =>
         useSpheroid
           .isInstanceOf[Literal] && isDistanceJoinOptimizable(leftShape, rightShape, distance)
+      case RS_DWithin(Seq(leftShape, rightShape, distance)) =>
+        isDistanceJoinOptimizable(leftShape, rightShape, distance)
 
       case _: LessThan | _: LessThanOrEqual =>
         val (smaller, larger) = (expression.children.head, expression.children(1))