How to use geoarrow::algorithm::geo::Intersect

[JosiahParry]

I'm trying to get a bit more familiar with the georrow-rs library. I'd like to be able to test the intersection between two geoarrow scalar geometries. I've tried to use the geoarrow::algorithm::geo::Intersect trait but I am unable to figure out how it works.

Is the intention to cast to geo with to_geo() and then use geo's algorithms in these cases?

For example finding the pairwise intersection between polygon arrays, this is the best i could figure out

use geo::Intersects;
use geoarrow::array::PolygonArray;
use geoarrow::trait_::{GeometryArrayAccessor, GeometryScalarTrait};

pub fn test_intersects(left: PolygonArray<i32>, right: PolygonArray<i32>) -> Vec<bool> {
    left.iter()
        .zip(right.iter())
        .map(|(l, r)| {
            let l = l.unwrap();
            let r = r.unwrap();

            l.to_geo().intersects(&r.to_geo())
        })
        .collect::<Vec<_>>()
}

[kylebarron]

I've tried to use the geoarrow::algorithm::geo::Intersect trait but I am unable to figure out how it works.

Assuming you're on the latest main, that's recently been refactored to use geo-index. There's a new set of Indexed arrays that stores an immutable RTree along with the array data. By defining new structs, we can give the user control over whether they want to use a spatial index with multiple operations. And in the future, since geo-index is ABI stable, we could potentially share an array with its index to Python or R.

You want to first construct an IndexedPolygonArray from each array with IndexedPolygonArray::new, and then you can pass them into the Intersects trait, which will do a pairwise intersection. Note that this uses intersection_candidates_with_other_tree, so it may have bugs right now, but in theory should be very efficient.

Is the intention to cast to geo with to_geo() and then use geo's algorithms in these cases?

Yes. The thought being that the user can control when they want to take a reference onto an existing geoarrow scalar (which is an index into the array and which itself implements geometry traits) and when they want to materialize a geo scalar for algorithms. In practice, every algorithm in geoarrow::algorithm::geo calls to_geo()

For example finding the pairwise intersection between polygon arrays, this is the best i could figure out

There are also helpers for pairwise boolean operations. Look at try_binary_boolean. So you can make your implementation shorter by doing

pub fn test_intersects(
    left_arr: &PolygonArray<i32>,
    right_arr: &PolygonArray<i32>,
) -> std::result::Result<BooleanArray, GeoArrowError> {
    left_arr.try_binary_boolean(right_arr, |left_scalar, right_scalar| {
        Ok(left_scalar.to_geo().intersects(&right_scalar.to_geo()))
    })
}

Additionally, try_binary_boolean will manage the nullability of values, so you don't need to have the unwrap.

binary_boolean and try_binary_boolean are effectively the same when the input arrays have no null values, but when there are null values, binary_boolean will still get run for all values, even when that row is actually null, the thought being that for fast operations, branching is worse than just doing the computation and hopefully getting auto-vectorized. But for something like intersects, you probably want to ensure you're not trying to still run it on null values.

Brute force intersection will also likely be somewhat slow. When the indexed intersection gets stabilized, that should hopefully be very fast.