paleolimbot commented on code in PR #1073: URL: https://github.com/apache/sedona-db/pull/1073#discussion_r3642572841
########## rust/sedona-spatial-join-raster/src/join_provider.rs: ########## @@ -0,0 +1,474 @@ +// 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, JoinType, Result}; +use datafusion_expr::ColumnarValue; +use sedona_common::{sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions}; +use sedona_expr::statistics::GeoStatistics; +use sedona_geometry::{transform::CrsEngine, wkb_factory::write_wkb_polygon}; +use sedona_proj::transform::with_global_proj_engine; +use sedona_raster::array::RasterStructArray; +use sedona_raster::traits::RasterRef; +use sedona_raster_functions::crs_utils::resolve_crs; +use sedona_raster_functions::footprint::{raster_footprint_corners, write_convexhull_wkb}; +use sedona_schema::{ + crs::{lnglat, 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). +#[derive(Debug)] +pub(crate) struct RasterJoinProvider { + default: DefaultSpatialJoinProvider, + target_crs: Crs, +} + +impl RasterJoinProvider { + pub(crate) fn new(target_crs: Crs) -> Self { + Self { + default: DefaultSpatialJoinProvider, + target_crs, + } + } +} + +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, + options, + join_type, + probe_threads_count, + metrics, + ) + } + + fn estimate_extra_memory_usage( + &self, + geo_stats: &GeoStatistics, + spatial_predicate: &SpatialPredicate, + options: &SpatialJoinOptions, + ) -> usize { + self.default + .estimate_extra_memory_usage(geo_stats, spatial_predicate, options) + } + + fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> { + Arc::new(RasterGeometryArrayFactory { + target_crs: self.target_crs.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 convex hull of 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. +/// +/// Cross-CRS raster footprints are indexed and refined as the convex hull of the +/// raster's four reprojected corners; projection curvature along the edges is not +/// modeled, so for large-extent rasters reprojected between very different CRSs +/// the hull can slightly under-cover the true footprint (rare missed matches at +/// the extreme edges). Same-CRS joins are exact. +#[derive(Debug)] +struct RasterGeometryArrayFactory { + /// The geometry operand's CRS: the common CRS this join compares in. Raster + /// footprints are reprojected into it. WGS84 when the geometry operand carries + /// no CRS — a missing CRS is assumed WGS84 (matching the `RS_*` kernel and + /// Sedona Spark), so the CRS-less geometry is compared as-is in a WGS84 frame. + target_crs: Crs, +} + +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); + + // A missing CRS on either side is assumed to be WGS84 (matching the `RS_*` + // kernel and Sedona Spark); resolve it once rather than per row. + let lnglat_crs = lnglat(); + let wgs84 = lnglat_crs.as_deref().expect("lnglat() is always Some"); + + with_global_proj_engine(|engine| { Review Comment: I believe we have access to the SedonaOptions when we create the array factory. If we do we should use that CrsEngine. ########## rust/sedona-spatial-join-raster/src/join_provider.rs: ########## @@ -0,0 +1,474 @@ +// 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, JoinType, Result}; +use datafusion_expr::ColumnarValue; +use sedona_common::{sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions}; +use sedona_expr::statistics::GeoStatistics; +use sedona_geometry::{transform::CrsEngine, wkb_factory::write_wkb_polygon}; +use sedona_proj::transform::with_global_proj_engine; +use sedona_raster::array::RasterStructArray; +use sedona_raster::traits::RasterRef; +use sedona_raster_functions::crs_utils::resolve_crs; +use sedona_raster_functions::footprint::{raster_footprint_corners, write_convexhull_wkb}; +use sedona_schema::{ + crs::{lnglat, 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). +#[derive(Debug)] +pub(crate) struct RasterJoinProvider { + default: DefaultSpatialJoinProvider, + target_crs: Crs, +} + +impl RasterJoinProvider { + pub(crate) fn new(target_crs: Crs) -> Self { + Self { + default: DefaultSpatialJoinProvider, + target_crs, + } + } +} + +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, + options, + join_type, + probe_threads_count, + metrics, + ) + } + + fn estimate_extra_memory_usage( + &self, + geo_stats: &GeoStatistics, + spatial_predicate: &SpatialPredicate, + options: &SpatialJoinOptions, + ) -> usize { + self.default + .estimate_extra_memory_usage(geo_stats, spatial_predicate, options) + } + + fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> { + Arc::new(RasterGeometryArrayFactory { + target_crs: self.target_crs.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 convex hull of 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. +/// +/// Cross-CRS raster footprints are indexed and refined as the convex hull of the +/// raster's four reprojected corners; projection curvature along the edges is not +/// modeled, so for large-extent rasters reprojected between very different CRSs +/// the hull can slightly under-cover the true footprint (rare missed matches at +/// the extreme edges). Same-CRS joins are exact. +#[derive(Debug)] +struct RasterGeometryArrayFactory { + /// The geometry operand's CRS: the common CRS this join compares in. Raster + /// footprints are reprojected into it. WGS84 when the geometry operand carries + /// no CRS — a missing CRS is assumed WGS84 (matching the `RS_*` kernel and + /// Sedona Spark), so the CRS-less geometry is compared as-is in a WGS84 frame. + target_crs: Crs, +} + +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); + + // A missing CRS on either side is assumed to be WGS84 (matching the `RS_*` + // kernel and Sedona Spark); resolve it once rather than per row. + let lnglat_crs = lnglat(); + let wgs84 = lnglat_crs.as_deref().expect("lnglat() is always Some"); + + with_global_proj_engine(|engine| { + 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, wgs84, engine, &mut builder)?; + rects.push(rect); + } + Ok(()) + })?; + + 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. `wgs84` is + /// the WGS84 CRS a missing side is assumed to be. + /// + /// CRS rules mirror the `RS_*` kernel: a missing CRS on either side is assumed + /// to be WGS84 (matching Sedona Spark), after which an equal CRS compares + /// directly and a genuine CRS difference reprojects the footprint's four + /// corners into the target CRS. The footprint is always the convex hull of the + /// four corners — projection curvature along the edges is not modeled (see the + /// type-level note on [`RasterGeometryArrayFactory`]). + fn append_footprint( + &self, + raster: &dyn RasterRef, + wgs84: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, + ) -> Result<Bounds2D> { + let raster_crs = resolve_crs(raster.crs())?; + let corners = raster_footprint_corners(raster); + + // Assume WGS84 for whichever side lacks a CRS, then compare in the target + // frame. Both-absent and same-CRS collapse to a direct comparison. + let raster_crs = raster_crs.as_deref().unwrap_or(wgs84); + let target_crs = self.target_crs.as_deref().unwrap_or(wgs84); Review Comment: If this is the current behaviour for any version of an RS predicate, we should change it. I think this is probably not as relevant in Sedona Spark since the GeoParquet reader properly attached CRSes (when there would have been a lot of crs-less values floating around that actually were WGS84). ########## rust/sedona-raster-functions/src/rs_spatial_predicates.rs: ########## @@ -318,24 +316,54 @@ fn evaluate_predicate_with_crs<Op: tg::BinaryPredicate>( from_a_to_b: bool, engine: &dyn CrsEngine, ) -> Result<bool> { - // If either side has no CRS, compare directly without transformation. - let (crs_a, crs_b) = match (crs_a, crs_b) { - (Some(a), Some(b)) => (a, b), - (None, None) => return evaluate_predicate::<Op>(wkb_a, wkb_b), - (Some(_), None) => { - return exec_err!( - "Cannot evaluate spatial predicate: \ - left geometry has CRS but right geometry does not" + match (crs_a, crs_b) { + (Some(crs_a), Some(crs_b)) => { + compare_in_crs::<Op>(wkb_a, crs_a, wkb_b, crs_b, from_a_to_b, engine) + } + // Neither side has a CRS: both are in the same (unknown) frame, so compare + // directly without transformation. + (None, None) => evaluate_predicate::<Op>(wkb_a, wkb_b), + // Exactly one side has a CRS: assume the missing side is WGS84 (matching + // Sedona Spark) and compare in the resulting CRS pair. Only these arms + // materialize the lnglat CRS, so the common paths avoid the allocation. + (Some(crs_a), None) => { + let lnglat_crs = lnglat().expect("lnglat() should always return Some"); Review Comment: The assumption of a missing side having a WGS84 CRS is a footgun and the current behaviour is better. We can fix this in Sedona Spark now that the GeoParquet reader attaches CRSes. ########## rust/sedona-raster-functions/src/footprint.rs: ########## @@ -0,0 +1,58 @@ +// 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::wkb_factory::write_wkb_polygon; +use sedona_raster::affine_transformation::to_world_coordinate; +use sedona_raster::traits::RasterRef; + +/// 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 of the raster footprint. +/// +/// The ring is the four [`raster_footprint_corners`] closed back to the +/// upper-left corner. This can be used to build Binary arrays, as the arrow-rs +/// `BinaryBuilder` implements [`std::io::Write`]. +pub fn write_convexhull_wkb(raster: &dyn RasterRef, out: &mut impl std::io::Write) -> Result<()> { Review Comment: This should densify the edges at least a little bit (around 10 points per edge is common) to ensure that when projected into another geometry's space, the calculated intersection is more likely to be accurate (many transformations turn straight lines into curved ones). The logic to do this is already in the repo for `st_segmentize`...you can move any utilities to sedona-geometry that are needed here. A more robustly parameterized method would transform the convex hull and the geometry into a geography (with a tolerance in meters), then compare it. If somebody needs that level of accuracy they can write a join that does that fairly easily. ########## rust/sedona-spatial-join-raster/src/join_provider.rs: ########## @@ -0,0 +1,474 @@ +// 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, JoinType, Result}; +use datafusion_expr::ColumnarValue; +use sedona_common::{sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions}; +use sedona_expr::statistics::GeoStatistics; +use sedona_geometry::{transform::CrsEngine, wkb_factory::write_wkb_polygon}; +use sedona_proj::transform::with_global_proj_engine; +use sedona_raster::array::RasterStructArray; +use sedona_raster::traits::RasterRef; +use sedona_raster_functions::crs_utils::resolve_crs; +use sedona_raster_functions::footprint::{raster_footprint_corners, write_convexhull_wkb}; +use sedona_schema::{ + crs::{lnglat, 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). +#[derive(Debug)] +pub(crate) struct RasterJoinProvider { + default: DefaultSpatialJoinProvider, + target_crs: Crs, +} + +impl RasterJoinProvider { + pub(crate) fn new(target_crs: Crs) -> Self { + Self { + default: DefaultSpatialJoinProvider, + target_crs, + } + } +} + +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, + options, + join_type, + probe_threads_count, + metrics, + ) + } + + fn estimate_extra_memory_usage( + &self, + geo_stats: &GeoStatistics, + spatial_predicate: &SpatialPredicate, + options: &SpatialJoinOptions, + ) -> usize { + self.default + .estimate_extra_memory_usage(geo_stats, spatial_predicate, options) + } + + fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> { + Arc::new(RasterGeometryArrayFactory { + target_crs: self.target_crs.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 convex hull of 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. +/// +/// Cross-CRS raster footprints are indexed and refined as the convex hull of the +/// raster's four reprojected corners; projection curvature along the edges is not +/// modeled, so for large-extent rasters reprojected between very different CRSs +/// the hull can slightly under-cover the true footprint (rare missed matches at +/// the extreme edges). Same-CRS joins are exact. +#[derive(Debug)] +struct RasterGeometryArrayFactory { + /// The geometry operand's CRS: the common CRS this join compares in. Raster + /// footprints are reprojected into it. WGS84 when the geometry operand carries + /// no CRS — a missing CRS is assumed WGS84 (matching the `RS_*` kernel and + /// Sedona Spark), so the CRS-less geometry is compared as-is in a WGS84 frame. + target_crs: Crs, +} + +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); + + // A missing CRS on either side is assumed to be WGS84 (matching the `RS_*` + // kernel and Sedona Spark); resolve it once rather than per row. + let lnglat_crs = lnglat(); + let wgs84 = lnglat_crs.as_deref().expect("lnglat() is always Some"); + + with_global_proj_engine(|engine| { + 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, wgs84, engine, &mut builder)?; + rects.push(rect); + } + Ok(()) + })?; + + 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. `wgs84` is + /// the WGS84 CRS a missing side is assumed to be. + /// + /// CRS rules mirror the `RS_*` kernel: a missing CRS on either side is assumed + /// to be WGS84 (matching Sedona Spark), after which an equal CRS compares + /// directly and a genuine CRS difference reprojects the footprint's four + /// corners into the target CRS. The footprint is always the convex hull of the + /// four corners — projection curvature along the edges is not modeled (see the + /// type-level note on [`RasterGeometryArrayFactory`]). + fn append_footprint( + &self, + raster: &dyn RasterRef, + wgs84: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, + ) -> Result<Bounds2D> { + let raster_crs = resolve_crs(raster.crs())?; + let corners = raster_footprint_corners(raster); + + // Assume WGS84 for whichever side lacks a CRS, then compare in the target + // frame. Both-absent and same-CRS collapse to a direct comparison. + let raster_crs = raster_crs.as_deref().unwrap_or(wgs84); + let target_crs = self.target_crs.as_deref().unwrap_or(wgs84); + + if raster_crs.crs_equals(target_crs) { + // Same (possibly assumed-WGS84) CRS: compare directly. The footprint is + // byte-identical to the one the `RS_*` kernel builds. + write_convexhull_wkb(raster, builder)?; + builder.append_value([]); + Ok(bounds_from_coords(&corners)) + } else { + // Genuine CRS difference: reproject the four corners into the target CRS + // and emit their convex hull. + append_reprojected_footprint(corners, raster_crs, target_crs, engine, builder) + } + } +} + +/// Reproject a raster's four `corners` from `from_crs` to `to_crs`, append the +/// convex hull of the reprojected corners to `builder`, and return its bounding +/// rectangle. +fn append_reprojected_footprint( + corners: [(f64, f64); 4], + from_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + to_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, +) -> Result<Bounds2D> { + let transform = engine + .get_transform_crs_to_crs(&from_crs.to_crs_string(), &to_crs.to_crs_string(), None, "") + .map_err(|e| exec_datafusion_err!("CRS transform error: {e}"))?; + + let mut hull = corners; + for corner in hull.iter_mut() { + transform + .transform_coord(corner) + .map_err(|e| exec_datafusion_err!("Transform error: {e}"))?; + } + + // Closed ring: the four reprojected corners plus the first repeated. + write_wkb_polygon( + builder, + [hull[0], hull[1], hull[2], hull[3], hull[0]].into_iter(), + ) + .map_err(|e| exec_datafusion_err!("Failed to write footprint WKB: {e}"))?; Review Comment: Is there another utility in this PR that does this? There may also be a utility that does this already in wkb_factory. ########## rust/sedona-spatial-join-raster/src/join_provider.rs: ########## @@ -0,0 +1,474 @@ +// 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, JoinType, Result}; +use datafusion_expr::ColumnarValue; +use sedona_common::{sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions}; +use sedona_expr::statistics::GeoStatistics; +use sedona_geometry::{transform::CrsEngine, wkb_factory::write_wkb_polygon}; +use sedona_proj::transform::with_global_proj_engine; +use sedona_raster::array::RasterStructArray; +use sedona_raster::traits::RasterRef; +use sedona_raster_functions::crs_utils::resolve_crs; +use sedona_raster_functions::footprint::{raster_footprint_corners, write_convexhull_wkb}; +use sedona_schema::{ + crs::{lnglat, 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). +#[derive(Debug)] +pub(crate) struct RasterJoinProvider { + default: DefaultSpatialJoinProvider, + target_crs: Crs, +} + +impl RasterJoinProvider { + pub(crate) fn new(target_crs: Crs) -> Self { + Self { + default: DefaultSpatialJoinProvider, + target_crs, + } + } +} + +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, + options, + join_type, + probe_threads_count, + metrics, + ) + } + + fn estimate_extra_memory_usage( + &self, + geo_stats: &GeoStatistics, + spatial_predicate: &SpatialPredicate, + options: &SpatialJoinOptions, + ) -> usize { + self.default + .estimate_extra_memory_usage(geo_stats, spatial_predicate, options) + } + + fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> { + Arc::new(RasterGeometryArrayFactory { + target_crs: self.target_crs.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 convex hull of 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. +/// +/// Cross-CRS raster footprints are indexed and refined as the convex hull of the +/// raster's four reprojected corners; projection curvature along the edges is not +/// modeled, so for large-extent rasters reprojected between very different CRSs +/// the hull can slightly under-cover the true footprint (rare missed matches at +/// the extreme edges). Same-CRS joins are exact. +#[derive(Debug)] +struct RasterGeometryArrayFactory { + /// The geometry operand's CRS: the common CRS this join compares in. Raster + /// footprints are reprojected into it. WGS84 when the geometry operand carries + /// no CRS — a missing CRS is assumed WGS84 (matching the `RS_*` kernel and + /// Sedona Spark), so the CRS-less geometry is compared as-is in a WGS84 frame. + target_crs: Crs, +} + +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); + + // A missing CRS on either side is assumed to be WGS84 (matching the `RS_*` + // kernel and Sedona Spark); resolve it once rather than per row. + let lnglat_crs = lnglat(); + let wgs84 = lnglat_crs.as_deref().expect("lnglat() is always Some"); + + with_global_proj_engine(|engine| { + 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, wgs84, engine, &mut builder)?; + rects.push(rect); + } + Ok(()) + })?; + + 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. `wgs84` is + /// the WGS84 CRS a missing side is assumed to be. + /// + /// CRS rules mirror the `RS_*` kernel: a missing CRS on either side is assumed + /// to be WGS84 (matching Sedona Spark), after which an equal CRS compares + /// directly and a genuine CRS difference reprojects the footprint's four + /// corners into the target CRS. The footprint is always the convex hull of the + /// four corners — projection curvature along the edges is not modeled (see the + /// type-level note on [`RasterGeometryArrayFactory`]). + fn append_footprint( + &self, + raster: &dyn RasterRef, + wgs84: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, + ) -> Result<Bounds2D> { + let raster_crs = resolve_crs(raster.crs())?; + let corners = raster_footprint_corners(raster); + + // Assume WGS84 for whichever side lacks a CRS, then compare in the target + // frame. Both-absent and same-CRS collapse to a direct comparison. + let raster_crs = raster_crs.as_deref().unwrap_or(wgs84); + let target_crs = self.target_crs.as_deref().unwrap_or(wgs84); + + if raster_crs.crs_equals(target_crs) { + // Same (possibly assumed-WGS84) CRS: compare directly. The footprint is + // byte-identical to the one the `RS_*` kernel builds. + write_convexhull_wkb(raster, builder)?; + builder.append_value([]); + Ok(bounds_from_coords(&corners)) + } else { + // Genuine CRS difference: reproject the four corners into the target CRS + // and emit their convex hull. + append_reprojected_footprint(corners, raster_crs, target_crs, engine, builder) + } + } +} + +/// Reproject a raster's four `corners` from `from_crs` to `to_crs`, append the +/// convex hull of the reprojected corners to `builder`, and return its bounding +/// rectangle. +fn append_reprojected_footprint( + corners: [(f64, f64); 4], + from_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + to_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, +) -> Result<Bounds2D> { + let transform = engine + .get_transform_crs_to_crs(&from_crs.to_crs_string(), &to_crs.to_crs_string(), None, "") + .map_err(|e| exec_datafusion_err!("CRS transform error: {e}"))?; + + let mut hull = corners; + for corner in hull.iter_mut() { + transform + .transform_coord(corner) + .map_err(|e| exec_datafusion_err!("Transform error: {e}"))?; + } + + // Closed ring: the four reprojected corners plus the first repeated. + write_wkb_polygon( + builder, + [hull[0], hull[1], hull[2], hull[3], hull[0]].into_iter(), + ) + .map_err(|e| exec_datafusion_err!("Failed to write footprint WKB: {e}"))?; + builder.append_value([]); + Ok(bounds_from_coords(&hull)) +} + +/// Bounding rectangle of a set of coordinates. [`Bounds2D::new`] enlarges the +/// f32 bounds outward so the rectangle conservatively contains every input +/// coordinate. +fn bounds_from_coords(coords: &[(f64, f64)]) -> Bounds2D { Review Comment: Unless this is a benchmarked performance optimization, using the existing utilities (Interval, WkbBounder2D, other lower level bounding utilities) would be better for this ########## rust/sedona-spatial-join-raster/src/physical_planner.rs: ########## @@ -0,0 +1,322 @@ +// 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::physical_plan::ExecutionPlan; +use datafusion_common::Result; +use datafusion_physical_expr::PhysicalExpr; +use sedona_query_planner::{ + spatial_join_physical_planner::{PlanSpatialJoinArgs, SpatialJoinPhysicalPlanner}, + spatial_predicate::{RelationPredicate, SpatialPredicate, SpatialRelationType}, +}; +use sedona_schema::{ + crs::{lnglat, 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()) + }; Review Comment: Do we want to do anything here given some knowledge of which side is the raster side and which side is the vector side? The raster side will almost always be smaller and have fewer vertices. I would guess that it would always be better to index the geometry side but I don't know this. ########## rust/sedona-spatial-join-raster/src/join_provider.rs: ########## @@ -0,0 +1,474 @@ +// 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, JoinType, Result}; +use datafusion_expr::ColumnarValue; +use sedona_common::{sedona_internal_datafusion_err, sedona_internal_err, SpatialJoinOptions}; +use sedona_expr::statistics::GeoStatistics; +use sedona_geometry::{transform::CrsEngine, wkb_factory::write_wkb_polygon}; +use sedona_proj::transform::with_global_proj_engine; +use sedona_raster::array::RasterStructArray; +use sedona_raster::traits::RasterRef; +use sedona_raster_functions::crs_utils::resolve_crs; +use sedona_raster_functions::footprint::{raster_footprint_corners, write_convexhull_wkb}; +use sedona_schema::{ + crs::{lnglat, 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). +#[derive(Debug)] +pub(crate) struct RasterJoinProvider { + default: DefaultSpatialJoinProvider, + target_crs: Crs, +} + +impl RasterJoinProvider { + pub(crate) fn new(target_crs: Crs) -> Self { + Self { + default: DefaultSpatialJoinProvider, + target_crs, + } + } +} + +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, + options, + join_type, + probe_threads_count, + metrics, + ) + } + + fn estimate_extra_memory_usage( + &self, + geo_stats: &GeoStatistics, + spatial_predicate: &SpatialPredicate, + options: &SpatialJoinOptions, + ) -> usize { + self.default + .estimate_extra_memory_usage(geo_stats, spatial_predicate, options) + } + + fn evaluated_array_factory(&self) -> Arc<dyn EvaluatedGeometryArrayFactory> { + Arc::new(RasterGeometryArrayFactory { + target_crs: self.target_crs.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 convex hull of 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. +/// +/// Cross-CRS raster footprints are indexed and refined as the convex hull of the +/// raster's four reprojected corners; projection curvature along the edges is not +/// modeled, so for large-extent rasters reprojected between very different CRSs +/// the hull can slightly under-cover the true footprint (rare missed matches at +/// the extreme edges). Same-CRS joins are exact. +#[derive(Debug)] +struct RasterGeometryArrayFactory { + /// The geometry operand's CRS: the common CRS this join compares in. Raster + /// footprints are reprojected into it. WGS84 when the geometry operand carries + /// no CRS — a missing CRS is assumed WGS84 (matching the `RS_*` kernel and + /// Sedona Spark), so the CRS-less geometry is compared as-is in a WGS84 frame. + target_crs: Crs, +} + +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); + + // A missing CRS on either side is assumed to be WGS84 (matching the `RS_*` + // kernel and Sedona Spark); resolve it once rather than per row. + let lnglat_crs = lnglat(); + let wgs84 = lnglat_crs.as_deref().expect("lnglat() is always Some"); + + with_global_proj_engine(|engine| { + 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, wgs84, engine, &mut builder)?; + rects.push(rect); + } + Ok(()) + })?; + + 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. `wgs84` is + /// the WGS84 CRS a missing side is assumed to be. + /// + /// CRS rules mirror the `RS_*` kernel: a missing CRS on either side is assumed + /// to be WGS84 (matching Sedona Spark), after which an equal CRS compares + /// directly and a genuine CRS difference reprojects the footprint's four + /// corners into the target CRS. The footprint is always the convex hull of the + /// four corners — projection curvature along the edges is not modeled (see the + /// type-level note on [`RasterGeometryArrayFactory`]). + fn append_footprint( + &self, + raster: &dyn RasterRef, + wgs84: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, + ) -> Result<Bounds2D> { + let raster_crs = resolve_crs(raster.crs())?; + let corners = raster_footprint_corners(raster); + + // Assume WGS84 for whichever side lacks a CRS, then compare in the target + // frame. Both-absent and same-CRS collapse to a direct comparison. + let raster_crs = raster_crs.as_deref().unwrap_or(wgs84); + let target_crs = self.target_crs.as_deref().unwrap_or(wgs84); + + if raster_crs.crs_equals(target_crs) { + // Same (possibly assumed-WGS84) CRS: compare directly. The footprint is + // byte-identical to the one the `RS_*` kernel builds. + write_convexhull_wkb(raster, builder)?; + builder.append_value([]); + Ok(bounds_from_coords(&corners)) + } else { + // Genuine CRS difference: reproject the four corners into the target CRS + // and emit their convex hull. + append_reprojected_footprint(corners, raster_crs, target_crs, engine, builder) + } + } +} + +/// Reproject a raster's four `corners` from `from_crs` to `to_crs`, append the +/// convex hull of the reprojected corners to `builder`, and return its bounding +/// rectangle. +fn append_reprojected_footprint( + corners: [(f64, f64); 4], + from_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + to_crs: &(dyn CoordinateReferenceSystem + Send + Sync), + engine: &dyn CrsEngine, + builder: &mut BinaryBuilder, +) -> Result<Bounds2D> { + let transform = engine + .get_transform_crs_to_crs(&from_crs.to_crs_string(), &to_crs.to_crs_string(), None, "") + .map_err(|e| exec_datafusion_err!("CRS transform error: {e}"))?; + + let mut hull = corners; + for corner in hull.iter_mut() { + transform + .transform_coord(corner) + .map_err(|e| exec_datafusion_err!("Transform error: {e}"))?; + } + + // Closed ring: the four reprojected corners plus the first repeated. + write_wkb_polygon( + builder, + [hull[0], hull[1], hull[2], hull[3], hull[0]].into_iter(), + ) + .map_err(|e| exec_datafusion_err!("Failed to write footprint WKB: {e}"))?; + builder.append_value([]); + Ok(bounds_from_coords(&hull)) +} + +/// Bounding rectangle of a set of coordinates. [`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 min_x = f64::INFINITY; + let mut max_x = f64::NEG_INFINITY; + let mut min_y = f64::INFINITY; + let mut max_y = f64::NEG_INFINITY; + for &(x, y) in coords { + min_x = min_x.min(x); + max_x = max_x.max(x); + min_y = min_y.min(y); + max_y = max_y.max(y); + } + Bounds2D::new((min_x, max_x), (min_y, max_y)) +} + +#[cfg(test)] +mod test { + use super::*; + + #[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 arrow_array::BinaryArray; + use sedona_raster::builder::RasterBuilder; + use sedona_raster::traits::{BandMetadata, RasterMetadata}; + use sedona_schema::crs::{deserialize_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 { Review Comment: Does the RasterSpec help this at all? -- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. To unsubscribe, e-mail: [email protected] For queries about this service, please contact Infrastructure at: [email protected]
