feat(vortex-geo): add intersection scalar function

vortex-geo/Cargo.toml

@@ -78,5 +78,9 @@ harness = false
 name = "convex_hull"
 harness = false
 
+[[bench]]
+name = "intersection"
+harness = false
+
 [lints]
 workspace = true

vortex-geo/benches/intersection.rs

@@ -0,0 +1,107 @@
+// SPDX-License-Identifier: Apache-2.0
+// SPDX-FileCopyrightText: Copyright the Vortex contributors
+
+//! Microbenchmarks for native `ST_Intersection` over polygon pairs.
+//!
+//! The cases cover simple building-like rectangles, more detailed boundaries, and strict null
+//! propagation. Inputs overlap because SpatialBench Q9 prefilters pairs with `ST_Intersects`.
+//!
+//! Run with `cargo bench -p vortex-geo --bench intersection`.
+
+#![expect(clippy::unwrap_used)]
+
+use std::f64::consts::TAU;
+use std::sync::LazyLock;
+
+use divan::Bencher;
+use divan::counter::ItemsCount;
+use mimalloc::MiMalloc;
+use vortex_array::ArrayRef;
+use vortex_array::Canonical;
+use vortex_array::ExecutionCtx;
+use vortex_array::IntoArray;
+use vortex_array::VortexSessionExecute;
+use vortex_array::arrays::MaskedArray;
+use vortex_array::validity::Validity;
+use vortex_geo::scalar_fn::intersection::GeoIntersection;
+use vortex_geo::test_harness::geo_session;
+use vortex_geo::test_harness::polygon_column;
+use vortex_session::VortexSession;
+
+// Scalar function execution allocates its output inside the timed region, so use the vendored
+// allocator instead of measuring glibc differences between CodSpeed runner images.
+#[global_allocator]
+static GLOBAL: MiMalloc = MiMalloc;
+
+static SESSION: LazyLock<VortexSession> = LazyLock::new(geo_session);
+
+const ROWS: usize = 512;
+
+fn main() {
+    divan::main();
+}
+
+fn regular_polygon(cx: f64, cy: f64, radius: f64, vertices: usize) -> Vec<(f64, f64)> {
+    (0..=vertices)
+        .map(|vertex| {
+            let angle = TAU * (vertex % vertices) as f64 / vertices as f64;
+            (cx + radius * angle.cos(), cy + radius * angle.sin())
+        })
+        .collect()
+}
+
+fn polygon_pairs(vertices: usize) -> (ArrayRef, ArrayRef) {
+    let left = polygon_column(
+        (0..ROWS)
+            .map(|row| vec![regular_polygon(row as f64, 0.0, 1.0, vertices)])
+            .collect(),
+    )
+    .unwrap();
+    let right = polygon_column(
+        (0..ROWS)
+            .map(|row| vec![regular_polygon(row as f64 + 0.5, 0.0, 1.0, vertices)])
+            .collect(),
+    )
+    .unwrap();
+    (left, right)
+}
+
+fn intersections(left: &ArrayRef, right: &ArrayRef, ctx: &mut ExecutionCtx) -> ArrayRef {
+    GeoIntersection::try_new_array(left.clone(), right.clone())
+        .unwrap()
+        .into_array()
+        .execute::<Canonical>(ctx)
+        .unwrap()
+        .into_array()
+}
+
+fn bench_intersections(bencher: Bencher, left: ArrayRef, right: ArrayRef) {
+    let mut ctx = SESSION.create_execution_ctx();
+    bencher
+        .counter(ItemsCount::new(ROWS))
+        .bench_local(|| intersections(&left, &right, &mut ctx));
+}
+
+#[divan::bench]
+fn rectangles(bencher: Bencher) {
+    let (left, right) = polygon_pairs(4);
+    bench_intersections(bencher, left, right);
+}
+
+#[divan::bench]
+fn thirty_two_vertex_boundaries(bencher: Bencher) {
+    let (left, right) = polygon_pairs(32);
+    bench_intersections(bencher, left, right);
+}
+
+#[divan::bench]
+fn nullable_rectangles(bencher: Bencher) {
+    let (left, right) = polygon_pairs(4);
+    let left = MaskedArray::try_new(
+        left,
+        Validity::from_iter((0..ROWS).map(|row| !row.is_multiple_of(8))),
+    )
+    .unwrap()
+    .into_array();
+    bench_intersections(bencher, left, right);
+}

vortex-geo/src/extension/multipolygon.rs

@@ -13,9 +13,11 @@ use arrow_schema::Field;
 use arrow_schema::extension::ExtensionType;
 use geo_traits::to_geo::ToGeoGeometry;
 use geo_types::Geometry;
+use geoarrow::array::GeoArrowArray;
 use geoarrow::array::GeoArrowArrayAccessor;
 use geoarrow::array::IntoArrow;
 use geoarrow::array::MultiPolygonArray;
+use geoarrow::array::MultiPolygonBuilder;
 use geoarrow::datatypes::CoordType;
 use geoarrow::datatypes::MultiPolygonType;
 use prost::Message;
@@ -24,6 +26,7 @@ use vortex_array::ExecutionCtx;
 use vortex_array::IntoArray;
 use vortex_array::arrays::ExtensionArray;
 use vortex_array::arrays::extension::ExtensionArrayExt;
+use vortex_array::builtins::ArrayBuiltins;
 use vortex_array::dtype::DType;
 use vortex_array::dtype::Nullability;
 use vortex_array::dtype::extension::ExtDType;
@@ -117,6 +120,27 @@ fn multipolygon_type(geo_metadata: &GeoMetadata, dimension: Dimension) -> MultiP
     MultiPolygonType::new(dimension.into(), geoarrow_metadata(geo_metadata))
 }
 
+/// Build a native 2-D [`MultiPolygon`] array from row-oriented `geo_types` multipolygons.
+pub(crate) fn build_multipolygon_array(
+    multipolygons: &[Option<geo_types::MultiPolygon<f64>>],
+    metadata: GeoMetadata,
+    nullability: Nullability,
+) -> VortexResult<ArrayRef> {
+    let multipolygons = MultiPolygonBuilder::from_nullable_multi_polygons(
+        multipolygons,
+        multipolygon_type(&metadata, Dimension::Xy),
+    )
+    .finish();
+    let storage_dtype = multipolygon_storage_dtype(Dimension::Xy, nullability);
+    let storage = ArrayRef::from_arrow(
+        multipolygons.to_array_ref().as_ref(),
+        nullability == Nullability::Nullable,
+    )?
+    .cast(storage_dtype.clone())?;
+    let ext_dtype = ExtDType::<MultiPolygon>::try_new(metadata, storage_dtype)?;
+    Ok(ExtensionArray::try_new(ext_dtype.erased(), storage)?.into_array())
+}
+
 /// Decode storage to `geo_types` for the geo scalar functions (CRS is irrelevant to planar ops).
 pub(crate) fn multipolygon_geometries(
     storage: &ArrayRef,

vortex-geo/src/lib.rs

@@ -27,6 +27,7 @@ use crate::scalar_fn::contains::GeoContains;
 use crate::scalar_fn::convex_hull::GeoConvexHull;
 use crate::scalar_fn::distance::GeoDistance;
 use crate::scalar_fn::envelope::GeoEnvelope;
+use crate::scalar_fn::intersection::GeoIntersection;
 use crate::scalar_fn::intersects::GeoIntersects;
 use crate::scalar_fn::length::GeoLength;
 use crate::scalar_fn::make_line::GeoMakeLine;
@@ -76,6 +77,7 @@ pub fn initialize(session: &VortexSession) {
     session.scalar_fns().register(GeoContains);
     session.scalar_fns().register(GeoDistance);
     session.scalar_fns().register(GeoIntersects);
+    session.scalar_fns().register(GeoIntersection);
     session.scalar_fns().register(GeoLength);
     session.scalar_fns().register(GeoMakeLine);