You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sedona.apache.org by ja...@apache.org on 2026/07/25 08:51:04 UTC
(sedona-db) branch main updated: feat(rust/sedona-spatial-join-raster): raster–vector spatial join (RS_Intersects/Contains/Within) (#1073)
This is an automated email from the ASF dual-hosted git repository.
james-willis pushed a commit to branch main
in repository https://gitbox.apache.org/repos/asf/sedona-db.git
The following commit(s) were added to refs/heads/main by this push:
new 3d07ddb41 feat(rust/sedona-spatial-join-raster): raster–vector spatial join (RS_Intersects/Contains/Within) (#1073)
3d07ddb41 is described below
commit 3d07ddb411367097e6eab0115f7186f7a00acdd0
Author: James Willis <ja...@users.noreply.github.com>
AuthorDate: Sat Jul 25 01:50:58 2026 -0700
feat(rust/sedona-spatial-join-raster): raster–vector spatial join (RS_Intersects/Contains/Within) (#1073)
---
Cargo.lock | 26 +
Cargo.toml | 2 +
docs/reference/sql/rs_intersects.qmd | 5 +
rust/sedona-functions/src/st_segmentize.rs | 19 +-
rust/sedona-geometry/src/interpolate.rs | 73 +++
rust/sedona-geometry/src/lib.rs | 1 +
.../sedona-query-planner/src/spatial_expr_utils.rs | 29 +-
rust/sedona-query-planner/src/spatial_predicate.rs | 43 +-
rust/sedona-raster-functions/src/footprint.rs | 165 ++++++
rust/sedona-raster-functions/src/lib.rs | 1 +
rust/sedona-raster-functions/src/rs_convexhull.rs | 29 +-
.../src/rs_spatial_predicates.rs | 65 ++-
rust/sedona-spatial-join-raster/Cargo.toml | 55 ++
.../src/join_provider.rs | 612 +++++++++++++++++++++
rust/sedona-spatial-join-raster/src/lib.rs | 36 ++
.../src/physical_planner.rs | 314 +++++++++++
.../tests/spatial_join_integration.rs | 599 ++++++++++++++++++++
rust/sedona/Cargo.toml | 3 +-
rust/sedona/src/context.rs | 15 +
rust/sedona/tests/raster_spatial_join.rs | 202 +++++++
20 files changed, 2221 insertions(+), 73 deletions(-)
diff --git a/Cargo.lock b/Cargo.lock
index 4818fa196..3ca1945d6 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -5763,6 +5763,7 @@ dependencies = [
"sedona-spatial-join",
"sedona-spatial-join-geography",
"sedona-spatial-join-gpu",
+ "sedona-spatial-join-raster",
"sedona-testing",
"sedona-tg",
"serde",
@@ -6439,6 +6440,31 @@ dependencies = [
"wkt 0.14.0",
]
+[[package]]
+name = "sedona-spatial-join-raster"
+version = "0.4.0"
+dependencies = [
+ "arrow-array",
+ "arrow-schema",
+ "datafusion",
+ "datafusion-common",
+ "datafusion-expr",
+ "datafusion-physical-expr",
+ "datafusion-physical-plan",
+ "rstest",
+ "sedona-common",
+ "sedona-expr",
+ "sedona-geometry",
+ "sedona-proj",
+ "sedona-query-planner",
+ "sedona-raster",
+ "sedona-raster-functions",
+ "sedona-schema",
+ "sedona-spatial-join",
+ "sedona-testing",
+ "tokio",
+]
+
[[package]]
name = "sedona-testing"
version = "0.4.0"
diff --git a/Cargo.toml b/Cargo.toml
index 5491848e4..9be2f203c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -45,6 +45,7 @@ members = [
"rust/sedona-spatial-join",
"rust/sedona-spatial-join-geography",
"rust/sedona-spatial-join-gpu",
+ "rust/sedona-spatial-join-raster",
"rust/sedona-query-planner",
"rust/sedona-testing",
"rust/sedona",
@@ -162,6 +163,7 @@ sedona-schema = { version = "0.4.0", path = "rust/sedona-schema" }
sedona-spatial-join = { version = "0.4.0", path = "rust/sedona-spatial-join" }
sedona-spatial-join-gpu = { version = "0.4.0", path = "rust/sedona-spatial-join-gpu" }
sedona-spatial-join-geography = { version = "0.4.0", path = "rust/sedona-spatial-join-geography" }
+sedona-spatial-join-raster = { version = "0.4.0", path = "rust/sedona-spatial-join-raster" }
sedona-query-planner = { version = "0.4.0", path = "rust/sedona-query-planner" }
sedona-testing = { version = "0.4.0", path = "rust/sedona-testing" }
diff --git a/docs/reference/sql/rs_intersects.qmd b/docs/reference/sql/rs_intersects.qmd
index f5633b980..755ddeef3 100644
--- a/docs/reference/sql/rs_intersects.qmd
+++ b/docs/reference/sql/rs_intersects.qmd
@@ -40,6 +40,11 @@ transformation fails, the extent of both sides are transformed to WGS 84 as a fa
If either argument has no CRS set, the comparison is performed directly without
CRS transformation.
+Across CRSes a reprojected raster footprint is densified with about 10 points
+per edge before it is reprojected, so it follows the curve of each edge instead
+of chording straight through the four reprojected corners. Same-CRS comparisons
+are exact.
+
## Examples
```sql
diff --git a/rust/sedona-functions/src/st_segmentize.rs b/rust/sedona-functions/src/st_segmentize.rs
index 5a8f1f15a..bced9bdad 100644
--- a/rust/sedona-functions/src/st_segmentize.rs
+++ b/rust/sedona-functions/src/st_segmentize.rs
@@ -31,6 +31,7 @@ use sedona_expr::{
};
use sedona_geometry::{
error::SedonaGeometryError,
+ interpolate::lerp,
wkb_factory::{
write_wkb_coord_trait, write_wkb_geometrycollection_header, write_wkb_linestring_header,
write_wkb_multilinestring_header, write_wkb_multipolygon_header, write_wkb_polygon_header,
@@ -312,8 +313,8 @@ fn write_interpolated_2d<C: CoordTrait<T = f64>>(
for i in 1..=num_segments {
let t = i as f64 / num_segments as f64;
- writer.write_all(&(x1 + t * (x2 - x1)).to_le_bytes())?;
- writer.write_all(&(y1 + t * (y2 - y1)).to_le_bytes())?;
+ writer.write_all(&lerp(x1, x2, t).to_le_bytes())?;
+ writer.write_all(&lerp(y1, y2, t).to_le_bytes())?;
}
Ok(())
}
@@ -334,9 +335,9 @@ fn write_interpolated_3d<C: CoordTrait<T = f64>>(
for i in 1..=num_segments {
let t = i as f64 / num_segments as f64;
- writer.write_all(&(x1 + t * (x2 - x1)).to_le_bytes())?;
- writer.write_all(&(y1 + t * (y2 - y1)).to_le_bytes())?;
- writer.write_all(&(d1 + t * (d2 - d1)).to_le_bytes())?;
+ writer.write_all(&lerp(x1, x2, t).to_le_bytes())?;
+ writer.write_all(&lerp(y1, y2, t).to_le_bytes())?;
+ writer.write_all(&lerp(d1, d2, t).to_le_bytes())?;
}
Ok(())
}
@@ -359,10 +360,10 @@ fn write_interpolated_4d<C: CoordTrait<T = f64>>(
for i in 1..=num_segments {
let t = i as f64 / num_segments as f64;
- writer.write_all(&(x1 + t * (x2 - x1)).to_le_bytes())?;
- writer.write_all(&(y1 + t * (y2 - y1)).to_le_bytes())?;
- writer.write_all(&(z1 + t * (z2 - z1)).to_le_bytes())?;
- writer.write_all(&(m1 + t * (m2 - m1)).to_le_bytes())?;
+ writer.write_all(&lerp(x1, x2, t).to_le_bytes())?;
+ writer.write_all(&lerp(y1, y2, t).to_le_bytes())?;
+ writer.write_all(&lerp(z1, z2, t).to_le_bytes())?;
+ writer.write_all(&lerp(m1, m2, t).to_le_bytes())?;
}
Ok(())
}
diff --git a/rust/sedona-geometry/src/interpolate.rs b/rust/sedona-geometry/src/interpolate.rs
new file mode 100644
index 000000000..efd7b7f53
--- /dev/null
+++ b/rust/sedona-geometry/src/interpolate.rs
@@ -0,0 +1,73 @@
+// 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.
+
+//! Linear interpolation helpers shared by geometry densification paths
+//! (`ST_Segmentize`, reprojected raster footprints).
+
+/// Linearly interpolate between `a` and `b` at parameter `t`.
+///
+/// `lerp(a, b, 0.0) == a` and `lerp(a, b, 1.0) == b`. This is the single-scalar
+/// primitive behind every densification path so the interpolation math is not
+/// re-derived per caller (per-dimension `ST_Segmentize` interpolation, per-axis
+/// footprint densification).
+#[inline]
+pub fn lerp(a: f64, b: f64, t: f64) -> f64 {
+ a + t * (b - a)
+}
+
+/// Append `n` evenly spaced interior points of the segment from `start` to `end`
+/// to `out` — the points at parameters `t = 1/(n+1), 2/(n+1), …, n/(n+1)`.
+///
+/// Neither endpoint is emitted; the caller owns `start` and `end`. Densifying a
+/// polygon edge this way makes the polygon follow the curved image of that
+/// straight edge under a nonlinear reprojection instead of chording across it.
+/// `out` is appended to (not cleared) so a ring can be assembled edge by edge
+/// into a caller-managed scratch buffer.
+pub fn densify_segment(start: (f64, f64), end: (f64, f64), n: usize, out: &mut Vec<(f64, f64)>) {
+ for i in 1..=n {
+ let t = i as f64 / (n + 1) as f64;
+ out.push((lerp(start.0, end.0, t), lerp(start.1, end.1, t)));
+ }
+}
+
+#[cfg(test)]
+mod test {
+ use super::*;
+
+ #[test]
+ fn lerp_hits_endpoints_and_midpoint() {
+ assert_eq!(lerp(2.0, 6.0, 0.0), 2.0);
+ assert_eq!(lerp(2.0, 6.0, 1.0), 6.0);
+ assert_eq!(lerp(2.0, 6.0, 0.5), 4.0);
+ assert_eq!(lerp(2.0, 6.0, 0.25), 3.0);
+ }
+
+ #[test]
+ fn densify_segment_emits_evenly_spaced_interior_points() {
+ let mut out = Vec::new();
+ densify_segment((0.0, 0.0), (4.0, 8.0), 3, &mut out);
+ // Interior points at t = 1/4, 2/4, 3/4 (endpoints excluded).
+ assert_eq!(out, vec![(1.0, 2.0), (2.0, 4.0), (3.0, 6.0)]);
+ }
+
+ #[test]
+ fn densify_segment_zero_points_is_noop() {
+ let mut out = vec![(9.0, 9.0)];
+ densify_segment((0.0, 0.0), (1.0, 1.0), 0, &mut out);
+ assert_eq!(out, vec![(9.0, 9.0)]);
+ }
+}
diff --git a/rust/sedona-geometry/src/lib.rs b/rust/sedona-geometry/src/lib.rs
index 5de9f65e2..b50ab36a9 100644
--- a/rust/sedona-geometry/src/lib.rs
+++ b/rust/sedona-geometry/src/lib.rs
@@ -19,6 +19,7 @@ pub mod bounding_box;
pub mod bounds;
pub mod error;
pub mod ewkb_factory;
+pub mod interpolate;
pub mod interval;
pub mod is_empty;
pub mod point_count;
diff --git a/rust/sedona-query-planner/src/spatial_expr_utils.rs b/rust/sedona-query-planner/src/spatial_expr_utils.rs
index 0aa581791..32d1d5b3c 100644
--- a/rust/sedona-query-planner/src/spatial_expr_utils.rs
+++ b/rust/sedona-query-planner/src/spatial_expr_utils.rs
@@ -71,6 +71,9 @@ pub(crate) fn collect_spatial_predicate_names(expr: &Expr) -> HashSet<String> {
| "st_equals"
| "st_dwithin"
| "st_knn"
+ | "rs_intersects"
+ | "rs_contains"
+ | "rs_within"
) {
acc.insert(func_name);
}
@@ -627,7 +630,7 @@ mod tests {
use datafusion_physical_expr::{PhysicalExpr, ScalarFunctionExpr};
use datafusion_physical_plan::joins::utils::ColumnIndex;
use datafusion_physical_plan::joins::utils::JoinFilter;
- use sedona_schema::datatypes::WKB_GEOMETRY;
+ use sedona_schema::datatypes::{RASTER, WKB_GEOMETRY};
use std::sync::Arc;
// Helper function to create a test schema
@@ -2295,6 +2298,30 @@ mod tests {
assert!(names.is_empty());
}
+ #[test]
+ fn test_collect_spatial_predicate_names_raster() {
+ // Raster-geometry predicates are recognized as spatial predicates so their joins
+ // route through the spatial-join planner (and fall back to nested-loop join).
+ for name in ["rs_intersects", "rs_contains", "rs_within"] {
+ let udf = Arc::new(ScalarUDF::from(SimpleScalarUDF::new(
+ name,
+ vec![
+ RASTER.storage_type().clone(),
+ WKB_GEOMETRY.storage_type().clone(),
+ ],
+ DataType::Boolean,
+ datafusion_expr::Volatility::Immutable,
+ Arc::new(|_| Ok(ColumnarValue::Scalar(ScalarValue::Boolean(Some(true))))),
+ )));
+ let expr = Expr::ScalarFunction(datafusion_expr::expr::ScalarFunction {
+ func: udf,
+ args: vec![col("raster"), col("geom")],
+ });
+ let names = collect_spatial_predicate_names(&expr);
+ assert_eq!(names, HashSet::from([name.to_string()]));
+ }
+ }
+
#[test]
fn test_find_st_knn_call_simple() {
let st_knn_udf = create_dummy_st_knn_udf();
diff --git a/rust/sedona-query-planner/src/spatial_predicate.rs b/rust/sedona-query-planner/src/spatial_predicate.rs
index 63ccc4718..0a9605364 100644
--- a/rust/sedona-query-planner/src/spatial_predicate.rs
+++ b/rust/sedona-query-planner/src/spatial_predicate.rs
@@ -190,16 +190,19 @@ impl SpatialRelationType {
/// Converts a function name string to a SpatialRelationType.
///
/// # Arguments
- /// * `name` - The spatial function name (e.g., "st_intersects", "st_contains")
+ /// * `name` - The spatial function name (e.g., "st_intersects", "st_contains").
+ /// Raster-geometry predicates (`rs_intersects`, `rs_contains`, `rs_within`)
+ /// map to the same relation variants; `invert()` handles the Contains/Within
+ /// swap when the join arguments are on opposite sides.
///
/// # Returns
/// * `Some(SpatialRelationType)` if the name is recognized
/// * `None` if the name is not a valid spatial relation function
pub fn from_name(name: &str) -> Option<Self> {
match name {
- "st_intersects" => Some(SpatialRelationType::Intersects),
- "st_contains" => Some(SpatialRelationType::Contains),
- "st_within" => Some(SpatialRelationType::Within),
+ "st_intersects" | "rs_intersects" => Some(SpatialRelationType::Intersects),
+ "st_contains" | "rs_contains" => Some(SpatialRelationType::Contains),
+ "st_within" | "rs_within" => Some(SpatialRelationType::Within),
"st_covers" => Some(SpatialRelationType::Covers),
"st_coveredby" | "st_covered_by" => Some(SpatialRelationType::CoveredBy),
"st_touches" => Some(SpatialRelationType::Touches),
@@ -531,6 +534,38 @@ mod tests {
assert_eq!(col.index(), index);
}
+ #[test]
+ fn from_name_recognizes_raster_predicates() {
+ // Raster-geometry predicates reuse the geometry relation variants; `invert()`
+ // handles the Contains/Within swap when the join arguments are on opposite sides.
+ assert_eq!(
+ SpatialRelationType::from_name("rs_intersects"),
+ Some(SpatialRelationType::Intersects)
+ );
+ assert_eq!(
+ SpatialRelationType::from_name("rs_contains"),
+ Some(SpatialRelationType::Contains)
+ );
+ assert_eq!(
+ SpatialRelationType::from_name("rs_within"),
+ Some(SpatialRelationType::Within)
+ );
+
+ // Matches the ST_ variants they mirror.
+ assert_eq!(
+ SpatialRelationType::from_name("rs_intersects"),
+ SpatialRelationType::from_name("st_intersects")
+ );
+ assert_eq!(
+ SpatialRelationType::from_name("rs_contains"),
+ SpatialRelationType::from_name("st_contains")
+ );
+ assert_eq!(
+ SpatialRelationType::from_name("rs_within"),
+ SpatialRelationType::from_name("st_within")
+ );
+ }
+
#[test]
fn relation_rewrite_success() -> Result<()> {
let on = SpatialPredicate::Relation(RelationPredicate {
diff --git a/rust/sedona-raster-functions/src/footprint.rs b/rust/sedona-raster-functions/src/footprint.rs
new file mode 100644
index 000000000..208e75626
--- /dev/null
+++ b/rust/sedona-raster-functions/src/footprint.rs
@@ -0,0 +1,165 @@
+// 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.
+
+//! Raster footprint helpers shared by the raster spatial-predicate kernels and
+//! the optimized raster spatial join.
+//!
+//! A raster's footprint is the convex hull of its four corners in world
+//! coordinates. Because the affine geotransform may include skew/rotation, each
+//! corner is computed individually rather than assumed axis-aligned.
+
+use datafusion_common::{DataFusionError, Result};
+use sedona_geometry::interpolate::densify_segment;
+use sedona_geometry::wkb_factory::write_wkb_polygon;
+use sedona_raster::affine_transformation::to_world_coordinate;
+use sedona_raster::traits::RasterRef;
+
+/// Number of interior points interpolated along each footprint edge when a
+/// footprint is densified before reprojection.
+///
+/// Reprojecting a footprint bends its straight edges into curves, so a footprint
+/// described only by its four corners chords across each curve and under-/over-
+/// covers the true extent. Densifying each edge with a handful of interior points
+/// (~10 per edge is a common, cheap choice) and reprojecting every point makes
+/// the reprojected footprint follow the curve. A footprint left in its own CRS
+/// has exact straight edges and needs no densification.
+pub const FOOTPRINT_POINTS_PER_EDGE: usize = 10;
+
+/// The four corners of a raster's footprint in world coordinates.
+///
+/// Returned in ring order: upper-left `(0, 0)`, upper-right `(width, 0)`,
+/// lower-right `(width, height)`, lower-left `(0, height)`.
+pub fn raster_footprint_corners(raster: &dyn RasterRef) -> [(f64, f64); 4] {
+ let width = raster.metadata().width();
+ let height = raster.metadata().height();
+
+ [
+ to_world_coordinate(raster, 0, 0),
+ to_world_coordinate(raster, width, 0),
+ to_world_coordinate(raster, width, height),
+ to_world_coordinate(raster, 0, height),
+ ]
+}
+
+/// Write WKB for the convex-hull polygon through four footprint `corners`.
+///
+/// `corners` are in ring order (upper-left, upper-right, lower-right,
+/// lower-left, as produced by [`raster_footprint_corners`]); the ring is closed
+/// back to the first corner. Shared by the native footprint (corners in the
+/// raster's own CRS) and the reprojected footprint (corners transformed into
+/// another CRS), so both paths emit byte-identical polygon WKB. This can be used
+/// to build Binary arrays, as the arrow-rs `BinaryBuilder` implements
+/// [`std::io::Write`].
+pub fn write_footprint_wkb(corners: [(f64, f64); 4], out: &mut impl std::io::Write) -> Result<()> {
+ let [ul, ur, lr, ll] = corners;
+
+ write_wkb_polygon(out, [ul, ur, lr, ll, ul].into_iter())
+ .map_err(|e| DataFusionError::External(e.into()))?;
+
+ Ok(())
+}
+
+/// Write WKB for the convex-hull polygon of the raster footprint.
+///
+/// The ring is the four [`raster_footprint_corners`] closed back to the
+/// upper-left corner.
+pub fn write_convexhull_wkb(raster: &dyn RasterRef, out: &mut impl std::io::Write) -> Result<()> {
+ write_footprint_wkb(raster_footprint_corners(raster), out)
+}
+
+/// Fill `ring` with a closed, densified footprint ring through four `corners`.
+///
+/// `corners` are in ring order (upper-left, upper-right, lower-right, lower-left,
+/// as produced by [`raster_footprint_corners`]). Each edge contributes its start
+/// corner followed by [`FOOTPRINT_POINTS_PER_EDGE`] interior points, and the ring
+/// is closed back to the upper-left corner. `ring` is cleared first and is a
+/// caller-managed scratch buffer so densifying many rasters allocates once.
+///
+/// Densify in the CRS where the footprint has exact straight edges (its own),
+/// then reproject every point of the returned ring: the reprojected ring then
+/// follows the curved image of each edge instead of chording across it.
+pub fn densify_footprint_ring(corners: [(f64, f64); 4], ring: &mut Vec<(f64, f64)>) {
+ ring.clear();
+ let [ul, ur, lr, ll] = corners;
+ for (start, end) in [(ul, ur), (ur, lr), (lr, ll), (ll, ul)] {
+ ring.push(start);
+ densify_segment(start, end, FOOTPRINT_POINTS_PER_EDGE, ring);
+ }
+ ring.push(ul);
+}
+
+/// Write WKB for the polygon through a pre-closed footprint `ring` (last point
+/// equal to the first), as produced by [`densify_footprint_ring`].
+///
+/// Companion to [`write_footprint_wkb`] for the variable-length densified ring;
+/// both emit a single-ring polygon and can build Binary arrays, as the arrow-rs
+/// `BinaryBuilder` implements [`std::io::Write`].
+pub fn write_footprint_ring_wkb(ring: &[(f64, f64)], out: &mut impl std::io::Write) -> Result<()> {
+ write_wkb_polygon(out, ring.iter().copied())
+ .map_err(|e| DataFusionError::External(e.into()))?;
+
+ Ok(())
+}
+
+#[cfg(test)]
+mod test {
+ use super::*;
+
+ /// Index of each corner in a densified ring: the start of edge `k` sits after
+ /// `k` corners and `k * FOOTPRINT_POINTS_PER_EDGE` interior points.
+ fn corner_index(k: usize) -> usize {
+ k * (FOOTPRINT_POINTS_PER_EDGE + 1)
+ }
+
+ #[test]
+ fn densify_footprint_ring_layout_and_coords() {
+ // Axis-aligned 10x10 square in ring order (UL, UR, LR, LL).
+ let corners = [(0.0, 10.0), (10.0, 10.0), (10.0, 0.0), (0.0, 0.0)];
+ let mut ring = Vec::new();
+ densify_footprint_ring(corners, &mut ring);
+
+ // 4 corners + 4 * FOOTPRINT_POINTS_PER_EDGE interior + 1 closing point.
+ assert_eq!(ring.len(), 4 * (FOOTPRINT_POINTS_PER_EDGE + 1) + 1);
+
+ // Corners land exactly on the interpolated ring, in order, and the ring
+ // closes back to the upper-left corner.
+ assert_eq!(ring[corner_index(0)], corners[0]);
+ assert_eq!(ring[corner_index(1)], corners[1]);
+ assert_eq!(ring[corner_index(2)], corners[2]);
+ assert_eq!(ring[corner_index(3)], corners[3]);
+ assert_eq!(*ring.last().unwrap(), corners[0]);
+
+ // The first edge (UL -> UR) runs along y = 10 with x strictly increasing
+ // through the interior points, spaced at t = i / (points_per_edge + 1).
+ for i in 1..=FOOTPRINT_POINTS_PER_EDGE {
+ let (x, y) = ring[corner_index(0) + i];
+ let t = i as f64 / (FOOTPRINT_POINTS_PER_EDGE + 1) as f64;
+ assert!((x - 10.0 * t).abs() < 1e-12, "x={x} t={t}");
+ assert_eq!(y, 10.0);
+ }
+ }
+
+ #[test]
+ fn densify_footprint_ring_reuses_scratch() {
+ let mut ring = vec![(42.0, 42.0)];
+ let corners = [(0.0, 1.0), (1.0, 1.0), (1.0, 0.0), (0.0, 0.0)];
+ densify_footprint_ring(corners, &mut ring);
+ // Pre-existing scratch content is cleared, not appended to.
+ assert_eq!(ring.len(), 4 * (FOOTPRINT_POINTS_PER_EDGE + 1) + 1);
+ assert_eq!(ring[0], corners[0]);
+ }
+}
diff --git a/rust/sedona-raster-functions/src/lib.rs b/rust/sedona-raster-functions/src/lib.rs
index b650fee01..aa53e4d48 100644
--- a/rust/sedona-raster-functions/src/lib.rs
+++ b/rust/sedona-raster-functions/src/lib.rs
@@ -18,6 +18,7 @@
pub mod crs_utils;
mod executor;
pub use executor::RasterExecutor;
+pub mod footprint;
pub mod register;
pub mod rs_band_accessors;
pub mod rs_bandpath;
diff --git a/rust/sedona-raster-functions/src/rs_convexhull.rs b/rust/sedona-raster-functions/src/rs_convexhull.rs
index 172c84d4b..c6bcf1dc6 100644
--- a/rust/sedona-raster-functions/src/rs_convexhull.rs
+++ b/rust/sedona-raster-functions/src/rs_convexhull.rs
@@ -17,16 +17,14 @@
use std::sync::Arc;
use crate::executor::RasterExecutor;
+use crate::footprint::write_convexhull_wkb;
use arrow_array::builder::{BinaryBuilder, StringViewBuilder};
-use datafusion_common::DataFusionError;
use datafusion_common::Result;
use datafusion_common::ScalarValue;
use datafusion_expr::{ColumnarValue, Volatility};
use sedona_expr::item_crs::make_item_crs;
use sedona_expr::scalar_udf::{SedonaScalarKernel, SedonaScalarUDF};
use sedona_geometry::types::Edges;
-use sedona_geometry::wkb_factory::write_wkb_polygon;
-use sedona_raster::affine_transformation::to_world_coordinate;
use sedona_raster::traits::RasterRef;
use sedona_schema::{datatypes::SedonaType, matchers::ArgMatcher};
@@ -101,31 +99,6 @@ impl SedonaScalarKernel for RsConvexHull {
}
}
-/// Write WKB for a convex hull polygon for the raster
-///
-/// For a raster, the convex hull is the polygon formed by the four corners
-/// of the raster in world coordinates. Due to skew/rotation in the affine
-/// transformation, each corner must be computed individually.
-fn write_convexhull_wkb(raster: &dyn RasterRef, out: &mut impl std::io::Write) -> Result<()> {
- let width = raster.metadata().width();
- let height = raster.metadata().height();
-
- // Compute the four corners in pixel coordinates:
- // Upper-left (0, 0), Upper-right (width, 0), Lower-right (width, height), Lower-left (0, height)
- let (ulx, uly) = to_world_coordinate(raster, 0, 0);
- let (urx, ury) = to_world_coordinate(raster, width, 0);
- let (lrx, lry) = to_world_coordinate(raster, width, height);
- let (llx, lly) = to_world_coordinate(raster, 0, height);
-
- write_wkb_polygon(
- out,
- [(ulx, uly), (urx, ury), (lrx, lry), (llx, lly), (ulx, uly)].into_iter(),
- )
- .map_err(|e| DataFusionError::External(e.into()))?;
-
- Ok(())
-}
-
#[cfg(test)]
mod tests {
use super::*;
diff --git a/rust/sedona-raster-functions/src/rs_spatial_predicates.rs b/rust/sedona-raster-functions/src/rs_spatial_predicates.rs
index 7c5650f9c..5e8b0ab21 100644
--- a/rust/sedona-raster-functions/src/rs_spatial_predicates.rs
+++ b/rust/sedona-raster-functions/src/rs_spatial_predicates.rs
@@ -33,19 +33,18 @@ use std::sync::Arc;
use crate::crs_utils::crs_transform_wkb;
use crate::crs_utils::resolve_crs;
+use crate::crs_utils::with_crs_engine;
use crate::executor::RasterExecutor;
+use crate::footprint::write_convexhull_wkb;
use arrow_array::builder::BooleanBuilder;
use arrow_schema::DataType;
+use datafusion_common::config::ConfigOptions;
use datafusion_common::exec_datafusion_err;
use datafusion_common::exec_err;
-use datafusion_common::DataFusionError;
use datafusion_common::Result;
use datafusion_expr::{ColumnarValue, Volatility};
use sedona_expr::scalar_udf::{SedonaScalarKernel, SedonaScalarUDF};
use sedona_geometry::transform::CrsEngine;
-use sedona_geometry::wkb_factory::write_wkb_polygon;
-use sedona_proj::transform::with_global_proj_engine;
-use sedona_raster::affine_transformation::to_world_coordinate;
use sedona_raster::traits::RasterRef;
use sedona_schema::crs::{lnglat, CrsRef};
use sedona_schema::{datatypes::SedonaType, matchers::ArgMatcher};
@@ -166,11 +165,28 @@ impl<Op: tg::BinaryPredicate + Send + Sync> SedonaScalarKernel for RsSpatialPred
&self,
arg_types: &[SedonaType],
args: &[ColumnarValue],
+ ) -> Result<ColumnarValue> {
+ self.invoke_batch_from_args(
+ arg_types,
+ args,
+ &SedonaType::Arrow(DataType::Boolean),
+ 0,
+ None,
+ )
+ }
+
+ fn invoke_batch_from_args(
+ &self,
+ arg_types: &[SedonaType],
+ args: &[ColumnarValue],
+ _return_type: &SedonaType,
+ _num_rows: usize,
+ config_options: Option<&ConfigOptions>,
) -> Result<ColumnarValue> {
match self.arg_order {
- ArgOrder::RasterGeom => self.invoke_raster_geom(arg_types, args),
- ArgOrder::GeomRaster => self.invoke_geom_raster(arg_types, args),
- ArgOrder::RasterRaster => self.invoke_raster_raster(arg_types, args),
+ ArgOrder::RasterGeom => self.invoke_raster_geom(arg_types, args, config_options),
+ ArgOrder::GeomRaster => self.invoke_geom_raster(arg_types, args, config_options),
+ ArgOrder::RasterRaster => self.invoke_raster_raster(arg_types, args, config_options),
}
}
}
@@ -181,6 +197,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
&self,
arg_types: &[SedonaType],
args: &[ColumnarValue],
+ config_options: Option<&ConfigOptions>,
) -> Result<ColumnarValue> {
// Ensure executor always sees (raster, geom)
let exec_arg_types = vec![arg_types[0].clone(), arg_types[1].clone()];
@@ -189,7 +206,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
let mut builder = BooleanBuilder::with_capacity(executor.num_iterations());
let mut raster_wkb = Vec::with_capacity(CONVEXHULL_WKB_SIZE);
- with_global_proj_engine(|engine| {
+ with_crs_engine(config_options, |engine| {
executor.execute_raster_wkb_crs_void(|raster_opt, maybe_wkb, geom_crs| {
match (raster_opt, maybe_wkb) {
(Some(raster), Some(geom_wkb)) => {
@@ -221,6 +238,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
&self,
arg_types: &[SedonaType],
args: &[ColumnarValue],
+ config_options: Option<&ConfigOptions>,
) -> Result<ColumnarValue> {
// Reorder so executor always sees (raster, geom)
let exec_arg_types = vec![arg_types[1].clone(), arg_types[0].clone()];
@@ -229,7 +247,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
let mut builder = BooleanBuilder::with_capacity(executor.num_iterations());
let mut raster_wkb = Vec::with_capacity(CONVEXHULL_WKB_SIZE);
- with_global_proj_engine(|engine| {
+ with_crs_engine(config_options, |engine| {
executor.execute_raster_wkb_crs_void(|raster_opt, maybe_wkb, geom_crs| {
match (raster_opt, maybe_wkb) {
(Some(raster), Some(geom_wkb)) => {
@@ -262,6 +280,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
&self,
arg_types: &[SedonaType],
args: &[ColumnarValue],
+ config_options: Option<&ConfigOptions>,
) -> Result<ColumnarValue> {
// Ensure executor always sees (raster, raster)
let exec_arg_types = vec![arg_types[0].clone(), arg_types[1].clone()];
@@ -271,7 +290,7 @@ impl<Op: tg::BinaryPredicate + Send + Sync> RsSpatialPredicate<Op> {
let mut wkb0 = Vec::with_capacity(CONVEXHULL_WKB_SIZE);
let mut wkb1 = Vec::with_capacity(CONVEXHULL_WKB_SIZE);
- with_global_proj_engine(|engine| {
+ with_crs_engine(config_options, |engine| {
executor.execute_raster_raster_void(|_i, r0_opt, r1_opt| {
match (r0_opt, r1_opt) {
(Some(r0), Some(r1)) => {
@@ -375,32 +394,15 @@ fn evaluate_predicate<Op: tg::BinaryPredicate>(wkb_a: &[u8], wkb_b: &[u8]) -> Re
/// = 9 + 4 + 80 = 93
const CONVEXHULL_WKB_SIZE: usize = 93;
-/// Create WKB for a convex hull polygon for the raster
-fn write_convexhull_wkb(raster: &dyn RasterRef, out: &mut impl std::io::Write) -> Result<()> {
- let width = raster.metadata().width();
- let height = raster.metadata().height();
-
- let (ulx, uly) = to_world_coordinate(raster, 0, 0);
- let (urx, ury) = to_world_coordinate(raster, width, 0);
- let (lrx, lry) = to_world_coordinate(raster, width, height);
- let (llx, lly) = to_world_coordinate(raster, 0, height);
-
- write_wkb_polygon(
- out,
- [(ulx, uly), (urx, ury), (lrx, lry), (llx, lly), (ulx, uly)].into_iter(),
- )
- .map_err(|e| DataFusionError::External(e.into()))?;
-
- Ok(())
-}
-
#[cfg(test)]
mod tests {
use super::*;
use arrow_array::{create_array, ArrayRef};
+ use datafusion_common::DataFusionError;
use datafusion_expr::ScalarUDF;
use rstest::rstest;
use sedona_geometry::types::Edges;
+ use sedona_proj::transform::{with_global_proj_engine, LazyProjEngine};
use sedona_raster::builder::RasterBuilder;
use sedona_raster::traits::{BandMetadata, RasterMetadata};
use sedona_schema::crs::deserialize_crs;
@@ -515,7 +517,10 @@ mod tests {
fn rs_intersects_raster_geom_crs_mismatch() {
let udf = rs_intersects_udf();
let geom_type = SedonaType::Wkb(Edges::Planar, deserialize_crs("EPSG:3857").unwrap());
- let tester = ScalarUdfTester::new(udf.into(), vec![RASTER, geom_type.clone()]);
+ // This overload reprojects (EPSG:3857 geometry vs the raster's OGC:CRS84),
+ // so the session must carry a real CRS engine.
+ let tester = ScalarUdfTester::new(udf.into(), vec![RASTER, geom_type.clone()])
+ .with_crs_engine(Arc::new(LazyProjEngine));
let rasters = generate_test_rasters(3, Some(0)).unwrap();
let (x, y) = with_global_proj_engine(|engine| {
diff --git a/rust/sedona-spatial-join-raster/Cargo.toml b/rust/sedona-spatial-join-raster/Cargo.toml
new file mode 100644
index 000000000..5862c4c10
--- /dev/null
+++ b/rust/sedona-spatial-join-raster/Cargo.toml
@@ -0,0 +1,55 @@
+# 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]
+name = "sedona-spatial-join-raster"
+version.workspace = true
+license.workspace = true
+keywords.workspace = true
+categories.workspace = true
+homepage.workspace = true
+repository.workspace = true
+description.workspace = true
+readme.workspace = true
+edition.workspace = true
+rust-version.workspace = true
+
+[dependencies]
+arrow-array = { workspace = true }
+arrow-schema = { workspace = true }
+datafusion-common = { workspace = true }
+datafusion-expr = { workspace = true }
+datafusion-physical-expr = { workspace = true }
+datafusion-physical-plan = { workspace = true }
+sedona-common = { workspace = true }
+sedona-expr = { workspace = true }
+sedona-geometry = { workspace = true }
+sedona-proj = { workspace = true }
+sedona-query-planner = { workspace = true }
+sedona-raster = { workspace = true }
+sedona-raster-functions = { workspace = true }
+sedona-schema = { workspace = true }
+sedona-spatial-join = { workspace = true }
+
+[dev-dependencies]
+arrow-array = { workspace = true }
+datafusion = { workspace = true, features = ["sql"] }
+datafusion-common = { workspace = true }
+datafusion-expr = { workspace = true }
+rstest = { workspace = true }
+sedona-proj = { workspace = true, features = ["proj-sys"] }
+sedona-testing = { workspace = true }
+tokio = { workspace = true, features = ["macros", "rt-multi-thread"] }
diff --git a/rust/sedona-spatial-join-raster/src/join_provider.rs b/rust/sedona-spatial-join-raster/src/join_provider.rs
new file mode 100644
index 000000000..7499b1ba4
--- /dev/null
+++ b/rust/sedona-spatial-join-raster/src/join_provider.rs
@@ -0,0 +1,612 @@
+// 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.
+
+use std::sync::Arc;
+
+use arrow_array::{builder::BinaryBuilder, Array, ArrayRef, StructArray};
+use datafusion_common::{exec_datafusion_err, exec_err, JoinType, Result};
+use datafusion_expr::ColumnarValue;
+use sedona_common::{
+ sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions, SpatialLibrary,
+};
+use sedona_expr::statistics::GeoStatistics;
+use sedona_geometry::{
+ interval::{Interval, IntervalTrait},
+ transform::CrsEngine,
+};
+use sedona_raster::array::RasterStructArray;
+use sedona_raster::traits::RasterRef;
+use sedona_raster_functions::crs_utils::resolve_crs;
+use sedona_raster_functions::footprint::{
+ densify_footprint_ring, raster_footprint_corners, write_footprint_ring_wkb, write_footprint_wkb,
+};
+use sedona_schema::{
+ crs::{CoordinateReferenceSystem, Crs},
+ datatypes::SedonaType,
+ datatypes::WKB_GEOMETRY,
+};
+use sedona_spatial_join::{
+ index::{spatial_index_builder::SpatialJoinBuildMetrics, SpatialIndexBuilder},
+ join_provider::{DefaultSpatialJoinProvider, SpatialJoinProvider},
+ operand_evaluator::{EvaluatedGeometryArray, EvaluatedGeometryArrayFactory},
+ utils::bounds::Bounds2D,
+ SpatialPredicate,
+};
+
+/// [`SpatialJoinProvider`] for raster/geometry spatial joins.
+///
+/// The R-tree index builder and the memory estimate delegate to the default
+/// provider; only the operand evaluator is raster-aware. The factory it produces
+/// reprojects raster footprints into `target_crs` (the geometry operand's CRS)
+/// using the session's [`CrsEngine`], captured at plan time.
+#[derive(Debug)]
+pub(crate) struct RasterJoinProvider {
+ default: DefaultSpatialJoinProvider,
+ target_crs: Crs,
+ engine: Arc<dyn CrsEngine + Send + Sync>,
+}
+
+impl RasterJoinProvider {
+ pub(crate) fn new(target_crs: Crs, engine: Arc<dyn CrsEngine + Send + Sync>) -> Self {
+ Self {
+ default: DefaultSpatialJoinProvider,
+ target_crs,
+ engine,
+ }
+ }
+}
+
+/// Pin the join's refiner to `tg`, the engine the `RS_*` predicate kernel uses,
+/// so the accelerated join and the kernel resolve boundary/touching cases
+/// identically regardless of the session's `spatial_join.spatial_library`.
+fn pin_refiner_to_tg(mut options: SpatialJoinOptions) -> SpatialJoinOptions {
+ options.spatial_library = SpatialLibrary::Tg;
+ options
+}
+
+impl SpatialJoinProvider for RasterJoinProvider {
+ fn try_new_spatial_index_builder(
+ &self,
+ schema: arrow_schema::SchemaRef,
+ spatial_predicate: SpatialPredicate,
+ options: SpatialJoinOptions,
+ join_type: JoinType,
+ probe_threads_count: usize,
+ metrics: SpatialJoinBuildMetrics,
+ ) -> Result<Box<dyn SpatialIndexBuilder>> {
+ // Footprints are ordinary planar WKB polygons, so the default R-tree
+ // builder and WKB refiner apply unchanged.
+ self.default.try_new_spatial_index_builder(
+ schema,
+ spatial_predicate,
+ pin_refiner_to_tg(options),
+ join_type,
+ probe_threads_count,
+ metrics,
+ )
+ }
+
+ fn estimate_extra_memory_usage(
+ &self,
+ geo_stats: &GeoStatistics,
+ spatial_predicate: &SpatialPredicate,
+ options: &SpatialJoinOptions,
+ ) -> usize {
+ // Match the refiner the join actually uses (pinned to `tg`) so the memory
+ // estimate reflects it.
+ self.default.estimate_extra_memory_usage(
+ geo_stats,
+ spatial_predicate,
+ &pin_refiner_to_tg(options.clone()),
+ )
+ }
+
+ fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> {
+ Arc::new(RasterGeometryArrayFactory {
+ target_crs: self.target_crs.clone(),
+ engine: self.engine.clone(),
+ })
+ }
+}
+
+/// Evaluates the operands of a raster/geometry spatial predicate.
+///
+/// The same factory sees both operands of the join. A raster operand is turned
+/// into its footprint — the polygon through the raster's four corners
+/// (see [`RasterGeometryArrayFactory::evaluate_raster`]) — and a geometry operand
+/// is evaluated with the default planar behavior, since the geometry is already
+/// in the target CRS.
+///
+/// A cross-CRS raster footprint densifies each of its four edges
+/// ([`FOOTPRINT_POINTS_PER_EDGE`] interior points) in the raster's own CRS, where
+/// the edges are exact straight lines, then reprojects every densified point into
+/// the target CRS. The indexed and refined footprint therefore follows the
+/// curved image of each edge rather than chording straight across it. Same-CRS
+/// joins keep the exact four-corner footprint (no reprojection, no densification).
+#[derive(Debug)]
+struct RasterGeometryArrayFactory {
+ /// The geometry operand's CRS: the common CRS this join compares in. Raster
+ /// footprints are reprojected into it. `None` when the geometry operand
+ /// carries no CRS (the raster operand must then also be CRS-less).
+ target_crs: Crs,
+ /// The session's CRS engine, used to reproject raster footprints into
+ /// `target_crs`.
+ engine: Arc<dyn CrsEngine + Send + Sync>,
+}
+
+impl EvaluatedGeometryArrayFactory for RasterGeometryArrayFactory {
+ fn try_new_evaluated_array(
+ &self,
+ geometry_array: ArrayRef,
+ sedona_type: &SedonaType,
+ distance_columnar_value: Option<&ColumnarValue>,
+ ) -> Result<EvaluatedGeometryArray> {
+ // Relation predicates (Intersects/Contains/Within) carry no distance;
+ // this factory is only wired up for those by the raster planner.
+ if distance_columnar_value.is_some() {
+ return sedona_internal_err!(
+ "Raster spatial joins do not support a distance predicate"
+ );
+ }
+
+ match sedona_type {
+ SedonaType::Raster => self.evaluate_raster(geometry_array),
+ // The geometry operand is already in the target CRS, so its default
+ // planar evaluation (Cartesian WKB bounds + WKB) is exactly right.
+ _ => EvaluatedGeometryArray::try_new(geometry_array, sedona_type),
+ }
+ }
+}
+
+impl RasterGeometryArrayFactory {
+ /// Evaluate a raster struct array into footprint polygons plus bounding
+ /// rectangles, reprojecting each footprint into the target CRS.
+ fn evaluate_raster(&self, raster_array: ArrayRef) -> Result<EvaluatedGeometryArray> {
+ let struct_array = raster_array
+ .as_any()
+ .downcast_ref::<StructArray>()
+ .ok_or_else(|| {
+ sedona_internal_datafusion_err!("Expected StructArray for raster operand")
+ })?;
+ let rasters = RasterStructArray::try_new(struct_array)
+ .map_err(|e| exec_datafusion_err!("Failed to read raster array: {e}"))?;
+
+ let num_rows = rasters.len();
+ let mut builder = BinaryBuilder::with_capacity(num_rows, num_rows * 96);
+ let mut rects = Vec::with_capacity(num_rows);
+
+ // Reused across rows so densifying each cross-CRS footprint (below)
+ // allocates its point ring once, not per raster.
+ let mut ring = Vec::new();
+
+ let engine: &dyn CrsEngine = self.engine.as_ref();
+ for i in 0..num_rows {
+ // A null raster produces a null footprint that never matches.
+ if rasters.is_null(i) {
+ builder.append_null();
+ rects.push(Bounds2D::empty());
+ continue;
+ }
+
+ let raster = rasters
+ .get(i)
+ .map_err(|e| exec_datafusion_err!("Failed to read raster row {i}: {e}"))?;
+ let rect = self.append_footprint(&raster, engine, &mut ring, &mut builder)?;
+ rects.push(rect);
+ }
+
+ let footprint_array: ArrayRef = Arc::new(builder.finish());
+ EvaluatedGeometryArray::try_new_with_rects(footprint_array, rects, &WKB_GEOMETRY)
+ }
+
+ /// Append one raster's footprint WKB to `builder` and return its bounding
+ /// rectangle, reconciling the raster's CRS against the target CRS.
+ ///
+ /// CRS rules mirror the `RS_*` kernel: equal (or both-absent) CRS compares
+ /// directly; a genuine CRS difference densifies and reprojects the footprint
+ /// into the target CRS; a one-sided CRS is an error. Same-CRS footprints are
+ /// the exact four-corner polygon; cross-CRS footprints densify each edge
+ /// before reprojection so they follow its curve (see the type-level note on
+ /// [`RasterGeometryArrayFactory`]). `ring` is scratch space reused across rows
+ /// for the densified reprojected ring.
+ fn append_footprint(
+ &self,
+ raster: &dyn RasterRef,
+ engine: &dyn CrsEngine,
+ ring: &mut Vec<(f64, f64)>,
+ builder: &mut BinaryBuilder,
+ ) -> Result<Bounds2D> {
+ let raster_crs = resolve_crs(raster.crs())?;
+ let corners = raster_footprint_corners(raster);
+
+ match (raster_crs.as_deref(), self.target_crs.as_deref()) {
+ // No CRS on either side, or identical CRS: compare directly. The
+ // footprint is byte-identical to the one the `RS_*` kernel builds.
+ (None, None) => {
+ write_footprint_wkb(corners, builder)?;
+ builder.append_value([]);
+ Ok(bounds_from_coords(&corners))
+ }
+ (Some(raster_crs), Some(target_crs)) if raster_crs.crs_equals(target_crs) => {
+ write_footprint_wkb(corners, builder)?;
+ builder.append_value([]);
+ Ok(bounds_from_coords(&corners))
+ }
+ // Genuine CRS difference: densify each edge in the raster's own CRS
+ // and reproject every point into the target CRS.
+ (Some(raster_crs), Some(target_crs)) => {
+ append_reprojected_footprint(corners, raster_crs, target_crs, engine, ring, builder)
+ }
+ (Some(_), None) => {
+ exec_err!(
+ "Cannot evaluate spatial predicate: raster has a CRS but the geometry does not"
+ )
+ }
+ (None, Some(_)) => {
+ exec_err!(
+ "Cannot evaluate spatial predicate: geometry has a CRS but the raster does not"
+ )
+ }
+ }
+ }
+}
+
+/// Densify a raster's four `corners` in `from_crs`, reproject every densified
+/// point to `to_crs`, append the resulting footprint polygon to `builder`, and
+/// return its bounding rectangle.
+///
+/// The footprint is densified while it is still in the raster's own CRS, where
+/// each edge is an exact straight line, and every densified point is then
+/// reprojected — so the reprojected ring follows the curved image of each edge
+/// instead of chording straight across it. `ring` is caller-managed scratch
+/// space for the densified ring.
+///
+/// A transform that cannot be built, a transform that fails on any densified
+/// point, or a reprojected point that is non-finite yields a *null* footprint for
+/// this row (`append_null()` + [`Bounds2D::empty`], mirroring the null-raster
+/// path) so the row contributes no matches — rather than writing a garbage
+/// footprint/MBR or aborting the whole join. Recovering a failed reprojection via
+/// a WGS84 fallback (as the `RS_*` kernel does for the comparison as a whole) is a
+/// separate open question and is not attempted here.
+fn append_reprojected_footprint(
+ corners: [(f64, f64); 4],
+ from_crs: &(dyn CoordinateReferenceSystem + Send + Sync),
+ to_crs: &(dyn CoordinateReferenceSystem + Send + Sync),
+ engine: &dyn CrsEngine,
+ ring: &mut Vec<(f64, f64)>,
+ builder: &mut BinaryBuilder,
+) -> Result<Bounds2D> {
+ let Ok(transform) = engine.get_transform_crs_to_crs(
+ &from_crs.to_crs_string(),
+ &to_crs.to_crs_string(),
+ None,
+ "",
+ ) else {
+ builder.append_null();
+ return Ok(Bounds2D::empty());
+ };
+
+ // Densify the four-corner ring in the native CRS (exact straight edges),
+ // then reproject every point so the footprint follows the reprojected curve.
+ densify_footprint_ring(corners, ring);
+ for point in ring.iter_mut() {
+ if transform.transform_coord(point).is_err() {
+ builder.append_null();
+ return Ok(Bounds2D::empty());
+ }
+ }
+
+ // A single non-finite reprojected point (e.g. a point outside the target
+ // projection's valid domain) would poison the footprint WKB and its MBR, so
+ // emit a null footprint for this row instead of writing garbage.
+ if ring.iter().any(|(x, y)| !x.is_finite() || !y.is_finite()) {
+ builder.append_null();
+ return Ok(Bounds2D::empty());
+ }
+
+ write_footprint_ring_wkb(ring, builder)?;
+ builder.append_value([]);
+ Ok(bounds_from_coords(ring))
+}
+
+/// Bounding rectangle of a set of coordinates, accumulated as f64 [`Interval`]s
+/// (the same primitive the geometry operand's bounds use). [`Bounds2D::new`]
+/// enlarges the f32 bounds outward so the rectangle conservatively contains
+/// every input coordinate.
+fn bounds_from_coords(coords: &[(f64, f64)]) -> Bounds2D {
+ let mut x = Interval::empty();
+ let mut y = Interval::empty();
+ for &(cx, cy) in coords {
+ x.update_value(cx);
+ y.update_value(cy);
+ }
+ Bounds2D::new(x, y)
+}
+
+#[cfg(test)]
+mod test {
+ use super::*;
+ use sedona_proj::transform::LazyProjEngine;
+
+ /// The session's default CRS engine. The same-CRS / no-CRS / one-sided-CRS
+ /// paths exercised by these tests never invoke it (reprojection is tested
+ /// separately with a mock engine below), so any engine works here.
+ fn test_engine() -> Arc<dyn CrsEngine + Send + Sync> {
+ Arc::new(LazyProjEngine)
+ }
+
+ #[test]
+ fn bounds_from_coords_covers_corners() {
+ let corners = [(2.0, 3.08), (2.29, 2.4), (2.09, 3.0), (2.2, 2.48)];
+ let bounds = bounds_from_coords(&corners);
+ let ((min_x, max_x), (min_y, max_y)) = bounds.into_inner();
+
+ // f32 bounds must conservatively contain the f64 extent.
+ assert!((min_x as f64) <= 2.0);
+ assert!((max_x as f64) >= 2.29);
+ assert!((min_y as f64) <= 2.4);
+ assert!((max_y as f64) >= 3.08);
+ }
+
+ // --- Evaluator tests over real rasters -------------------------------
+
+ use std::rc::Rc;
+
+ use arrow_array::BinaryArray;
+ use sedona_geometry::bounding_box::BoundingBox;
+ use sedona_geometry::error::SedonaGeometryError;
+ use sedona_geometry::transform::CrsTransform;
+ use sedona_raster::builder::RasterBuilder;
+ use sedona_raster::traits::{BandMetadata, RasterMetadata};
+ use sedona_raster_functions::footprint::write_convexhull_wkb;
+ use sedona_schema::crs::lnglat;
+ use sedona_schema::datatypes::RASTER;
+ use sedona_schema::raster::{BandDataType, StorageType};
+ use sedona_testing::rasters::generate_test_rasters;
+
+ /// A 1x1 raster over world coords (0,0)-(1,1), with `crs` (or none).
+ fn build_unit_raster(crs: Option<&str>) -> arrow_array::StructArray {
+ let mut builder = RasterBuilder::new(1);
+ let metadata = RasterMetadata {
+ width: 1,
+ height: 1,
+ upperleft_x: 0.0,
+ upperleft_y: 1.0,
+ scale_x: 1.0,
+ scale_y: -1.0,
+ skew_x: 0.0,
+ skew_y: 0.0,
+ };
+ builder.start_raster(&metadata, crs).unwrap();
+ builder
+ .start_band(BandMetadata {
+ datatype: BandDataType::UInt8,
+ nodata_value: None,
+ storage_type: StorageType::InDb,
+ outdb_url: None,
+ outdb_band_id: None,
+ })
+ .unwrap();
+ builder.band_data_writer().append_value([0u8]);
+ builder.finish_band().unwrap();
+ builder.finish_raster().unwrap();
+ builder.finish().unwrap()
+ }
+
+ /// Same-CRS: the footprint is byte-identical to the `RS_*` kernel's convex
+ /// hull, its MBR covers the hand-computed corners of raster 1, and a null
+ /// raster yields a null footprint with an empty rect.
+ #[test]
+ fn same_crs_footprint_matches_kernel_and_bounds() {
+ let rasters = generate_test_rasters(3, Some(0)).unwrap();
+
+ // Expected footprint bytes for raster 1, straight from the shared kernel helper.
+ let mut expected_wkb = Vec::new();
+ {
+ let arr = RasterStructArray::try_new(&rasters).unwrap();
+ let raster1 = arr.get(1).unwrap();
+ write_convexhull_wkb(&raster1, &mut expected_wkb).unwrap();
+ }
+
+ let factory = RasterGeometryArrayFactory {
+ target_crs: lnglat(),
+ engine: test_engine(),
+ };
+ let evaluated = factory.evaluate_raster(Arc::new(rasters)).unwrap();
+ let footprints = evaluated
+ .geometry_array()
+ .as_any()
+ .downcast_ref::<BinaryArray>()
+ .unwrap();
+
+ // Null raster -> null footprint, empty rect.
+ assert!(footprints.is_null(0));
+ assert!(evaluated.rect(0).is_empty());
+
+ // Same-CRS footprint is byte-identical to the kernel's convex hull.
+ assert_eq!(footprints.value(1), expected_wkb.as_slice());
+
+ // Raster 1's footprint corners (from GDAL): (2.0, 3.0), (2.2, 3.08),
+ // (2.29, 2.48), (2.09, 2.4). The MBR must conservatively cover them.
+ let ((min_x, max_x), (min_y, max_y)) = evaluated.rect(1).clone().into_inner();
+ assert!((min_x as f64) <= 2.0);
+ assert!((max_x as f64) >= 2.29);
+ assert!((min_y as f64) <= 2.4);
+ assert!((max_y as f64) >= 3.08);
+ }
+
+ /// A raster with a CRS joined against a CRS-less geometry (target CRS
+ /// `None`) is a one-sided CRS: it must error rather than silently compare.
+ #[test]
+ fn raster_with_crs_but_crsless_target_errors() {
+ let rasters = build_unit_raster(Some("OGC:CRS84"));
+ let factory = RasterGeometryArrayFactory {
+ target_crs: None,
+ engine: test_engine(),
+ };
+ let err = factory.evaluate_raster(Arc::new(rasters)).err().unwrap();
+ assert!(err.message().contains("has a CRS but"), "unexpected: {err}");
+ }
+
+ /// A CRS-less raster joined against a geometry that has a CRS is also a
+ /// one-sided CRS and must error.
+ #[test]
+ fn crsless_raster_with_crs_target_errors() {
+ let rasters = build_unit_raster(None);
+ let factory = RasterGeometryArrayFactory {
+ target_crs: lnglat(),
+ engine: test_engine(),
+ };
+ let err = factory.evaluate_raster(Arc::new(rasters)).err().unwrap();
+ assert!(err.message().contains("has a CRS but"), "unexpected: {err}");
+ }
+
+ /// A CRS-less raster with a CRS-less geometry compares directly (no error,
+ /// no reprojection): the footprint is the native convex hull.
+ #[test]
+ fn crsless_both_sides_compares_directly() {
+ let rasters = build_unit_raster(None);
+ let mut expected_wkb = Vec::new();
+ {
+ let arr = RasterStructArray::try_new(&rasters).unwrap();
+ let raster0 = arr.get(0).unwrap();
+ write_convexhull_wkb(&raster0, &mut expected_wkb).unwrap();
+ }
+
+ let factory = RasterGeometryArrayFactory {
+ target_crs: None,
+ engine: test_engine(),
+ };
+ let evaluated = factory
+ .try_new_evaluated_array(Arc::new(rasters), &RASTER, None)
+ .unwrap();
+ let footprints = evaluated
+ .geometry_array()
+ .as_any()
+ .downcast_ref::<BinaryArray>()
+ .unwrap();
+ assert_eq!(footprints.value(0), expected_wkb.as_slice());
+ }
+
+ /// How a [`MockTransform`] treats each corner it is asked to reproject.
+ #[derive(Debug, Clone, Copy)]
+ enum CornerFate {
+ /// Leave the coordinate unchanged (a well-behaved transform).
+ Identity,
+ /// Reproject the corner to a non-finite coordinate.
+ NonFinite,
+ /// Fail while transforming the corner.
+ Error,
+ }
+
+ #[derive(Debug)]
+ struct MockTransform(CornerFate);
+
+ impl CrsTransform for MockTransform {
+ fn transform_coord(&self, coord: &mut (f64, f64)) -> Result<(), SedonaGeometryError> {
+ match self.0 {
+ CornerFate::Identity => Ok(()),
+ CornerFate::NonFinite => {
+ coord.1 = f64::INFINITY;
+ Ok(())
+ }
+ CornerFate::Error => Err(SedonaGeometryError::Invalid("boom".to_string())),
+ }
+ }
+ }
+
+ /// A [`CrsEngine`] that hands back a [`MockTransform`], or fails to build a
+ /// transform at all when `fate` is `None`.
+ #[derive(Debug)]
+ struct MockEngine {
+ fate: Option<CornerFate>,
+ }
+
+ impl CrsEngine for MockEngine {
+ fn get_transform_crs_to_crs(
+ &self,
+ _from: &str,
+ _to: &str,
+ _area_of_interest: Option<BoundingBox>,
+ _options: &str,
+ ) -> Result<Rc<dyn CrsTransform>, SedonaGeometryError> {
+ match self.fate {
+ Some(fate) => Ok(Rc::new(MockTransform(fate))),
+ None => Err(SedonaGeometryError::Invalid("no transform".to_string())),
+ }
+ }
+
+ fn get_transform_pipeline(
+ &self,
+ _pipeline: &str,
+ _options: &str,
+ ) -> Result<Rc<dyn CrsTransform>, SedonaGeometryError> {
+ Err(SedonaGeometryError::Unknown)
+ }
+
+ fn to_projjson(&self, _crs_string: &str) -> Result<String, SedonaGeometryError> {
+ Err(SedonaGeometryError::Unknown)
+ }
+ }
+
+ /// A footprint reprojection that fails — the transform can't be built, a
+ /// densified point's transform errors, or a reprojected point is non-finite —
+ /// yields a null footprint with an empty MBR for that row (mirroring the
+ /// null-raster path) rather than writing a garbage footprint/MBR or aborting
+ /// the whole join. A transform that succeeds with finite points still produces
+ /// a real, non-null footprint.
+ #[test]
+ fn failed_reprojection_yields_null_footprint() {
+ let crs = lnglat();
+ let crs = crs.as_deref().unwrap();
+ let corners = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)];
+
+ for fate in [None, Some(CornerFate::Error), Some(CornerFate::NonFinite)] {
+ let engine = MockEngine { fate };
+ let mut ring = Vec::new();
+ let mut builder = BinaryBuilder::new();
+ let bounds =
+ append_reprojected_footprint(corners, crs, crs, &engine, &mut ring, &mut builder)
+ .unwrap();
+ let footprints = builder.finish();
+ assert!(
+ footprints.is_null(0),
+ "a failed reprojection ({fate:?}) must produce a null footprint"
+ );
+ assert!(
+ bounds.is_empty(),
+ "a null footprint ({fate:?}) must have an empty MBR so it never matches"
+ );
+ }
+
+ // A finite transform still produces a real footprint with a non-empty MBR.
+ let engine = MockEngine {
+ fate: Some(CornerFate::Identity),
+ };
+ let mut ring = Vec::new();
+ let mut builder = BinaryBuilder::new();
+ let bounds =
+ append_reprojected_footprint(corners, crs, crs, &engine, &mut ring, &mut builder)
+ .unwrap();
+ let footprints = builder.finish();
+ assert!(
+ !footprints.is_null(0),
+ "a finite reprojection must produce a real footprint"
+ );
+ assert!(!bounds.is_empty());
+ }
+}
diff --git a/rust/sedona-spatial-join-raster/src/lib.rs b/rust/sedona-spatial-join-raster/src/lib.rs
new file mode 100644
index 000000000..f0a1b84b7
--- /dev/null
+++ b/rust/sedona-spatial-join-raster/src/lib.rs
@@ -0,0 +1,36 @@
+// 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.
+
+//! Optimized spatial joins with a raster operand.
+//!
+//! A raster/geometry spatial predicate (`RS_Intersects`, `RS_Contains`,
+//! `RS_Within`) is accelerated by evaluating each raster into its footprint (the
+//! polygon through the raster's four corners) as a WKB polygon plus a bounding
+//! rectangle. The footprint is reprojected into the geometry operand's CRS so the
+//! R-tree filter and the WKB refiner — both unchanged from the default planar
+//! spatial join — compare footprints and geometries in a common CRS.
+//!
+//! A cross-CRS raster footprint densifies each of its four edges (~10 interior
+//! points) in the raster's own CRS, where the edges are exact straight lines,
+//! then reprojects every densified point into the target CRS. The indexed and
+//! refined footprint therefore follows the curved image of each edge rather than
+//! chording straight across it. Same-CRS joins keep the exact four-corner
+//! footprint (no reprojection, no densification). Antimeridian crossings and
+//! geodesic edges are not modeled.
+
+mod join_provider;
+pub mod physical_planner;
diff --git a/rust/sedona-spatial-join-raster/src/physical_planner.rs b/rust/sedona-spatial-join-raster/src/physical_planner.rs
new file mode 100644
index 000000000..abd4ba025
--- /dev/null
+++ b/rust/sedona-spatial-join-raster/src/physical_planner.rs
@@ -0,0 +1,314 @@
+// 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.
+
+use std::sync::Arc;
+
+use arrow_schema::Schema;
+use datafusion_common::Result;
+use datafusion_physical_expr::PhysicalExpr;
+use datafusion_physical_plan::ExecutionPlan;
+use sedona_common::SedonaOptions;
+use sedona_geometry::transform::CrsEngine;
+use sedona_proj::transform::LazyProjEngine;
+use sedona_query_planner::{
+ spatial_join_physical_planner::{PlanSpatialJoinArgs, SpatialJoinPhysicalPlanner},
+ spatial_predicate::{RelationPredicate, SpatialPredicate, SpatialRelationType},
+};
+use sedona_schema::{crs::Crs, datatypes::SedonaType, matchers::ArgMatcher};
+use sedona_spatial_join::{
+ physical_planner::{repartition_probe_side, should_swap_join_order},
+ SpatialJoinExec,
+};
+
+use crate::join_provider::RasterJoinProvider;
+
+/// [`SpatialJoinPhysicalPlanner`] implementation for raster/geometry spatial joins.
+///
+/// Handles `RS_Intersects`/`RS_Contains`/`RS_Within` predicates where exactly one
+/// operand is a raster and the other is a planar geometry, producing an optimized
+/// [`SpatialJoinExec`] backed by a [`RasterJoinProvider`]. All other predicates
+/// (including raster/raster, for which there is no fixed common CRS to compare in)
+/// are declined with `None` so they fall back to the nested-loop join.
+#[derive(Debug)]
+pub struct RasterSpatialJoinPhysicalPlanner;
+
+impl RasterSpatialJoinPhysicalPlanner {
+ /// Create a new raster join planner.
+ pub fn new() -> Self {
+ Self
+ }
+}
+
+impl Default for RasterSpatialJoinPhysicalPlanner {
+ fn default() -> Self {
+ Self::new()
+ }
+}
+
+impl SpatialJoinPhysicalPlanner for RasterSpatialJoinPhysicalPlanner {
+ fn plan_spatial_join(
+ &self,
+ args: &PlanSpatialJoinArgs<'_>,
+ ) -> Result<Option<Arc<dyn ExecutionPlan>>> {
+ let Some(target_crs) = raster_geometry_target_crs(
+ args.spatial_predicate,
+ &args.physical_left.schema(),
+ &args.physical_right.schema(),
+ )?
+ else {
+ return Ok(None);
+ };
+
+ let should_swap = args.join_type.supports_swap()
+ && should_swap_join_order(
+ args.join_options,
+ args.physical_left.as_ref(),
+ args.physical_right.as_ref(),
+ )?;
+
+ // Repartition the probe side when enabled, mirroring the default planner.
+ let (physical_left, physical_right) = if args.join_options.repartition_probe_side {
+ repartition_probe_side(
+ args.physical_left.clone(),
+ args.physical_right.clone(),
+ args.spatial_predicate,
+ should_swap,
+ )?
+ } else {
+ (args.physical_left.clone(), args.physical_right.clone())
+ };
+
+ // Capture the session's CRS engine at plan time so raster footprint
+ // reprojection honors a user-injected engine. Absent a SedonaOptions
+ // extension, fall back to the process-global PROJ engine via
+ // LazyProjEngine (the same default the session registers).
+ let engine: Arc<dyn CrsEngine + Send + Sync> = args
+ .options
+ .extensions
+ .get::<SedonaOptions>()
+ .map(|opts| opts.runtime.crs_engine().clone())
+ .unwrap_or_else(|| Arc::new(LazyProjEngine));
+
+ let exec = SpatialJoinExec::try_new(
+ physical_left,
+ physical_right,
+ args.spatial_predicate.clone(),
+ args.remainder.cloned(),
+ args.join_type,
+ None,
+ args.join_options,
+ )?
+ .with_spatial_join_provider(Arc::new(RasterJoinProvider::new(target_crs, engine)));
+
+ if should_swap {
+ exec.swap_inputs().map(Some)
+ } else {
+ Ok(Some(Arc::new(exec) as Arc<dyn ExecutionPlan>))
+ }
+ }
+}
+
+/// Return the target CRS (the geometry operand's CRS) when `spatial_predicate`
+/// is a raster/geometry relation this planner can accelerate, or `None` when it
+/// should fall through to another planner.
+///
+/// Handled: a `Relation` predicate of `Intersects`/`Contains`/`Within` with
+/// exactly one raster operand and one planar-geometry operand. The returned
+/// [`Crs`] is the geometry operand's schema CRS, which is `None` when the
+/// geometry carries no CRS.
+fn raster_geometry_target_crs(
+ spatial_predicate: &SpatialPredicate,
+ left_schema: &Schema,
+ right_schema: &Schema,
+) -> Result<Option<Crs>> {
+ let SpatialPredicate::Relation(RelationPredicate {
+ left,
+ right,
+ relation_type,
+ }) = spatial_predicate
+ else {
+ return Ok(None);
+ };
+
+ if !matches!(
+ relation_type,
+ SpatialRelationType::Intersects
+ | SpatialRelationType::Contains
+ | SpatialRelationType::Within
+ ) {
+ return Ok(None);
+ }
+
+ let left_type = operand_sedona_type(left, left_schema)?;
+ let right_type = operand_sedona_type(right, right_schema)?;
+
+ let is_geometry = ArgMatcher::is_geometry();
+ let target = match (
+ matches!(left_type, SedonaType::Raster),
+ matches!(right_type, SedonaType::Raster),
+ ) {
+ // Exactly one raster operand; the other must be a planar geometry.
+ (true, false) if is_geometry.match_type(&right_type) => Some(right_type.crs().clone()),
+ (false, true) if is_geometry.match_type(&left_type) => Some(left_type.crs().clone()),
+ // raster/raster (no fixed common CRS), geometry/geometry, geography, or
+ // item-crs operands are left to other planners / the nested-loop fallback.
+ _ => None,
+ };
+
+ Ok(target)
+}
+
+/// Resolve the [`SedonaType`] an operand expression evaluates to against `schema`.
+fn operand_sedona_type(expr: &Arc<dyn PhysicalExpr>, schema: &Schema) -> Result<SedonaType> {
+ let return_field = expr.return_field(schema)?;
+ SedonaType::from_storage_field(&return_field)
+}
+
+#[cfg(test)]
+mod test {
+ use std::sync::Arc;
+
+ use arrow_schema::Schema;
+ use datafusion_physical_expr::{expressions::Column, PhysicalExpr};
+ use sedona_geometry::types::Edges;
+ use sedona_query_planner::spatial_predicate::{
+ RelationPredicate, SpatialPredicate, SpatialRelationType,
+ };
+ use sedona_schema::crs::{deserialize_crs, lnglat};
+ use sedona_schema::datatypes::{SedonaType, RASTER, WKB_GEOMETRY};
+
+ use super::raster_geometry_target_crs;
+
+ fn schema_with(field_name: &str, sedona_type: &SedonaType) -> Arc<Schema> {
+ Arc::new(Schema::new(vec![sedona_type
+ .to_storage_field(field_name, true)
+ .unwrap()]))
+ }
+
+ fn col(name: &str) -> Arc<dyn PhysicalExpr> {
+ Arc::new(Column::new(name, 0))
+ }
+
+ fn relation(
+ left: Arc<dyn PhysicalExpr>,
+ right: Arc<dyn PhysicalExpr>,
+ relation_type: SpatialRelationType,
+ ) -> SpatialPredicate {
+ SpatialPredicate::Relation(RelationPredicate::new(left, right, relation_type))
+ }
+
+ #[test]
+ fn raster_geometry_resolves_geometry_crs_as_target() {
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let raster_schema = schema_with("raster", &RASTER);
+ let geom_schema = schema_with("geom", &geom_type);
+
+ // (raster, geometry): target is the geometry's CRS.
+ let pred = relation(col("raster"), col("geom"), SpatialRelationType::Intersects);
+ let target = raster_geometry_target_crs(&pred, &raster_schema, &geom_schema)
+ .unwrap()
+ .expect("raster/geometry should be handled");
+ assert!(target
+ .as_deref()
+ .unwrap()
+ .crs_equals(lnglat().as_deref().unwrap()));
+
+ // (geometry, raster): still resolves to the geometry's CRS.
+ let pred = relation(col("geom"), col("raster"), SpatialRelationType::Contains);
+ let target = raster_geometry_target_crs(&pred, &geom_schema, &raster_schema)
+ .unwrap()
+ .expect("geometry/raster should be handled");
+ assert!(target
+ .as_deref()
+ .unwrap()
+ .crs_equals(lnglat().as_deref().unwrap()));
+ }
+
+ #[test]
+ fn crs_less_geometry_yields_crs_less_target() {
+ let raster_schema = schema_with("raster", &RASTER);
+ let geom_schema = schema_with("geom", &WKB_GEOMETRY);
+
+ let pred = relation(col("raster"), col("geom"), SpatialRelationType::Within);
+ let target = raster_geometry_target_crs(&pred, &raster_schema, &geom_schema)
+ .unwrap()
+ .expect("raster/geometry should be handled");
+ assert!(
+ target.is_none(),
+ "CRS-less geometry means a CRS-less target"
+ );
+ }
+
+ #[test]
+ fn different_geometry_crs_is_carried_through() {
+ let geom_type = SedonaType::Wkb(Edges::Planar, deserialize_crs("EPSG:3857").unwrap());
+ let raster_schema = schema_with("raster", &RASTER);
+ let geom_schema = schema_with("geom", &geom_type);
+
+ let pred = relation(col("geom"), col("raster"), SpatialRelationType::Intersects);
+ let target = raster_geometry_target_crs(&pred, &geom_schema, &raster_schema)
+ .unwrap()
+ .expect("raster/geometry should be handled");
+ assert!(target
+ .as_deref()
+ .unwrap()
+ .crs_equals(deserialize_crs("EPSG:3857").unwrap().as_deref().unwrap()));
+ }
+
+ #[test]
+ fn raster_raster_is_declined() {
+ let raster_schema = schema_with("raster", &RASTER);
+ let pred = relation(
+ col("raster"),
+ col("raster"),
+ SpatialRelationType::Intersects,
+ );
+ assert!(
+ raster_geometry_target_crs(&pred, &raster_schema, &raster_schema)
+ .unwrap()
+ .is_none(),
+ "raster/raster has no fixed common CRS and must be declined"
+ );
+ }
+
+ #[test]
+ fn geometry_geometry_is_declined() {
+ let geom_schema = schema_with("geom", &WKB_GEOMETRY);
+ let pred = relation(col("geom"), col("geom"), SpatialRelationType::Intersects);
+ assert!(
+ raster_geometry_target_crs(&pred, &geom_schema, &geom_schema)
+ .unwrap()
+ .is_none(),
+ "pure geometry predicates are left to the default planner"
+ );
+ }
+
+ #[test]
+ fn unsupported_relation_type_is_declined() {
+ let geom_type = SedonaType::Wkb(Edges::Planar, lnglat());
+ let raster_schema = schema_with("raster", &RASTER);
+ let geom_schema = schema_with("geom", &geom_type);
+
+ // Equals is not an RS_ predicate.
+ let pred = relation(col("raster"), col("geom"), SpatialRelationType::Equals);
+ assert!(
+ raster_geometry_target_crs(&pred, &raster_schema, &geom_schema)
+ .unwrap()
+ .is_none()
+ );
+ }
+}
diff --git a/rust/sedona-spatial-join-raster/tests/spatial_join_integration.rs b/rust/sedona-spatial-join-raster/tests/spatial_join_integration.rs
new file mode 100644
index 000000000..9c72be667
--- /dev/null
+++ b/rust/sedona-spatial-join-raster/tests/spatial_join_integration.rs
@@ -0,0 +1,599 @@
+// 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.
+
+//! Integration tests for raster/geometry spatial joins.
+//!
+//! The strict-parity tests run the same query two ways over the same data:
+//!
+//! * **Optimized** — the raster extension planner turns the `RS_*` join into a
+//! `SpatialJoinExec` that evaluates each raster into its (reprojected)
+//! footprint polygon.
+//! * **Nested-loop** — with the spatial-join optimizer/planner absent, the join
+//! stays a `NestedLoopJoinExec` that evaluates the `RS_*` UDF kernel directly.
+//!
+//! When the rasters and the geometry share a CRS (as they do here — both are
+//! lng/lat) neither side reprojects, and both refine with the same `tg` engine,
+//! so the two row sets must be byte-identical. This asserts the accelerated path
+//! does not silently diverge from the established kernel.
+//!
+//! `proj-sys` is enabled as a dev-dependency, so the `RS_*` kernel's global PROJ
+//! engine initializes for the nested-loop side.
+
+use std::sync::Arc;
+
+use arrow_array::{Int32Array, RecordBatch, StructArray};
+use arrow_schema::{Field, Schema};
+use datafusion::{
+ catalog::MemTable,
+ execution::SessionStateBuilder,
+ physical_plan::{displayable, ExecutionPlan},
+ prelude::{SessionConfig, SessionContext},
+};
+use datafusion_common::cast::as_int32_array;
+use datafusion_common::tree_node::{TreeNode, TreeNodeRecursion};
+use datafusion_common::Result;
+use rstest::rstest;
+use sedona_common::option::SpatialJoinOptions;
+use sedona_common::SedonaOptions;
+use sedona_geometry::transform::CrsEngine;
+use sedona_geometry::types::Edges;
+use sedona_proj::transform::{with_global_proj_engine, LazyProjEngine};
+use sedona_query_planner::{
+ optimizer::register_spatial_join_logical_optimizer, query_planner::SedonaQueryPlanner,
+};
+use sedona_raster::affine_transformation::to_world_coordinate;
+use sedona_raster::array::RasterStructArray;
+use sedona_raster::traits::RasterRef;
+use sedona_raster_functions::footprint::{densify_footprint_ring, FOOTPRINT_POINTS_PER_EDGE};
+use sedona_schema::crs::lnglat;
+use sedona_schema::datatypes::{SedonaType, RASTER};
+use sedona_schema::raster::BandDataType;
+use sedona_spatial_join::SpatialJoinExec;
+use sedona_spatial_join_raster::physical_planner::RasterSpatialJoinPhysicalPlanner;
+use sedona_testing::create::create_array;
+use sedona_testing::raster_spec::RasterSpec;
+use sedona_testing::rasters::generate_test_rasters;
+
+/// Geometry type sharing the CRS that `generate_test_rasters` stamps on its
+/// rasters (lng/lat), so the predicate compares in a common CRS.
+fn geom_type() -> SedonaType {
+ SedonaType::Wkb(Edges::Planar, lnglat())
+}
+
+/// Register a raster table `r(rid INT, raster RASTER)`.
+///
+/// `generate_test_rasters(3, Some(0))` yields: rid 0 a null raster, rid 1 a
+/// footprint over x∈[2.0, 2.29] y∈[2.4, 3.08], rid 2 over x∈[3.0, 3.84]
+/// y∈[2.4, 4.24]. All non-null rasters carry the lng/lat CRS.
+fn register_raster_table(ctx: &SessionContext) -> Result<()> {
+ let rasters = generate_test_rasters(3, Some(0)).unwrap();
+ let rid = Int32Array::from(vec![0, 1, 2]);
+
+ let schema = Arc::new(Schema::new(vec![
+ Field::new("rid", arrow_schema::DataType::Int32, false),
+ RASTER.to_storage_field("raster", true)?,
+ ]));
+ let batch = RecordBatch::try_new(schema.clone(), vec![Arc::new(rid), Arc::new(rasters)])?;
+ let table = MemTable::try_new(schema, vec![vec![batch]])?;
+ ctx.register_table("r", Arc::new(table))?;
+ Ok(())
+}
+
+/// Register a geometry table `g(gid INT, geom GEOMETRY)` covering centers,
+/// corners/edges, an enclosing polygon, a non-intersecting zone, and a null.
+fn register_geom_table(ctx: &SessionContext) -> Result<()> {
+ let geom_type = geom_type();
+ let gid = Int32Array::from(vec![10, 20, 30, 40, 50, 60, 70]);
+ let geom = create_array(
+ &[
+ Some("POINT (2.15 2.75)"), // inside raster 1
+ Some("POINT (3.4 3.2)"), // inside raster 2
+ Some("POINT (2.0 3.0)"), // on a raster-1 corner (boundary/edge)
+ Some("POINT (0 0)"), // outside all rasters
+ Some("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"), // encloses rasters 1 and 2
+ Some("POLYGON ((2.0 2.4, 2.29 2.4, 2.29 3.08, 2.0 3.08, 2.0 2.4))"), // around raster 1
+ None, // null geometry never matches
+ ],
+ &geom_type,
+ );
+
+ let schema = Arc::new(Schema::new(vec![
+ Field::new("gid", arrow_schema::DataType::Int32, false),
+ geom_type.to_storage_field("geom", true)?,
+ ]));
+ let batch = RecordBatch::try_new(schema.clone(), vec![Arc::new(gid), geom])?;
+ let table = MemTable::try_new(schema, vec![vec![batch]])?;
+ ctx.register_table("g", Arc::new(table))?;
+ Ok(())
+}
+
+/// Build a context with the raster UDFs registered but no tables. When
+/// `optimized`, the spatial-join logical optimizer and the raster physical
+/// planner are registered (and spatial join enabled), producing a
+/// `SpatialJoinExec`. Otherwise the join stays a `NestedLoopJoinExec`.
+fn build_context(optimized: bool) -> Result<SessionContext> {
+ let mut session_config = SessionConfig::from_env()?.with_batch_size(16);
+ session_config = session_config.with_option_extension(SedonaOptions::default());
+
+ // Install a real CRS engine, mirroring the engine a normal session
+ // registers. Both the accelerated join and the RS_* kernel read the engine
+ // from the session options to reproject cross-CRS footprints/geometries.
+ {
+ let opts = session_config
+ .options_mut()
+ .extensions
+ .get_mut::<SedonaOptions>()
+ .unwrap();
+ opts.runtime = opts.runtime.with_crs_engine(Arc::new(LazyProjEngine));
+ }
+
+ let mut state_builder = SessionStateBuilder::new();
+ if optimized {
+ state_builder = register_spatial_join_logical_optimizer(state_builder)?;
+ state_builder = state_builder.with_query_planner(Arc::new(
+ SedonaQueryPlanner::new().with_spatial_join_physical_planner(Arc::new(
+ RasterSpatialJoinPhysicalPlanner::new(),
+ )),
+ ));
+ let opts = session_config
+ .options_mut()
+ .extensions
+ .get_mut::<SedonaOptions>()
+ .unwrap();
+ opts.spatial_join = SpatialJoinOptions::default();
+ }
+
+ let state = state_builder.with_config(session_config).build();
+ let ctx = SessionContext::new_with_state(state);
+
+ // Register the raster functions (RS_Intersects/RS_Contains/RS_Within).
+ let function_set = sedona_raster_functions::register::default_function_set();
+ function_set.scalar_udfs().for_each(|udf| {
+ ctx.register_udf(udf.clone().into());
+ });
+
+ Ok(ctx)
+}
+
+/// Build a context with the same-CRS `r`/`g` tables registered.
+fn setup_context(optimized: bool) -> Result<SessionContext> {
+ let ctx = build_context(optimized)?;
+ register_raster_table(&ctx)?;
+ register_geom_table(&ctx)?;
+ Ok(ctx)
+}
+
+fn count_spatial_join_execs(plan: &Arc<dyn ExecutionPlan>) -> Result<usize> {
+ let mut count = 0;
+ plan.apply(|node| {
+ if node.as_any().downcast_ref::<SpatialJoinExec>().is_some() {
+ count += 1;
+ }
+ Ok(TreeNodeRecursion::Continue)
+ })?;
+ Ok(count)
+}
+
+async fn run_query(optimized: bool, sql: &str) -> Result<RecordBatch> {
+ let ctx = setup_context(optimized)?;
+ let df = ctx.sql(sql).await?;
+ let schema = df.schema().as_arrow().clone();
+ let plan = df.clone().create_physical_plan().await?;
+
+ let num_spatial_joins = count_spatial_join_execs(&plan)?;
+ let physical_str = displayable(plan.as_ref()).indent(true).to_string();
+ if optimized {
+ assert_eq!(
+ num_spatial_joins, 1,
+ "expected an optimized SpatialJoinExec, got:\n{physical_str}"
+ );
+ } else {
+ assert_eq!(
+ num_spatial_joins, 0,
+ "nested-loop baseline should not contain SpatialJoinExec, got:\n{physical_str}"
+ );
+ assert!(
+ physical_str.contains("NestedLoopJoinExec"),
+ "nested-loop baseline should use NestedLoopJoinExec, got:\n{physical_str}"
+ );
+ }
+
+ let batches = df.collect().await?;
+ Ok(datafusion::arrow::compute::concat_batches(
+ &Arc::new(schema),
+ &batches,
+ )?)
+}
+
+/// Assert the optimized `SpatialJoinExec` returns byte-identical rows to the
+/// nested-loop kernel for a same-CRS raster/geometry join. Covers all three
+/// `RS_*` predicates, both operand orderings, and inner/left/right joins.
+#[rstest]
+#[tokio::test]
+async fn same_crs_strict_parity(
+ #[values("RS_Intersects", "RS_Contains", "RS_Within")] func: &str,
+ #[values("r.raster, g.geom", "g.geom, r.raster")] operands: &str,
+ #[values("INNER JOIN", "LEFT OUTER JOIN", "RIGHT OUTER JOIN")] join_type: &str,
+) -> Result<()> {
+ let sql = format!(
+ "SELECT r.rid, g.gid FROM r {join_type} g ON {func}({operands}) ORDER BY r.rid, g.gid"
+ );
+
+ let optimized = run_query(true, &sql).await?;
+ let nested_loop = run_query(false, &sql).await?;
+
+ assert_eq!(
+ optimized, nested_loop,
+ "optimized SpatialJoinExec must match the nested-loop kernel exactly for {sql}"
+ );
+ Ok(())
+}
+
+/// Anchor the parity comparison to hand-computed rows so a bug shared by both
+/// paths cannot pass. `RS_Intersects(raster, geom)` inner join: raster 1 meets
+/// its center point (gid 10), the enclosing polygon (gid 50), and its bbox
+/// polygon (gid 60), plus the corner point (gid 30); raster 2 meets its center
+/// (gid 20) and the enclosing polygon (gid 50). The null raster (rid 0) and the
+/// null geometry (gid 70) never match.
+#[tokio::test]
+async fn intersects_produces_expected_rows() -> Result<()> {
+ let sql = "SELECT r.rid, g.gid FROM r JOIN g ON RS_Intersects(r.raster, g.geom) \
+ ORDER BY r.rid, g.gid";
+ let result = run_query(true, sql).await?;
+
+ datafusion::assert_batches_eq!(
+ [
+ "+-----+-----+",
+ "| rid | gid |",
+ "+-----+-----+",
+ "| 1 | 10 |",
+ "| 1 | 30 |",
+ "| 1 | 50 |",
+ "| 1 | 60 |",
+ "| 2 | 20 |",
+ "| 2 | 50 |",
+ "+-----+-----+",
+ ],
+ &[result]
+ );
+ Ok(())
+}
+
+/// Cross-CRS: verify the accelerated path against an independent,
+/// construction-based reference. The raster is skewed in EPSG:3857; the geometry
+/// operand is lng/lat.
+///
+/// The accelerated footprint densifies each edge of the raster's four-corner quad
+/// in the native CRS and reprojects every point into lng/lat, so the footprint
+/// follows the curve of each reprojected edge rather than chording across it. The
+/// reference builds that same densified, reprojected ring, then places test
+/// points whose membership is known by construction: the centroid and the
+/// halfway-to-each-corner points are strictly inside the footprint, while the
+/// reflection of the centroid through each corner is strictly outside. These
+/// points sit well away from the edges, so the expected result is exact (no
+/// tolerance) no matter how finely the edges are modeled. The reference
+/// reprojects in the same 3857 -> lng/lat direction as the accelerated path, so a
+/// wrong transform direction places the accelerated footprint elsewhere and none
+/// of the (correctly placed) inside points match — this validates transform
+/// direction and the footprint/refine pipeline against PROJ (the same engine a
+/// rasterio/geopandas reference would use).
+#[tokio::test]
+async fn cross_crs_matches_constructed_reference() -> Result<()> {
+ let raster = build_skewed_raster_3857();
+
+ // Footprint corners in the raster's native CRS (EPSG:3857). These are exact.
+ let native_corners = skewed_raster_native_corners(&raster);
+
+ // Densify the four-corner footprint in the native CRS and reproject every
+ // point into lng/lat: this is exactly the footprint the accelerated path
+ // builds (each edge follows its reprojected curve, not a straight chord).
+ let lnglat_crs = lnglat().unwrap().to_crs_string();
+ let ref_ring: Vec<(f64, f64)> = with_global_proj_engine(|engine| {
+ let transform = engine
+ .get_transform_crs_to_crs("EPSG:3857", &lnglat_crs, None, "")
+ .map_err(|e| datafusion_common::exec_datafusion_err!("{e}"))?;
+ let mut ring = Vec::new();
+ densify_footprint_ring(native_corners, &mut ring);
+ for point in ring.iter_mut() {
+ transform
+ .transform_coord(point)
+ .map_err(|e| datafusion_common::exec_datafusion_err!("{e}"))?;
+ }
+ Ok(ring)
+ })?;
+
+ // The four reprojected corners are the densified ring's vertices at the start
+ // of each edge. Test points are placed relative to these corners and the
+ // centroid — deep inside or far outside — so their membership is exact
+ // regardless of how finely the edges are modeled.
+ let corner_stride = FOOTPRINT_POINTS_PER_EDGE + 1;
+ let hull: [(f64, f64); 4] = [
+ ref_ring[0],
+ ref_ring[corner_stride],
+ ref_ring[2 * corner_stride],
+ ref_ring[3 * corner_stride],
+ ];
+ let cx = hull.iter().map(|c| c.0).sum::<f64>() / 4.0;
+ let cy = hull.iter().map(|c| c.1).sum::<f64>() / 4.0;
+
+ // lng/lat test points with membership known by construction against the
+ // convex chord-quad `hull`.
+ let mut projected: Vec<(i32, (f64, f64), bool)> = vec![(0, (cx, cy), true)];
+ for (i, &(x, y)) in hull.iter().enumerate() {
+ // Halfway from centroid to a vertex: strictly inside a convex polygon.
+ projected.push((
+ 10 + i as i32,
+ (cx + 0.5 * (x - cx), cy + 0.5 * (y - cy)),
+ true,
+ ));
+ // Reflection of the centroid through a vertex: strictly outside.
+ projected.push((
+ 20 + i as i32,
+ (cx + 2.0 * (x - cx), cy + 2.0 * (y - cy)),
+ false,
+ ));
+ }
+
+ // Expected matches: the gids of the inside points.
+ let mut expected_inside: Vec<i32> = projected
+ .iter()
+ .filter(|&&(_, _, inside)| inside)
+ .map(|&(gid, _, _)| gid)
+ .collect();
+ expected_inside.sort_unstable();
+
+ // Build the lng/lat geometry table.
+ let geom_type = geom_type();
+ let gids: Vec<i32> = projected.iter().map(|&(gid, _, _)| gid).collect();
+ let wkts: Vec<String> = projected
+ .iter()
+ .map(|&(_, (lon, lat), _)| format!("POINT ({lon} {lat})"))
+ .collect();
+ let wkt_opts: Vec<Option<&str>> = wkts.iter().map(|w| Some(w.as_str())).collect();
+ let geom = create_array(&wkt_opts, &geom_type);
+ let geom_schema = Arc::new(Schema::new(vec![
+ Field::new("gid", arrow_schema::DataType::Int32, false),
+ geom_type.to_storage_field("geom", true)?,
+ ]));
+ let geom_batch = RecordBatch::try_new(
+ geom_schema.clone(),
+ vec![Arc::new(Int32Array::from(gids)), geom],
+ )?;
+
+ // Register the 3857 raster + lng/lat geometry tables in an optimized context.
+ let ctx = build_context(true)?;
+ let raster_schema = Arc::new(Schema::new(vec![RASTER.to_storage_field("raster", true)?]));
+ let raster_batch = RecordBatch::try_new(raster_schema.clone(), vec![Arc::new(raster)])?;
+ ctx.register_table(
+ "r",
+ Arc::new(MemTable::try_new(raster_schema, vec![vec![raster_batch]])?),
+ )?;
+ ctx.register_table(
+ "g",
+ Arc::new(MemTable::try_new(geom_schema, vec![vec![geom_batch]])?),
+ )?;
+
+ let sql = "SELECT g.gid FROM r JOIN g ON RS_Intersects(r.raster, g.geom) ORDER BY g.gid";
+ let df = ctx.sql(sql).await?;
+ let plan = df.clone().create_physical_plan().await?;
+ assert_eq!(
+ count_spatial_join_execs(&plan)?,
+ 1,
+ "cross-CRS raster join should use the optimized SpatialJoinExec"
+ );
+ let batches = df.collect().await?;
+
+ let mut got: Vec<i32> = Vec::new();
+ for batch in &batches {
+ let col = as_int32_array(batch.column(0))?;
+ for i in 0..col.len() {
+ got.push(col.value(i));
+ }
+ }
+ got.sort_unstable();
+
+ assert_eq!(
+ got, expected_inside,
+ "cross-CRS accelerated join must match the constructed reference"
+ );
+
+ // Sanity: we exercised both inside (5) and outside (4) points.
+ assert_eq!(expected_inside.len(), 5);
+ assert_eq!(got.len(), 5);
+ assert_eq!(projected.len(), 9);
+ Ok(())
+}
+
+/// A 4x4 raster in EPSG:3857 with skew, so its footprint is a non-axis-aligned
+/// quadrilateral whose edges curve when reprojected to lng/lat.
+fn build_skewed_raster_3857() -> StructArray {
+ RasterSpec::d2(4, 4)
+ .crs(Some("EPSG:3857"))
+ // GDAL geotransform [origin_x, scale_x, skew_x, origin_y, skew_y, scale_y].
+ .transform([0.0, 100_000.0, 30_000.0, 2_000_000.0, 20_000.0, -100_000.0])
+ .band(BandDataType::UInt8)
+ .build()
+}
+
+/// The four footprint corners of the single-row `raster` in its native CRS, in
+/// ring order (upper-left, upper-right, lower-right, lower-left).
+fn skewed_raster_native_corners(raster: &StructArray) -> [(f64, f64); 4] {
+ let rasters = RasterStructArray::try_new(raster).unwrap();
+ let r = rasters.get(0).unwrap();
+ let w = r.metadata().width();
+ let h = r.metadata().height();
+ [
+ to_world_coordinate(&r, 0, 0),
+ to_world_coordinate(&r, w, 0),
+ to_world_coordinate(&r, w, h),
+ to_world_coordinate(&r, 0, h),
+ ]
+}
+
+/// True if `p` lies inside (or on the boundary of) the convex quadrilateral
+/// `quad`, tested by requiring the cross product of each directed edge with the
+/// vector to `p` to keep a consistent sign.
+fn point_in_convex_quad(quad: &[(f64, f64); 4], p: (f64, f64)) -> bool {
+ let mut sign = 0.0f64;
+ for i in 0..4 {
+ let a = quad[i];
+ let b = quad[(i + 1) % 4];
+ let cross = (b.0 - a.0) * (p.1 - a.1) - (b.1 - a.1) * (p.0 - a.0);
+ if cross != 0.0 {
+ if sign == 0.0 {
+ sign = cross.signum();
+ } else if cross.signum() != sign {
+ return false;
+ }
+ }
+ }
+ true
+}
+
+/// Densification catches a point in the sliver between a footprint edge's true
+/// reprojected curve and the straight chord through its two reprojected corners.
+///
+/// Reprojecting the skewed EPSG:3857 raster into lng/lat bows at least one edge
+/// outward: its true (curved) boundary bulges past the straight chord between the
+/// two reprojected corners. A point placed just inside that curve is inside the
+/// true footprint but outside the four-corner chord hull, so a join that only
+/// reprojected the corners would miss it. Because the accelerated join densifies
+/// each edge before reprojection, the footprint follows the curve and the point
+/// matches — the accuracy gain densification buys.
+#[tokio::test]
+async fn densification_catches_curved_edge_sliver() -> Result<()> {
+ let raster = build_skewed_raster_3857();
+ let native_corners = skewed_raster_native_corners(&raster);
+
+ // Reproject the four corners (the chord hull) and each edge's native midpoint
+ // (a point on the true reprojected curve) into lng/lat.
+ let lnglat_crs = lnglat().unwrap().to_crs_string();
+ let mut hull = native_corners;
+ let mut curve_mids = [(0.0, 0.0); 4];
+ with_global_proj_engine(|engine| {
+ let transform = engine
+ .get_transform_crs_to_crs("EPSG:3857", &lnglat_crs, None, "")
+ .map_err(|e| datafusion_common::exec_datafusion_err!("{e}"))?;
+ for corner in hull.iter_mut() {
+ transform
+ .transform_coord(corner)
+ .map_err(|e| datafusion_common::exec_datafusion_err!("{e}"))?;
+ }
+ for (e, mid) in curve_mids.iter_mut().enumerate() {
+ let a = native_corners[e];
+ let b = native_corners[(e + 1) % 4];
+ let mut m = ((a.0 + b.0) / 2.0, (a.1 + b.1) / 2.0);
+ transform
+ .transform_coord(&mut m)
+ .map_err(|e| datafusion_common::exec_datafusion_err!("{e}"))?;
+ *mid = m;
+ }
+ Ok(())
+ })?;
+
+ let cx = hull.iter().map(|c| c.0).sum::<f64>() / 4.0;
+ let cy = hull.iter().map(|c| c.1).sum::<f64>() / 4.0;
+
+ // Find the edge whose true curve bows farthest outward (its midpoint lies
+ // outside the chord hull) and build a point just inside that curve.
+ let mut best: Option<(f64, (f64, f64))> = None; // (outward depth, test point)
+ for e in 0..4 {
+ let a = hull[e];
+ let b = hull[(e + 1) % 4];
+ let curve_mid = curve_mids[e];
+ let dx = b.0 - a.0;
+ let dy = b.1 - a.1;
+ let len = (dx * dx + dy * dy).sqrt();
+ // Unit left normal of the chord a -> b.
+ let nx = -dy / len;
+ let ny = dx / len;
+ // Signed perpendicular offset of the curve midpoint and the centroid from
+ // the chord line. Opposite signs => the curve bows away from the centroid,
+ // i.e. outward past the chord hull.
+ let perp_curve = (curve_mid.0 - a.0) * nx + (curve_mid.1 - a.1) * ny;
+ let perp_centroid = (cx - a.0) * nx + (cy - a.1) * ny;
+ if perp_curve == 0.0 || perp_curve.signum() == perp_centroid.signum() {
+ continue; // straight or bows inward: no outward sliver on this edge
+ }
+ let depth = perp_curve.abs();
+ // Step from the true curve halfway back toward the interior: still outside
+ // the chord hull (< depth away), but inside the densified footprint (the
+ // densified edge hugs the true curve to within ~depth/100).
+ let inward = perp_centroid.signum();
+ let test_pt = (
+ curve_mid.0 + 0.5 * depth * nx * inward,
+ curve_mid.1 + 0.5 * depth * ny * inward,
+ );
+ let better = match best {
+ Some((d, _)) => depth > d,
+ None => true,
+ };
+ if better {
+ best = Some((depth, test_pt));
+ }
+ }
+
+ let (_, test_pt) = best.expect("reprojecting the skewed raster must bow an edge outward");
+
+ // The point is outside the four-corner chord hull: a join that reprojected
+ // only the corners would miss it.
+ assert!(
+ !point_in_convex_quad(&hull, test_pt),
+ "test point {test_pt:?} must lie outside the four-corner chord hull"
+ );
+
+ // Register the raster + the single sliver point in an optimized context.
+ let ctx = build_context(true)?;
+ let raster_schema = Arc::new(Schema::new(vec![RASTER.to_storage_field("raster", true)?]));
+ let raster_batch = RecordBatch::try_new(raster_schema.clone(), vec![Arc::new(raster)])?;
+ ctx.register_table(
+ "r",
+ Arc::new(MemTable::try_new(raster_schema, vec![vec![raster_batch]])?),
+ )?;
+
+ let geom_type = geom_type();
+ let wkt = format!("POINT ({} {})", test_pt.0, test_pt.1);
+ let geom = create_array(&[Some(wkt.as_str())], &geom_type);
+ let geom_schema = Arc::new(Schema::new(vec![
+ Field::new("gid", arrow_schema::DataType::Int32, false),
+ geom_type.to_storage_field("geom", true)?,
+ ]));
+ let geom_batch = RecordBatch::try_new(
+ geom_schema.clone(),
+ vec![Arc::new(Int32Array::from(vec![1])), geom],
+ )?;
+ ctx.register_table(
+ "g",
+ Arc::new(MemTable::try_new(geom_schema, vec![vec![geom_batch]])?),
+ )?;
+
+ let sql = "SELECT g.gid FROM r JOIN g ON RS_Intersects(r.raster, g.geom)";
+ let df = ctx.sql(sql).await?;
+ let plan = df.clone().create_physical_plan().await?;
+ assert_eq!(
+ count_spatial_join_execs(&plan)?,
+ 1,
+ "cross-CRS raster join should use the optimized SpatialJoinExec"
+ );
+ let batches = df.collect().await?;
+ let matched: usize = batches.iter().map(|b| b.num_rows()).sum();
+ assert_eq!(
+ matched, 1,
+ "the densified footprint must match the sliver point that the four-corner chord hull misses"
+ );
+ Ok(())
+}
diff --git a/rust/sedona/Cargo.toml b/rust/sedona/Cargo.toml
index a868009de..1cf62bcfa 100644
--- a/rust/sedona/Cargo.toml
+++ b/rust/sedona/Cargo.toml
@@ -44,7 +44,7 @@ proj = ["sedona-proj/proj-sys"]
gdal = ["sedona-gdal/gdal-sys"]
# Enable bindgen for GDAL versions without pre-built bindings (e.g., GDAL 3.13+).
bindgen = ["sedona-gdal/bindgen"]
-spatial-join = ["dep:sedona-spatial-join"]
+spatial-join = ["dep:sedona-spatial-join", "dep:sedona-spatial-join-raster"]
gpu = ["dep:sedona-spatial-join-gpu", "sedona-spatial-join-gpu/gpu"]
s2geography = ["dep:sedona-s2geography", "dep:sedona-spatial-join-geography"]
@@ -91,6 +91,7 @@ sedona-schema = { workspace = true }
sedona-spatial-join = { workspace = true, optional = true }
sedona-spatial-join-gpu = { workspace = true, optional = true }
sedona-spatial-join-geography = { workspace = true, optional = true }
+sedona-spatial-join-raster = { workspace = true, optional = true }
sedona-query-planner = { workspace = true }
sedona-s2geography = { workspace = true, optional = true }
sedona-testing = { workspace = true }
diff --git a/rust/sedona/src/context.rs b/rust/sedona/src/context.rs
index 8970aab4b..ae0b69efc 100644
--- a/rust/sedona/src/context.rs
+++ b/rust/sedona/src/context.rs
@@ -178,6 +178,21 @@ impl SedonaContext {
GpuSpatialJoinPhysicalPlanner::new(),
));
}
+
+ // Register the raster join last so it is consulted before the GPU
+ // planner (planners are tried in reverse registration order). It
+ // recognizes raster/geometry predicates and declines everything else
+ // with `None`, so a raster predicate is handled here rather than
+ // reaching the GPU planner — which would otherwise hard-error on a
+ // raster operand when `GpuOptions.fallback_to_cpu` is false. Anything
+ // it declines falls through to the GPU/geography/default planners.
+ {
+ use sedona_spatial_join_raster::physical_planner::RasterSpatialJoinPhysicalPlanner;
+
+ planner = planner.with_spatial_join_physical_planner(Arc::new(
+ RasterSpatialJoinPhysicalPlanner::new(),
+ ));
+ }
}
// Register the geography bounder for spherical edge types. This enables geography
diff --git a/rust/sedona/tests/raster_spatial_join.rs b/rust/sedona/tests/raster_spatial_join.rs
new file mode 100644
index 000000000..de2f1e948
--- /dev/null
+++ b/rust/sedona/tests/raster_spatial_join.rs
@@ -0,0 +1,202 @@
+// 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.
+
+//! End-to-end test for raster-geometry optimized spatial joins.
+//!
+//! A join predicated on `RS_Intersects(raster, geom)` is recognized as a spatial
+//! predicate (so it becomes a `SpatialJoin` extension node) and is planned as the
+//! optimized `SpatialJoinExec` by the raster extension planner, which evaluates
+//! each raster into its footprint polygon rather than falling back to a
+//! `NestedLoopJoinExec`.
+//!
+//! The row-correctness case is gated behind the `proj` feature because the raster
+//! footprint evaluator (like the `RS_Intersects` kernel) initializes the global
+//! PROJ engine (which needs `sedona-proj/proj-sys`); the plan-shape case only
+//! plans and so runs without PROJ.
+
+// The spatial-join planner/optimizer are only wired in when this feature is enabled.
+#![cfg(feature = "spatial-join")]
+
+use std::sync::Arc;
+
+use arrow_array::{Int32Array, RecordBatch};
+use arrow_schema::{Field, Schema};
+use datafusion::catalog::MemTable;
+use datafusion::physical_plan::displayable;
+use sedona::context::SedonaContext;
+use sedona_geometry::types::Edges;
+use sedona_schema::crs::lnglat;
+use sedona_schema::datatypes::{SedonaType, RASTER};
+use sedona_testing::create::create_array;
+use sedona_testing::rasters::generate_test_rasters;
+
+/// Geometry type sharing the CRS that `generate_test_rasters` stamps on its rasters
+/// (lng/lat), so the predicate compares in a common CRS without reprojection.
+fn geom_type() -> SedonaType {
+ SedonaType::Wkb(Edges::Planar, lnglat())
+}
+
+/// Register a raster table `r(rid INT, raster RASTER)`.
+///
+/// `generate_test_rasters(3, Some(0))` yields: index 0 null, index 1 with a convex
+/// hull around x∈[2.0, 2.29] y∈[2.4, 3.08], index 2 around x∈[3.0, 3.84] y∈[2.4, 4.24].
+fn register_raster_table(ctx: &SedonaContext) {
+ let rasters = generate_test_rasters(3, Some(0)).unwrap();
+ let rid = Int32Array::from(vec![0, 1, 2]);
+
+ let raster_field = RASTER.to_storage_field("raster", true).unwrap();
+ let schema = Arc::new(Schema::new(vec![
+ Field::new("rid", arrow_schema::DataType::Int32, false),
+ raster_field,
+ ]));
+ let batch =
+ RecordBatch::try_new(schema.clone(), vec![Arc::new(rid), Arc::new(rasters)]).unwrap();
+ let table = MemTable::try_new(schema, vec![vec![batch]]).unwrap();
+ ctx.ctx.register_table("r", Arc::new(table)).unwrap();
+}
+
+/// Register a geometry table `g(gid INT, geom GEOMETRY)`.
+fn register_geom_table(ctx: &SedonaContext) {
+ let geom_type = geom_type();
+ let gid = Int32Array::from(vec![10, 20, 30]);
+ let geom = create_array(
+ &[
+ Some("POINT (2.15 2.75)"), // inside raster 1 only
+ Some("POLYGON ((0 0, 10 0, 10 10, 0 10, 0 0))"), // contains rasters 1 and 2
+ Some("POINT (0 0)"), // outside all rasters
+ ],
+ &geom_type,
+ );
+
+ let geom_field = geom_type.to_storage_field("geom", true).unwrap();
+ let schema = Arc::new(Schema::new(vec![
+ Field::new("gid", arrow_schema::DataType::Int32, false),
+ geom_field,
+ ]));
+ let batch = RecordBatch::try_new(schema.clone(), vec![Arc::new(gid), geom]).unwrap();
+ let table = MemTable::try_new(schema, vec![vec![batch]]).unwrap();
+ ctx.ctx.register_table("g", Arc::new(table)).unwrap();
+}
+
+const RASTER_JOIN_SQL: &str = "SELECT r.rid, g.gid \
+ FROM r JOIN g ON RS_Intersects(r.raster, g.geom) \
+ ORDER BY r.rid, g.gid";
+
+/// The `RS_Intersects` join is recognized as a spatial predicate (so the optimized
+/// logical plan contains the `SpatialJoin` extension node) and is planned as the
+/// optimized `SpatialJoinExec` by the raster extension planner, not the
+/// `NestedLoopJoinExec` fallback. Planning only, so no PROJ required.
+#[tokio::test]
+async fn raster_join_uses_optimized_spatial_join() {
+ let ctx = SedonaContext::new_local_interactive().await.unwrap();
+ register_raster_table(&ctx);
+ register_geom_table(&ctx);
+
+ let df = ctx.sql(RASTER_JOIN_SQL).await.unwrap();
+
+ let optimized = df.clone().into_optimized_plan().unwrap();
+ let logical_str = format!("{}", optimized.display_indent());
+ assert!(
+ logical_str.contains("SpatialJoin"),
+ "expected a SpatialJoin extension node in the optimized logical plan, got:\n{logical_str}"
+ );
+
+ let physical = df.create_physical_plan().await.unwrap();
+ let physical_str = displayable(physical.as_ref()).indent(true).to_string();
+ assert!(
+ physical_str.contains("SpatialJoinExec"),
+ "raster join should use the optimized SpatialJoinExec, got:\n{physical_str}"
+ );
+ assert!(
+ !physical_str.contains("NestedLoopJoinExec"),
+ "raster join must not fall back to NestedLoopJoinExec, got:\n{physical_str}"
+ );
+}
+
+/// The raster join executes correctly through the nested-loop fallback: only
+/// (raster 1, point-in-1), (raster 1, big-polygon), and (raster 2, big-polygon)
+/// intersect, and the null raster (rid 0) never matches.
+#[cfg(feature = "proj")]
+#[tokio::test]
+async fn raster_join_produces_correct_rows() {
+ let ctx = SedonaContext::new_local_interactive().await.unwrap();
+ register_raster_table(&ctx);
+ register_geom_table(&ctx);
+
+ let batches = ctx
+ .sql(RASTER_JOIN_SQL)
+ .await
+ .unwrap()
+ .collect()
+ .await
+ .unwrap();
+ datafusion::assert_batches_eq!(
+ [
+ "+-----+-----+",
+ "| rid | gid |",
+ "+-----+-----+",
+ "| 1 | 10 |",
+ "| 1 | 20 |",
+ "| 2 | 20 |",
+ "+-----+-----+",
+ ],
+ &batches
+ );
+}
+
+/// A geometry-only join is unaffected by the raster changes: it still uses the
+/// optimized `SpatialJoinExec` (no regression to the existing ST_* path).
+#[tokio::test]
+async fn geometry_join_still_uses_optimized_spatial_join() {
+ let ctx = SedonaContext::new_local_interactive().await.unwrap();
+ register_geom_table(&ctx);
+
+ // Self-join two geometry tables on ST_Intersects.
+ let geom_type = geom_type();
+ let gid = Int32Array::from(vec![100, 200]);
+ let geom = create_array(
+ &[
+ Some("POLYGON ((0 0, 5 0, 5 5, 0 5, 0 0))"),
+ Some("POINT (100 100)"),
+ ],
+ &geom_type,
+ );
+ let geom_field = geom_type.to_storage_field("geom", true).unwrap();
+ let schema = Arc::new(Schema::new(vec![
+ Field::new("gid2", arrow_schema::DataType::Int32, false),
+ geom_field,
+ ]));
+ let batch = RecordBatch::try_new(schema.clone(), vec![Arc::new(gid), geom]).unwrap();
+ let table = MemTable::try_new(schema, vec![vec![batch]]).unwrap();
+ ctx.ctx.register_table("g2", Arc::new(table)).unwrap();
+
+ let sql = "SELECT g.gid, g2.gid2 \
+ FROM g JOIN g2 ON ST_Intersects(g.geom, g2.geom)";
+ let df = ctx.sql(sql).await.unwrap();
+
+ let physical = df.clone().create_physical_plan().await.unwrap();
+ let physical_str = displayable(physical.as_ref()).indent(true).to_string();
+ assert!(
+ physical_str.contains("SpatialJoinExec"),
+ "geometry join should use the optimized SpatialJoinExec, got:\n{physical_str}"
+ );
+
+ // Sanity: the join still executes and produces at least the expected overlap row.
+ let batches = df.collect().await.unwrap();
+ let total_rows: usize = batches.iter().map(|b| b.num_rows()).sum();
+ assert!(total_rows >= 1, "expected at least one intersecting pair");
+}