Trait-based Geometries

[urschrei]

There's been some recent discussion around implementing geometry primitives in terms of traits. The idea isn't new (see e.g. #67), and has seen some concrete use in https://github.com/andelf/rust-postgis/blob/master/src/types.rs for instance, but we haven't pursued it in geo-types as yet.

The need for efficient (zero-copy) access to the geometry representations in other formats (e.g. Arrow) has brought the idea to the fore again.

This would require a significant overhaul of many of geos APIs, and we're nowhere near a concrete idea of what those changes would look like. For the moment, this discussion serves as a repository for ideas, sketches, experiments, with the possible goal of producing an RFC for discussion.

[michaelkirk]

This would require a significant overhaul of many of geos APIs

It seems like trait-based geo algorithms could be developed in a separate repository or crate. If it pans out, there might be a reason to have geo depend on it, or merge it into georust/geo if it seems beneficial.

This might help avoid a situation where we're rewriting (and frankly, probably complicating) already working geo code under the auspices of "maybe this will be useful to someone someday".

This is similar to what we did with Proj - which used to have a hard dependency on geo-types, but now we have a proj::Coord trait that represents the data type and a proj::Transform trait that represents the algorithm. In theory other people can rely on these Proj traits for their geometry types (though I don't know if anyone is).

[urschrei]

This is similar to what we did with Proj - which used to have a hard dependency on geo-types, but now we have a proj::Coord trait that represents the data type and a proj::Transform trait that represents the algorithm.

I had a similar thought and am definitely in favour of this approach.

[kylebarron]

It's taken me like 4 months of learning, but I think I'm finally to a point where I can contribute to the discussion 😄.

Motivation

I'm really excited about cross-language tools and workflows, where it's easy to pass data between Rust to Python to C++, and similarly in the browser between WebAssembly and JS. Today, there tends to be data silos where code was written for one memory model, and thus in order to use that data from another program, it needs to be serialized and deserialized.

In the geo world, that data silo tends to be GEOS/JTS, because so much of today's geo tools have been implemented around them. But with the advent of FlatGeobuf and GeoArrow, I think there's potential in having the memory model as the core specification, not the library.

I think the Georust community is really well poised to take advantage of this going forward because:

• Algorithms are implemented in pure Rust, meaning that it's easy to bind/compile to various targets and that algorithms are thread safe by default (I think...) for use in e.g. rayon
• Rust's trait system enables zero-cost abstractions over the actual memory layout

Potential zero-copy memory formats

GeoArrow

Apache Arrow is a columnar data format designed for cross-language support and zero-copy access. GeoArrow is an in-progress specification to extend Arrow to handle geospatial vector data. When using the arrow-native coordinate encoding, any coordinate of any geometry in the array can be accessed in constant time.

Arrow has already seen wide use outside of the Geo sphere, but I think there's a ton of potential for it in the geo world as well. I do a lot of geo data visualization in the browser, and I've written about how deck.gl (a WebGL-enabled map rendering library) can use GeoArrow data directly, saving a ton of CPU processing time. So you could imagine a data exploration workflow like:

georust algorithms -->
  Python memory -->
    JS memory (e.g. inside Jupyter Notebooks, one copy, but no serde) -->
      GPU memory for visualization

Being able to go through this entire process where the memory format is virtually unchanged in each step would be super fast and reduce memory overheads.

FlatGeobuf

FlatGeobuf is a performance-oriented vector geometry format. Until now, FlatGeobuf has been used almost exclusively as a file format, where language support for FlatGeobuf means loading the file to convert it into that language's native geometry type. But Flatbuffers, which FlatGeobuf is built upon, allows for zero-copy access. So FlatGeobuf could also be used as a memory format. Indeed, I believe this is why FlatGeobuf does not offer compression, so that it could be memory mapped. But if you can't operate on the original memory layout, then memory mapping loses some of its value.

GEOS?

I don't know the internals of GEOS well, but it looks like it is theoretically possible to read GEOS memory without a copy. Then georust could share memory with e.g. shapely in Python.

Trait Definition

I've been hacking on a trait definition on this branch, as of today at commit 0fb9433. (Due to Rust's trait orphan rules, it seems easier to hack on this on a branch in geo than in a separate repo). The trait implementation so far is in the traits folder.

Trait methods

One big question would be how many methods to require on the traits themselves, vs implementing on top of the traits. Fewer methods on the traits would of course make them easier to implement but harder to use. For now, I just define the most minimal methods to be able to access underlying data.

Geometry types

For now I only implemented the 6 core types: (Multi)Point, (Multi)LineString, and (Multi)Polygon. At least GeoArrow doesn't define other types like Rect or Triangle or Line. Down the road (if there's momentum behind traits), it would probably make sense to implement traits for them as well.

Ownership

I at first planned to exclusively use references in the traits, but ended up hitting lots of

cannot return value referencing local variable

when trying to pass in Arrow data. As I understand it, Arrow2 is designed to support copy-on-write semantics, so clones are really cheap and using owned objects in the traits seemed to make things easier for now, but this is certainly something to revisit.

Generic f32/f64

For now, for simplicity, I assume f64 in the traits, but this is something that I would expect to be relaxed in the future.

Geo function implementation

The biggest hurdle I faced early in geo was that all existing algorithms are defined as traits themselves. It's not clear to me today (especially without trait specialization) if there's a way to define a trait that takes any implementation of another trait as an argument, at least not without resorting to trait objects and dynamic dispatch.

As an aside: I don't really know what the overhead of dynamic dispatch would be, but I assume it could be significant. Especially in situations like GeoArrow and FlatGeobuf, often the metadata would declare that all geometries are a single geometry type, and so the downcasting on each geometry would be unnecessary overhead.

For now, in an effort to get something compiling, I resorted to creating a custom top-level area_polygon function which takes geom: &'a impl PolygonTrait<'a> as an argument.

It did hint at the level of work required: many functions access struct data directly, and so updating all geo algorithms to use only the methods defined on the traits would be a significant undertaking.

Writing

Theoretically you could imagine some way of also writing data from a function to a specific output format, but for now I'll consider this out of scope.

Trait Implementations

GeoArrow

The GeoArrow spec defines how to store a column of geometries in Arrow memory. The existing format spec is defined here, though there's also discussion about changing that to use a StructArray approach. Since Polars doesn't today support the Arrow FixedSizeList type, I went ahead and prototyped the StructArray approach. (They're very similar, just a question of whether the internal coordinate array is interleaved or separated.) I also try to give an overview of what logical/physical representations of GeoArrow memory are in a blog post here.

I have a WIP GeoArrow implementation of these traits here (permalink).

The idea is that a user will have a column of singly-typed memory, say a column of polygons. polygon_index will index into the input array and return a single polygon element ArrowPolygonScalar. The underlying data isn't copied, just the pointer on the Box. Then on that ArrowPolygonScalar, we can implement the PolygonTrait methods. We need this distinction between a column of geometries and a single geometry because Arrow is a columnar format while these traits are implemented on singular geometries.

Arrow2 also exposes dynamically typed Scalar objects, like StructScalar, but I decided not to go with those because I didn't want to be doing dynamic dispatch on every row.

I hacked enough together to get geopolars's area method to call the trait-enabled function directly. It compiles and runs, but I have a bug somewhere in the area function and all results are 0.0.

Conclusion

It's a big undertaking, but in the long term, even after all of this, I'm still motivated to get traits working and stable. I welcome harsh feedback though 🙂

[kylebarron]

I read that blog post a few times when it came out and I think need to read it a few more times to understand it 😅. The specific compiler error I was hitting there was

use std::slice::Iter;

impl<'a> LineStringTrait<'a> for LineString<'a> {
    type ItemType = Point<'a>;
    type Iter = Iter<'a, Self::ItemType>;

    fn points(&'a self) -> Self::Iter {
        (0..self.num_points()).into_iter().map(|i| self.point(i).unwrap())
    }
}

error[E0308]: mismatched types
  --> geopolars/geopolars-arrow/src/linestring/scalar.rs:29:9
   |
28 |     fn points(&'a self) -> Self::Iter {
   |                            ---------- expected `std::slice::Iter<'_, point::scalar::Point<'_>>` because of return type
29 |         (0..self.num_points()).into_iter().map(|i| self.point(i).unwrap())
   |         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ expected struct `std::slice::Iter`, found struct `std::iter::Map`
   |
   = note: expected struct `std::slice::Iter<'_, point::scalar::Point<'_>>`
              found struct `std::iter::Map<std::ops::Range<usize>, [closure@geopolars/geopolars-arrow/src/linestring/scalar.rs:29:48: 29:51]>`

and I couldn't figure out how to return an iterator using a closure 😕

[kylebarron]

So after reading the blog post a couple more times, my thought is that the GAT feature might be required in order to have the trait be generic over CoordNum? So that points() can return an iterator containing items generic over CoordNum? But the errors I get now with just f64 coords feel separate from GAT...?

[michaelkirk]

Yeah, it might not be helpful for the problem at hand. Sorry for the distraction!

[kylebarron]

No worries! I figure it's probably worth benchmarking what I have now just to see if accessing geometries from arrow and a vec of georust structs is in the same ballpark

[kylebarron]

I was able to fix this! The main issue is that I was trying to return a closure instead of a concrete type which implements Iterator. Learning from the awesome arrow2 crate, I was implement concrete iterators for each array type, and restore the iterator type in each trait.
https://github.com/geopolars/geoarrow/blob/ff1f840ff4b64a9f1854672f59a2bdc4acd8df8f/src/multipolygon/scalar.rs#L35-L37
https://github.com/geopolars/geoarrow/blob/ff1f840ff4b64a9f1854672f59a2bdc4acd8df8f/src/geo_traits/multipolygon.rs#L24-L26

[michaelkirk]

Thanks for fleshing out some working code to move this along!

It's not clear to me today (especially without trait specialization) if there's a way to define a trait that takes any implementation of another trait as an argument, at least not without resorting to trait objects and dynamic dispatch.

I may be misunderstanding, but would something like a blanket implementation solve the problem? e.g. how we impl Translate for anything that implements AffineOps: https://github.com/georust/geo/blob/main/geo/src/algorithm/translate.rs#L49

(How) do GeometryCollection's fit into this proposal? A heterogenous bag of potentially different geometry types seems like it could be hard to reconcile with a column store, but I'm not really familiar with the subject.

[kylebarron]

like a blanket implementation solve the problem

It's been a month or two since I was exploring the "trait on trait" approach. I got super lost in compiler errors at the time, so for now went with a functional approach so I could get something working. I'll try again and come back with the compiler errors I'm hitting.

(How) do GeometryCollection's fit into this proposal?

I think GeometryCollections do/should fit somewhere in this proposal, but haven't taken the time to prototype them. At a minimal level I could imagine something like

pub trait GeometryCollectionTrait<'a>: Send + Sync {
    type ItemType: // Union of all other geometry traits...?;
    type Iter: Iterator<Item = &'a Self::ItemType>;

    /// The number of geometries in this GeometryCollection
    fn num_geometries(&'a self) -> usize;

    /// Access to a specified interior ring in this Polygon
    /// Will return None if the provided index is out of bounds
    fn geometry(&'a self, i: usize) -> Option<Self::ItemType>;
}

I don't know off the top of my head; I would have to think about how Union of all other geometry traits would work... would it return a trait object that needs to be downcasted?

A heterogenous bag of potentially different geometry types seems like it could be hard to reconcile with a column store, but I'm not really familiar with the subject.

To be clear, my personal use case is columnar memory, but I don't think the traits themselves should have any knowledge of/be designed for a specific memory type. FlatGeobuf is not columnar (and supports GeometryCollection) and I want that to be well supported by these traits as well.

In the GeoArrow-specific case, there's a proposal for how to store a union type, so that there can be different geometry types across each row, but no one has discussed GeometryCollections yet.

[kylebarron]

I may be misunderstanding, but would something like a blanket implementation solve the problem? e.g. how we impl Translate for anything that implements AffineOps: main/geo/src/algorithm/translate.rs#L49

I'm not sure it's possible when you want to have multiple implementations over multiple traits

use geo_types::traits::{LineStringTrait, MultiLineStringTrait};

impl<'a, T, G> HaversineLength<T> for G
where
    T: CoordFloat + FromPrimitive,
    G: LineStringTrait<'a>,
{
    fn haversine_length(&self) -> T {
        todo!()
    }
}

impl<'a, T, G> HaversineLength<T> for G
where
    T: CoordFloat + FromPrimitive,
    G: MultiLineStringTrait<'a>,
{
    fn haversine_length(&self) -> T {
        todo!()
    }
}

gives

error[E0119]: conflicting implementations of trait `algorithm::haversine_length::HaversineLength<_>`
  --> geo/src/algorithm/haversine_length.rs:96:1
   |
86 | impl<'a, T, G> HaversineLength<T> for G
   | --------------------------------------- first implementation here
...
96 | impl<'a, T, G> HaversineLength<T> for G
   | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ conflicting implementation

And if I try to implement for both a concrete type and a blanket implementation, that won't work either

impl<'a, T, G> HaversineLength<T> for G
where
    T: CoordFloat + FromPrimitive,
    G: LineStringTrait<'a>,
{
    fn haversine_length(&self) -> T {
        todo!()
    }
}

impl<T> HaversineLength<T> for Line<T>
where
    T: CoordFloat + FromPrimitive,
{
    fn haversine_length(&self) -> T {
        let (start, end) = self.points();
        start.haversine_distance(&end)
    }
}

error[E0119]: conflicting implementations of trait `algorithm::haversine_length::HaversineLength<_>` for type `geo_types::Line<_>`
  --> geo/src/algorithm/haversine_length.rs:64:1
   |
53 | impl<'a, T, G> HaversineLength<T> for G
   | --------------------------------------- first implementation here
...
64 | impl<T> HaversineLength<T> for Line<T>
   | ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ conflicting implementation for `geo_types::Line<_>`
   |
   = note: upstream crates may add a new impl of trait `geo_types::traits::LineStringTrait<'_>` for type `geo_types::Line<_>` in future versions

[kylebarron]

(How) do GeometryCollection's fit into this proposal? A heterogenous bag of potentially different geometry types seems like it could be hard to reconcile with a column store, but I'm not really familiar with the subject.

I was able to reuse most of the rust-postgis implementation to get geometry collections to compile as a trait and GeometryType enum. If for now we focus on being able to read data through traits (i.e. no writing), then I can punt on how to store geometry collections in columnar memory, and create a fake geometry collection shim that still allows me to implement the full GeometryTrait type on geoarrow

[michaelkirk]

Generic f32/f64

For now, for simplicity, I assume f64 in the traits, but this is something that I would expect to be r

Also, I'm hopeful we could keep some kind of CoordNum based generic number type where possible. Integer types and "bring your own precision model" isn't something I've personally used yet in geo, but I think some people might appreciate the option, and might be kind of a novel thing that we can support better in rust than in, e.g. cpp, so I'm a little resistant to dropping it.

[kylebarron]

Not trying to suggest dropping support 🙂. I just need to look into the geo source more to understand how that would work. Would we need two different sets of geometry trait definitions? One over CoordFloat and another over CoordNum, given that some functions might only accept CoordFloat and not CoordNum? Still being relatively new to rust it's easy to get lost in the generics 🙂

[michaelkirk]

I don't think we'd need different trait definitions, since CoordFloat is a subtype of CoordNum.

I'd guess the geometry traits would have a generic scalar type (probably as an associated type rather than a full blown generic), and that we'd constrain the trait implementations based on this associated type.

For example, I have a super half-baked go at adding a CoordinateTrait here:
1b0b9e6#diff-add86e5b21c0020878489b2e6a7605c9fa7164b53dda3fb1d899be1146812cf7R42

Traits that work with integers are impl'd for CoordNum like Translate here:
1b0b9e6#diff-ab765888e422e7ac563a36df6a35fb0db704290cccfa4b69b9cae0e782c6b2eeR52

Where traits that need floats are impl'd for CoordFloat like Skew here:
1b0b9e6#diff-50e3a37c19483d58ef3fc2055660718a42ec177961684b6817a97d6cba732ad5R125

That branch was kind of a wild spike, so I'm not vouching for the implementation particulars, but I hope it's enough to serve as an answer to your question about different trait definitions.

[kylebarron]

Also, I'm hopeful we could keep some kind of CoordNum based generic number type

In main...kylebarron:geo:kyle/traits-v2 I have it working being generic over CoordNum. It does use a GAT for the base PointTrait and CoordTrait, as we need to use Self::T in i.e. x()

pub trait CoordTrait: Send + Sync {
    type T: CoordNum + Send + Sync;

    /// x component of this coord
    fn x(&self) -> Self::T;
}

[metasim]

In Discord, @kylebarron mentioned concern about "how to adapt georust algorithms to use traits in a backwards-compatible way". I responded there since that concern wasn't (initially at least) brought up here, but it was suggested I include my comment here. It's a minor one, and more of a side-thought, but perhaps worth considering.

metasim:

Your comment about backward-compatibility triggered a thought: In my past public projects (namely, RasterFrames), I've worked tirelessly for backward compatibility, treating it as sacrosanct. Often I would think "How can I make this change while minimizing others' irritation?". In hindsight--while I still believe that should be a consideration--I think I made regretful design compromises because I overemphasized the concerns of existing users (they seemed precious) at the expense of future users (they seemed rare).

Today, In addition to concern for existing users, I now ask: "How can I make this change that best sets the stage for future adoption?". I find it's rarely an either/or decision. Usually there's a continuum of options.

So I'm wondering if you interpreted "backwards-compatibility" as a continuum, are there ways to provide smooth pathways for libraries to adopt trait geometries while allowing you to keep your design goals intact? Something less than compile-time compatibility, but easier than "rewrite all the things".

I'm not sure my idea/question is a good one, and maybe there is a clear way toward backward-compatibility that the more experienced may offer, but I think trait-based geometries has a lot of power (for all the reasons you outlined in your post) and questioning the importance of backwards-compatibility vs. green-field vs. something in-between might lead to more creative and sustainable solutions.

kylebarron:

Mainly it's that I'm pretty new to the georust project (and rust in general) and so I definitely don't want to be coming in saying "let's change everything!" with a lack of respect for why things were architected this way in the first place.

So I'm wondering if you interpreted "backwards-compatibility" as a continuum

I take it as a given that whatever trait definitions are decided on, those traits will be implemented for the existing geo-types structs. So I think that if someone had a function that today accepted geom: geo_types::Polgon, that should be able to be changed to geom: impl PolygonTrait with no actual runtime changes, given that Polygon would implement PolygonTrait.

The main thing I ran into is that it seems easier to do that from a functional approach instead of a trait approach.

So what I could envision is having the trait delegate to public top-level functions that accept traits. Existing users could still use trait methods, while new users could opt-in to the functional approach if they want to use geometry access traits.

[RmStorm]

For everyone ending up here. This effort is alive and well, particularly in #1157.