Support for spherical geometry (S2) in Python

[jorisvandenbossche]

General issue to track the idea of improving the support for working with spherical geometries (instead of planar) in GeoPandas and more generally the Python geospatial ecosystem.

Some background links:

• The S2 library from google: https://s2geometry.io/
• R bindings for this: https://r-spatial.github.io/s2/
• Announcement from R's sf that they will use S2 by default: https://r-spatial.org/r/2020/06/17/s2.html#sf-10-goodbye-flat-earth-welcome-s2-spherical-geometry
• A geopandas issue about this: ENH: Support for spherical geometry geopandas#2098

[jorisvandenbossche]

And the previous community meeting, we had a discussion about this with @benbovy. Pasting the notes from that meeting here:

• see this thread for context
• https://github.com/benbovy/pys2index
• overview from Benoit:
  • lightweight wrapper for S2 point index with API similar to scipy.spatial.cKDTree
  • performance in benchmarks so far is quite fast
  • would like to have vectorized wrappers for S2
  • S2 appears to be actively maintained and about to get additional functionality soon
  • used python-xtensor and pybind11 to work with S2 and numpy arrays
• Two possible approaches to integrate:
  • own way to store geometries specfic to backend engine; convert geometries on the fly to S2 objects as part of specific operation (e.g., predicate)
    • R library converts on the fly to GEOS or S2 as needed
    • More info: https://r-spatial.github.io/sf/articles/sf7.html
  • Wrapper classes:
    • pygeos uses Python C extension wrapper for GEOS geometries so that GEOS objects are managed according to Python object lifecycles
• consider putting in a request for NumFOCUS small development grant to start building out some of this support (next cycle may open early 2022): https://numfocus.org/programs/small-development-grants

[benbovy]

Some updates on this:

I've tried using pybind11 and xtensor-python to implement some vectorized functions working on numpy object arrays where array elements are Python-wrapped s2 objects (i.e., an approach similar to the one used in pygeos). Without success so far, unfortunately (see benbovy/pys2index#4).

There has been some effort to support numpy object arrays and/or ufuncs in pybind11 (e.g., pybind/pybind11#1152 and also this fork). I think it would really help, although it is not yet clear to me if/how we can adapt it in this case. I left a comment there: pybind/pybind11#1152 (comment).

[brendan-ward]

Thanks for the update @benbovy !

I haven't been able to fully get my head around the use of pybind11 / xtensor-python parts here, but a couple thoughts came to mind:

1. a common use would be in GeoPandas, where a geometry array is 1-dimensional. So supporting n-dimensional ufuncs is maybe a nice-to-have rather than must have? Thus the essential bit is really just a loop over a 1D array of objects (pointers to wrapped S2 entities). It seems like we could leverage numpy to do the broadcasting before calling into the pybind11 / C++ tier? (not sure, haven't really thought this out)

2. maybe it would be easier to leverage numpy C API (i.e., meet it where it's at), as is done in pygeos, and then build this as a Python C extension a bit more manually (like pygeos) rather than via pybind11? To do this, I think we'd have to create a C API around S2, so in addition to the boilerplate required for Python C extension there would also be some maintenance burden on wrapping S2 as a C API. One advantage of this approach is that it should enable numpy ufuncs in the C tier, much like pygeos. This would need some investigation into how S2 might be best wrapped as a C API, but we could theoretically leverage some of the recent patterns worked out in GEOS (passing coordinate buffers instead of setting individually). Plus there might be some advantage of having similar patterns at the lower (C) level between the two implementations on top of underlying geometry engines (GEOS and S2); this might be able to better leverage contributors between the pygeos / shapely side and S2 bindings side. (though might be trying too hard to make this analogous to pygeos)

[martinfleis]

@benbovy on a less technical note, if you still have a capacity and an interest in applying for the NumFOCUS SDG, the timeline for the next round is following:

• Call for Proposals Announcement: February 4, 2022
• Proposal Submission Deadline: March 4, 2022
• Committee Selection Deadline: March 18, 2022
• Notification to Applicants Deadline: April 15, 2022

If we should go ahead, I suggest organising a call in the first half of February to discuss the details.

[benbovy]

Thanks for your feedback @brendan-ward !

A common use would be in GeoPandas, where a geometry array is 1-dimensional. So supporting n-dimensional ufuncs is maybe a nice-to-have rather than must have? Thus the essential bit is really just a loop over a 1D array of objects (pointers to wrapped S2 entities)

I agree, we could certainly start with 1-dimensional arrays. That said, longer term I'd be interested in exploring the possibility of reusing n-dimensional vectorized S2 Python wrappers (as well as pygeos) with Xarray, which is undergoing some heavy refactoring to allow custom indexes. Xarray + geospatial indexes would be great for processing outputs of climate / ocean / landscape evolution / earth system models.

I think we'd have to create a C API around S2

That's my main concern. Unfortunately, there's no s2_c.h equivalent of geos_c.h as pointed out by @paleolimbot (https://youtu.be/zqRhF2XM1CE?t=1197). See also r-spatial/s2#160. Long term it would be great to have a common S2 C API that can be maintained by both the Python and R (and maybe other) geospatial communities, but this represents quite a big task IMHO. Also, it's best if we can avoid vendoring S2 in a C/Python wrapper.

On the other hand, if we can make this possible without much effort, simply using pybind11 (and xtensor-python) to create the Python wrappers would allow us to develop new features / iterate very quickly and would avoid the burden of maintaining an intermediate C API + dealing with the low-level Numpy C API (+ vendoring S2). It not yet clear to me what is needed to make this approach work, though. It would be great to have thoughts from folks at QuantStack. @wolfv, have you ever thought about adding vectorized functions / ufuncs working with numpy arrays of wrapped CGAL objects (as np.object dtype) in scikit-geometry?

[benbovy]

@martinfleis that sounds good. I'm happy to discuss the details within the next two weeks.

Before going further with a NumFOCUS application, I'd like to clarify the technical challenges of the C API vs. pybind11 approaches. TBH I'm not sure to have the full capacity of going deep into either the design of a S2 C API or heavy customization of pybind11 internals, especially within the context of a NumFOCUS small development grant.

[brendan-ward]

@paleolimbot I'd be curious about your thoughts on a shared C API against S2 that could be leveraged for both Python, sf, and anything else that can use a C API. It looks like some of the code already in sf could be leveraged to get this started (per your thinking in r-spatial/s2#160).

Specifically, if this is something of interest to the broader community, should we try and take this up directly with S2? Long term, it seems ideal that the C API lives inside that project, as it does with GEOS, and given the approach with GEOS it seems also reasonable that it not be expected to have full coverage of the underlying C++ API, at least not initially.

However, that might lead to other challenges:

• slower iteration speed in downstream libraries (vs vendor + wrap in those libraries) since it means waiting on C API releases
• since our target audience here is also using GEOS for planar spatial operations, it might be ideal to make the C API around S2 a bit more similar to GEOS (if such is possible), in order to better simplify the calling patterns from Python / R / etc. As opposed to going with a C API that follows the S2 C++ API more closely.

Since it looks like sf has already taken the vendor and partially wrap with C API approach, we'd be curious about the maintenance effort there as we think through our approach here.

[jorisvandenbossche]

From watching the FOSS4G talk of @paleolimbot mentioned above (https://www.youtube.com/watch?v=zqRhF2XM1CE), it is my understanding that the S2 Geometry library itself it quite low level, and all the functionality to work on "simple features" (as what we get from GEOS for planar features) actually lives in the R s2 package.

If that's correct, we shouldn't put effort in directly interfacing with S2, but rather with the C or C++ layer in R's s2. Which would then indeed require making this part of s2 more standalone as outlined in r-spatial/s2#160 / mentioned by Brendan in the comment above.

[paleolimbot]

Thanks all for looping me in here! @jorisvandenbossche is right that anything resembling a GEOS-like interface lives in the R package. More specifically, I copied the BigQuery Geography API (which copies the PostGIS API more or less), because I could test it as I went.

I'm conflicted about a C API for S2 because it robs S2 of some of the things that make it so great (despite the fact that I opened that ticket for myself). That said, S2 C++ is complicated to use and an awesome library that nobody has access to in Python is not doing anybody any good. I have some ideas about how to organize a C API and will give it a try since it seems like there may be interest in helping if I get it started. The nice part about doing this within the S2 R package is that there's near 100% test coverage, although eventually a C API should live elsewhere.

Another approach I've been thinking about is one that operates on vectors of geometries (as ABI-stable C ArrowArrays) rather than individual geometries. Operating on chunks of geometries is nice because you don't have to sacrifice performance on points to support multipolygons (et al.)...you just send more features in a chunk and the implementation can loop without a bunch of tiny allocs or virtual method calls. This depends on the somewhat experimental stuff that Joris, I, and others have been working on with how geometry can be passed around as a struct ArrowArray.

Or a combination of both! I'll start my own C API ticket and see what happens.

I probably missed a few things from the thread...let me know if I didn't catch something important!

[brendan-ward]

one that operates on vectors of geometries (as ABI-stable C ArrowArrays)

This sounds very interesting! I've been wondering where some of the boundaries might be with those (mostly a separate topic), namely serialization to / from file storage vs on the compute side between a client library (e.g., something in Python or R) and the underlying geometry implementation (e.g., GEOS or S2).

Is the idea that the core object that is used in the client library becomes the ArrowArray, which is then passed to a layer in C/C++ that translates from those to GEOS or S2 objects? Is that limited to geometry construction / serialization from the underlying library (GEOS / S2), and otherwise the client holds on to wrapped pointers to GEOS or S2 objects?

you just send more features in a chunk and the implementation can loop without a bunch of tiny allocs or virtual method calls

I was thinking about this too in terms of trying to figure out which operations in a hypothetical C API for S2 would need to work on scalar geometries vs 1D arrays of geometries and have the underlying implementation do the looping rather than the caller (passing arrays seems very appealing). The part where I got hung up was that on the Python side (in pygeos/shapely), we're leveraging looping mechanics from numpy to give us the inner loop over multi-dimensional arrays and then operating on scalar geometries. I.e., it is the client that knows how to do the looping rather than the underlying implementation. Though our most common use cases so far are likely 1D arrays of geometries (though @benbovy 's comment re: n-dimensional arrays is well taken), and I might not be thinking about this correctly at all.

I have some ideas about how to organize a C API and will give it a try since it seems like there may be interest in helping if I get it started

Yes - interested! (may not have tons of capacity to throw at this, but interested in being involved). I was thinking of drafting a Google Doc or similar to outline some thoughts on the C API, but your thoughts on this are going to be a much better starting point. Is there someplace where we (community) can start collaborating around details about the API more specifically?

Is there interest in a meeting to discuss this topic, and shared needs in the S2 / ArrowArray space for both Python and R?

[benbovy]

From this discussion it seems that the best path forward would be to develop a common, GEOS-like C API for S2 that can be leveraged from different "client" ecosystems such as R or Python geospatial libraries. Although you guys certainly have much more experience than I have with GEOS, Arrow and/or the Numpy C API, I'd still be happy to help somehow if I'm able to do that, as I would really like to see more S2 capabilities accessible from within Python.

Some naive questions:

@paleolimbot how much effort do you think it would require to refactor bits from the S2 R package into a target-language agnostic C API? Do you think that short or mid-term we could already have something (even with a few features only) so that we can start building Python wrappers on it?

Perhaps a very naive question: would it be possible to separate the building of a GEOS-like C API from the building of the S2 library itself? @brendan-ward's comment re: ideally the C API would live inside the S2 project, although I don't know if there's interest in that from the S2 maintainers point of view and this would also probably further slower the iteration speed. Note: there are conda packages for the S2 library.

Could we easily reuse with Numpy arrays a C API that is working with C Arrow arrays? I.e., perform broadcasting, convert to Arrow array, call the C API and then convert back to Numpy array at little overhead cost? As a user, I indeed find the pygeos/shapely approach based on Numpy ufuncs quite appealing :), even though I acknowledge that it is probably overkill for 80% of use cases.

I'm conflicted about a C API for S2 because it robs S2 of some of the things that make it so great

I was also wondering about that, but I guess that S2-specific wrappers may live alongside (and more-or-less independently from) a GEOS-like "standardized" S2 C API, depending on what we want to achieve using Python, R,...?

[brendan-ward]

@benbovy (speaking for myself here) - I don't think we've mady any decisions about which direction to pursue or even recommend, we're still gathering info to weigh pros / cons. So please don't take the suggestions around a C API for S2 as trying to shut down your work on the pybind11 / xtensor-python approach!

It does seem that in elaborating how we'd want a C API to work would also inform how we'd want an API via pybind11/xtensor-python to work (though a good bit of that would be inside the C++ tier rather than API exposed to python).

Question: if for the short term we wanted to stick to the 1D approach to get this implemented most quickly, and manage all the looping ourselves as a simple loop over that array, would that sidestep the issues around how to handle object arrays in xtensor-python that you were bumping into? I.e., sidestep the ufunc related issues specifically. Or is the issue more at a fundamental level of trying to have an object array of S2 objects at any dimension? (sorry - not entirely following your box example linked above)

Yet another idea, no idea if it would work. I wonder if there is middle ground of using pybind11 but with a more C style interface rather than object-oriented interface (i.e., don't use pybind to expose C++ classes with methods).

So in pygeos our core geometry object is just a struct that holds a couple pointers and is defined as a type:

typedef struct {
  PyObject_HEAD void* ptr;
  void* ptr_prepared;
} GeometryObject;

(other mechanics around managing the python type wrapping around this omitted here)

Which are then always passed around as PyObject pointers for use in numpy arrays or scalars; numpy doesn't need to know that . The ptr just points to a GEOSGeometry, which (if I follow) is itself just a struct around a GEOS C++ geometry object.

So (if I follow correctly), it seems like it would be possible (in theory, not up on pybind11 yet) to do something similar: constructors return a pointer to a struct (let's call this PyS2Geometry for example) that is what is stored into numpy arrays and passed back to Python; this has has an underlying C++ class underneath (let's call this S2Geometry) that is not visible to pybind11 specifically, but is accessible in your C++ tier by following the pointer. Note that the S2Geometry could just be a shell that manages more specific types underneath or via inheritance.

Then I think you could expose functions like intersects(PyS2Geometry a, PyS2Geometry b) (scalar example) to Python, and in your C++ tier get back the underlying S2Geometry, and call methods or functions using those.

So in this purely hypothetical example, pybind11 is really just helping with the boilerplate around creating a Python C extension, but allowing you to more easily have a C++ implementation underneath rather than having to define a C API in the middle. No idea what that means for working with numpy at the C++ level though.

[benbovy]

No worries at all! I'm really open about which direction to pursue and I'm enjoying the discussion here! Sorry if my tone suggested otherwise.

The box example is also confusing, sorry. I was messing around with pybind11 (and dumped some messy code in a branch), trying to see if we can create vectorized functions with custom types, both exposed to Python. I defined a box class wrapping the S2 equivalent but it could have been any other custom type, not even related to S2.

Note that the S2Geometry could just be a shell that manages more specific types underneath or via inheritance.

That is what I had in mind too (instead of exposing more specific S2 derived types to Python).

pybind11 is really just helping with the boilerplate around creating a Python C extension, but allowing you to more easily have a C++ implementation underneath rather than having to define a C API in the middle.

Yes this is what I wanted to check. Taking your example, I was wondering if pybind11 could even help with the boilerplate around the PyS2Geometry struct and the vectorization, i.e.,

• expose the S2Geometry type in Python: py::class_<S2Geometry>(m, "S2Geometry")
• register it as Numpy object dtype, e.g., using this macro (not available yet in pybind11): PYBIND11_NUMPY_OBJECT_DTYPE(S2Geometry)
• expose intersects(const S2Geometry& a, const S2Geometry& b) as a vectorized Numpy (universal) function: m.def("intersects", py::vectorize(interstects)).
• or do something like m.def("intersects", [](py::array_t<S2Geometry> a, py::array_t<S2Geometry> b){ ... })

However, I couldn't get it to work with the box example above. It is likely that I'm thinking about it too naively and that I'm missing something important. For example, py::vectorize doesn't work with non-POD types. It is not clear either to me how we can separate Numpy stuff from Eigen stuff implemented in this PR and in this pybind11 fork.

[brendan-ward]

@benbovy I did some testing with pybind11 around returning numpy arrays of custom objects here.

Let me preface this by saying I am not at all up to speed on modern C++, esp. pointer semantics, still quite unfamiliar with pybind11, and basically don't know what. I'm doing here. I think this is working for a hypothetical example of wrapping a foreign object into a class that we then return as a pointer. I only tried the 1D array case and did not use xtensor-python.

The things that seemed to make this work were using PYBIND11_NUMPY_OBJECT_DTYPE from the pybind #1152 issue to register our custom type as an object dtype for numpy, and then returning a unique_ptr (should it be shared_ptr??) to our custom class as each element in the array (so that instances weren't copied when assigning into the array for return). Also had to register the custom class with pybind11, so that it knew how to __repr__ it, for example.

For the example I construct a custom Point object from x and y input arrays (looping ourselves), where a Point holds an instance of a foreign object, meant to simulate an S2 object, which allocates and deallocates memory (I haven't yet profiled for memory leaks); it seems like it is calling the destructors properly and should release it.

Point was just meant to stand in for a wrapper around S2 objects, not a specific concrete Point type. Not sure if this solves the specific issues you were bumping into above.