ENH: prepared geometry as additional (cached) attribute on GeometryObject

[jorisvandenbossche]

@brendan-ward @caspervdw I experimented today with a possible alternative approach on supporting prepared geometries, compared with #84

I added a second static attribute to the GeometryObject struct, which can optionally (for a prepared geom) be populated.

For me, the main reason I wanted to test this approach is that this keeps the geometry and its prepared version together. This means that the user doesn't need to keep track of a separate array with the prepared versions of their data (which could then get out of sync, etc): you just have your single array of geometries, that also can be prepared.

The ufuncs that I added here are just a quick hack to get something working. But from a quick test, it seems to be working:

import geopandas
import pygeos

df = geopandas.read_file(geopandas.datasets.get_path("naturalearth_lowres")) 
countries = pygeos.from_shapely(np.array(df.geometry))  
df2 = geopandas.read_file(geopandas.datasets.get_path("naturalearth_cities"))
cities = pygeos.from_shapely(np.array(df2.geometry)) 

# the current intersects result
In [10]: pygeos.lib.intersects(countries, cities[:, None])    
Out[10]: 
array([[False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       ...,
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False],
       [False, False, False, ..., False, False, False]])

# timing the existing function
In [14]: %timeit pygeos.lib.intersects(countries, cities[:, None])  
18.8 ms ± 91.4 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# timing the prepared version (but without prepared geoms -> just an extra if check)
In [16]: %timeit pygeos.lib.intersects_prepared(countries, cities[:, None]) 
18.9 ms ± 62 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

# preparing the geometries
In [18]: temp = pygeos.lib.prepare(countries) 

# calling the prepared version again, but now with prepared geoms
In [19]: %timeit pygeos.lib.intersects_prepared(countries, cities[:, None]) 
1.27 ms ± 3.16 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

So that shows a POC of having the prepared geometries in the same array, and also shows the nice speed-up that you can get from that!

Some additional thoughts on this:

• It gives some additional memory usage, also when not using prepared geometries (the extra static attribute). If your geometries have a certain size (eg for polygons), this will probably not be very significant, but for points it might be more substantial.

• I now made a separate intersects_prepared version, mainly to be able to compare with the existing intersects. In principle, this can be a single function that can handle both (it seems the extra check doesn't give much slowdown, at least on this test case)

• I added a ufunc to prepare the geometries (based on WIP: support for prepared geometries #84). We maybe also want to "unprepare" in case you want to get rid of some memory usage.
  But in principle, the predicates that support prepared geoms could also "prepare on the fly". That might make it even easier for the user, but gives a bit less control. Most of the time you probably just want that pygeos does it under the hood for you, but there might be cases where prepared geoms are less interesting (I am not too familiar on this, we should probably try to get some good test data of various kinds for such things).

• It modifies the Geometry objects inplace in some way. It doesn't change their meaning, though (it doesn't change coordinates), so I would say this doesn't violate the principle of an immutable geometry object.

[jorisvandenbossche]

I copy pasted some lines of code here to get it working, but in principle we don't need to return anything here (the geoms get modified in place). But I need to check if ufuncs can have zero return values.

(this is the part that can still give segfault after a while, as I suppose the return value here is not properly set up for GC)

[brendan-ward]

@jorisvandenbossche thanks for working on this! Also, very sorry not to have had a chance to make more progress digging into #84

What was the time for preparing the geometries? That might be an interesting data point in thinking through if preparing on the fly makes sense - though likely we'd need to run that on a large collection of polygons to get a good heuristic of the time impact.

Where I would be concerned / surprised is if a function that does preparing on the fly modifies the array I pass into it, as opposed to creating prepared geometries exclusively within the bounds of the predicate function. Were you thinking the former or the latter?

I've struggled with issues of how users should keep track of prepared geoms too, and I like your solution of keeping the prepared & original geometry together. I think in an ideal world, a user wouldn't need to worry about prepared geometries at all; they'd just be used under the hood in the predicate functions where applicable. Meaning: it is possible that we could convert all the predicate functions to prepare geometry on the fly. I liked the developer ergonomics of using prepared geometries on the fly in #87; one less concept to learn and keep track of, but much faster results.

I think the interesting cases to ponder here are the tradeoffs between:

1. cases where you want to prepare the geometry and use them in several predicates (either over the lifetime of the program, or as a chain of different predicates), if the time to prepare the geometries is non-negligible.
2. cases where you know preparing your geometry costs you more time than you gain from faster predicate functions.

It might be that both of those cases are rare enough that preparing on the fly within the predicate functions provides the best general API.

[jorisvandenbossche]

What was the time for preparing the geometries?

I pushed a version of the prepare ufunc with no return value, so that made it possible to time this (before the return value was not properly GCed, so doing that in a time loop segfaulted). And it seems that for this case, the preparing is only taking a tiny bit of time, while being very beneficial for the intersects:

In [7]: %timeit pygeos.lib.intersects(countries, cities[:, None])  
18.8 ms ± 291 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [8]: %timeit pygeos.lib.intersects_prepared(countries, cities[:, None]) 
18.9 ms ± 140 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [9]: %timeit pygeos.lib.prepare(countries) 
26.5 µs ± 101 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [10]: pygeos.lib.prepare(countries)  
Out[10]: ()

In [11]: %timeit pygeos.lib.intersects_prepared(countries, cities[:, None])   
1.23 ms ± 5.34 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

The current version of prepare still has a memory leak though, since I am just overwriting the PreparedGeometry pointers even if they already exist (so meaning they don't get cleaned up by python). Fixing this will make the function even faster when they are already prepared (avoiding you do it a second time), but then it's harder to time ;)

[brendan-ward]

@jorisvandenbossche
I did a bit of testing with this, with some more complex data from one of my projects: 2.5k polygons (waterbodies) and 284k linestrings (rivers / streams):

prepare():

• polygons: 0.0733s
• lines: 0.0026s

If I then do a spatial join to find potentially overlapping features (14.5k potential intersections):
intersects: 4.92s
intersects_prepared: 0.18s

Even one of my worst-case datasets at 2M lines can be prepared in ~ 1s.

The performance gain of using prepared intersections greatly outweighs the time it takes to prepare geometries.

Based on that, I think we should simply prepare on the fly in every predicate function, and also in the STRtree.query() as per #87, and not expose prepared geometry constructs into user-space.

I definitely don't want to negate any of your work here, but it seems like it would simplify the API and reduce number of concepts that users need to understand (i.e., when do I prepare? how do I unprepare?), while making it so that the normal predicate functions are blazing fast (compared to their unprepared counterparts).

Then, if we later see there is a compelling case NOT to prepare on the fly (e.g., prepare takes longer than predicate), perhaps we could expose that as a parameter?

[jorisvandenbossche]

Cool, thanks for testing!

Indeed, it seems the time it takes is small compared to the other typical operations (and the advantage big for those predicates).
Then there is still the other aspect: memory. If we do it on the fly / store them in the GeometryObject, without much control for the user on this, that means some additional memory usage. However, from looking a bit into the GEOS API, I think it keeps the original geometry as a reference (so it doesn't fully copy it, so it wouldn't lead to duplication if we have both in GeometryObject), and the preparation itself is mostly caching the envelope I think in many cases. So I think it shouldn't be a problem.

I definitely don't want to negate any of your work here, but it seems like it would simplify the API

No, no :) The explicit prepare ufunc was mostly copied from your PR, so no problem :)
The main thing I wanted to experiment with here was having the prepared geometry as an additional attribute in the GeometryObject (compared to the separate pyobject), which is needed for the "on the fly" preparation.

I am wondering if we still want some flag to enable this or not, if not for testing (as we are doing here).

---

BTW, are those test data you used open data? It would be good to gather a set of data examples with varying characteristics to test things like this (also for things like sjoin when we are going to add that))

[jorisvandenbossche]

Another thought: how about points, I suppose for those it doesn't matter much to prepare them? But does it harm?

[brendan-ward]

For points, it looks like very little overhead (if I'm reading the GEOS source correctly): https://github.com/libgeos/geos/blob/master/src/geom/prep/PreparedPoint.cpp

I don't think there is any harm there.

The data I have been using are mostly derived from open data sources:

• https://www.usgs.gov/core-science-systems/ngp/national-hydrography/nhdplus-high-resolution
• https://www.census.gov/cgi-bin/geo/shapefiles/index.php
• https://www.epa.gov/eco-research/level-iii-and-iv-ecoregions-continental-united-states

(several steps of preprocessing)

I think it would be OK to share those derivatives for testing purposes (given proper attribution, disclaimers, etc), though they are pretty big! Is there an existing library of datasets used by geopandas / etc that is not in the codebase? Happy to try and pull out useful extracts for testing...

[jorisvandenbossche]

OK, I updated this PR, cleaned up the code, to have something minimal that could be merged (to unblock other PRs to use prepared geometries or to further optimize this).
For now, I kept the prepare ufunc (but fixed its memory leak) and only converted intersects to use the prepared geom if available. But next step could be to always use the prepared version, and thus prepare the geom on the fly if needed.

[jorisvandenbossche]

@pramsey Sorry for pinging you here, but I wanted to ask a quick question for which you might have insight.

Based on some basic timings (using polgyons and points layer, with an intersects operation), using the prepared version of the predicate function gives a very clear speed boost (as you advised to us on twitter! :-)) And the overhead of doing the preparation is almost nihil.

But so the question we are asking ourselves: should we then just always use the prepared versions of the predicates? (and keep the GEOSPreparedGeometry cached with the GEOSGeometry)

Or, based on your experience with this, are you aware of cases where you might not want to use prepared geometries for the predicates? (in which case it could be useful to give some control to the end user on using prepared geometries or not)

[pramsey]

That's a hard call, because it's very pattern based... the overhead isn't really nil, it's around O(n), but once you've done the preparation the computation is around O(log(n)) so doing it the first time you get asked a predicate question is not a bad idea. The trouble is you're then assuming you'll be asked a second predicate question, regarding the same geometry.
In PostGIS world, we do prepare them the second time we are asked a predicate question for the same geometry, as that kind of indicates that we are in a repeating pattern, and will likely get a lot of reuse from the preparation.
But you're right, the different between answer and prepare-and-answer is not as large as one might think.

[jorisvandenbossche]

Thanks for the response!

Maybe a less vague question: for example in PostGIS, if you do a ST_Intersects(..), is this simply always using a prepared intersects? Or do you still use some kind of heuristic based on which you use the "normal" GEOSIntersects_r?

so doing it the first time you get asked a predicate question is not a bad idea. The trouble is you're then assuming you'll be asked a second predicate question, regarding the same geometry.

That's true. On the one hand, I assume that keeping the prepared geometry "cached" doesn't give much overhead, even if you didn't need it (there is some memory overhead though). But on the other hand, since the preparation itself is so fast, there is maybe also not much to be gained in keeping them cached after the first predicate question.

[jorisvandenbossche]

I did some timings with a variant where we always prepare the geom when needed (in the predicate functions) but without caching it.
For the case I was timing above (pygeos.lib.intersects(countries, cities[:, None])), this is actually quite a bit slower (~9ms, vs 1.3ms with caching as in this PR). However, what we do in this code snippet is a full cartesian product, and so for each polygon, we test it multiple times against one of the points. The ufunc broadcasting appraoch gives here a crude appraoch, where the caching clearly helps. But in this snippet we are actually doing a kind of a spatial join, which we would otherwise do with a spatial index or custom implementation re-using the prepared geoms. So I am not fully sure we should optimize for the "ufunc broadcasting" use case.

[brendan-ward]

@jorisvandenbossche thanks for running more timing tests on this! With the caching, are you doing that on the fly, on first instance of intersects for each geometry? Or is there time outside of that loop used for preparing?

Are there other possible use cases for ufunc broadcasting beyond a spatial join type operation that would be useful to consider here?

What is the memory overhead for prepared geometries? If it is more than a tiny amount, or variable based on the complexity of the geometry, that might be something for us to consider in terms of memory vs performance.

What about copying geometries to new arrays? Do the cached prepared geometries go with them, or need to be recreated?

[jorisvandenbossche]

With the caching, are you doing that on the fly, on first instance of intersects for each geometry? Or is there time outside of that loop used for preparing?

Sorry for not being very clear. So the two things I was comparing:

• The approach in this PR: cache the prepared geoms in the GeometryObject (whether they are manually pre-computed and cached, or computed and cached on first usage doesn't matter much for performance I think, but more for the user interface)
• Don't cache the prepared geom, but always create+destroy it when needed. If prepared geoms are then used in the broadcasting intersects ufunc as in the code example, this gives a noticeable overhead. See ec9818d for how the code looks for that test.

With this PR (using the variables of the top-post):

In [2]: %time pygeos.lib.prepare(countries)  
CPU times: user 89 µs, sys: 49 µs, total: 138 µs
Wall time: 143 µs

In [3]: %timeit pygeos.intersects(countries, cities[:, None])  
1.39 ms ± 79.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

With the alternative approach (second bullet point above), always recomputing the prepared geometry when needed:

In [4]: %timeit pygeos.intersects(countries, cities[:, None])  
9.66 ms ± 515 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

and as comparison, with master without prepared geometries at all:

In [2]: %timeit pygeos.intersects(countries, cities[:, None])  
19.4 ms ± 722 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

So for this "brute force" intersects of all combinations of geometries with such a broadcasting ufunc, caching the prepared geoms is clearly beneficial. The only thing I am not sure about is if that is not optimizing the wrong thing.
Because for this use case, ideally you would use a spatial index. Using your PR gives this (it's of course not exactly the same, as the output are indices and not a boolean mask)

In [2]: %%timeit 
   ...: tree = pygeos.STRtree(countries) 
   ...: [tree.query(geom, predicate="intersects") for geom in cities]
943 µs ± 15.3 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

which is even a bit faster than the caching of prepared geoms (and moving this loop over the points into C for an actual sjoin algo will make it even faster).

What about copying geometries to new arrays? Do the cached prepared geometries go with them, or need to be recreated?

That's an advantage of caching the prepared geometries, I think. When copying geometries to a new array (eg taking a subset of an array), those geometries are not actually copied, only new pointers to them are made. So the prepared geoms can be reused and don't need to be recreated. While an STRtree needs to be reconstructed when taking a subset of the array.

I should maybe try to combine our PR on STRtree and this one, to see if caching would be beneficial in such a situation (I assume it will not give much benefit there).

[brendan-ward]

@jorisvandenbossche sorry - I didn't respond to your question in #87. What you outlined in your comment there seems reasonable, and would allow us to get or create cached prepared geometries in that context.

I'd be curious to know how combining these two PRs performs differently. In theory, I'd hope that the performance of a single tree query with predicate (e.g., intersects) performs relatively the same with the existing approach or with the use of cached prepared geometries (when accounting for the time to create them in the comparsion).

Where I'd expect the cached approach to have an edge is repeated queries against the tree with different subsets and predicates, i.e., first query those that intersect, then of those, query those that are completely contained and query those that cross (I do this in another project, but I think maybe that's an edge case).

[brendan-ward]

New thought, likely half-baked. Instead of adding a parameter to predicate functions to enable / disable use of prepared geometries, so that user can disable if they know prepared geometries are a bad solution for their specific dataset, what about adding an attribute that could be set on the geometry to disable preparing / caching?

Not sure how users would set that, but it seems like something that needs to be handled at the geometry object level. Otherwise various operations that touch a geometry and request the prepared version all need to be called with the same parameter to disable use of prepared.

That could also give better granularity: you could disable prepared ops on some but not all of your geometries in an array.

[caspervdw]

I'd like to revive this discussion. After not thinking a long time about it, I now arrive at that this PR already has the best strategy:

• allocate an extra static attribute _ptr_prepared on the Geometry object
• the prepare function fills it
• all functions that can make use of it let them check if there is a prepared geometry available. They should not prepare geometries themselves.
• there should also be a destroy_prepared_geometries function that frees some memory.

Note that these prepare and destroy_prepared_geometries will change the Geometry inplace! I think this is OK as the geometry does not change.

[brendan-ward]

all functions that can make use of it let them check if there is a prepared geometry available. They should not prepare geometries themselves.

How do we want to approach this with STRtree with predicates? In the context of querying with predicates, we only really need the geometries whose bounding boxes have overlaps with tree geometries to be prepared in this case, not all input geometries to the query (esp query_bulk) operation. Calling prepare on all inputs could be overkill in this case. But, we likely get a major benefit from using prepared geometries in the tree - to the point where we could call get_geom_prepared to retrieve the cached prepared geometry if it exists, but otherwise create it ourselves in STRtree if missing, so that all predicate operations in the tree use prepared geometries.

prepare and destroy_prepared_geometries will change the Geometry inplace! I think this is OK as the geometry does not change.

Meaning that the GeometryObject gets modified by setting / unsetting a pointer only, and the original GEOSGeometry does not get touched at all, right?

Also:
prepared geometries may already exist when you call prepare. Presumably we can just check the pointer and skip those that already have prepared geometries, since there aren't operations that would mutate the geometry separately from the prepared geometry (set_coordinates would create new geometries with no prepared geometries).

[caspervdw]

How do we want to approach this with STRtree with predicates? In the context of querying with predicates, we only really need the geometries whose bounding boxes have overlaps with tree geometries to be prepared in this case, not all input geometries to the query (esp query_bulk) operation. Calling prepare on all inputs could be overkill in this case. But, we likely get a major benefit from using prepared geometries in the tree - to the point where we could call get_geom_prepared to retrieve the cached prepared geometry if it exists, but otherwise create it ourselves in STRtree if missing, so that all predicate operations in the tree use prepared geometries.

We could let the query_bulk operations prepare geometries where it thinks it is necessary, and even clean them up afterwards. This could be controlled by extra argumenst query_bulk(..., prepare=True, destroy_prepared=False). As this function is more high-level than the functions in pygeos.predicates, I am 👍 of sticking in these features.
Note that the STRTree edits should be out of scope of this PR, but it is good to consider this already.

Meaning that the GeometryObject gets modified by setting / unsetting a pointer only, and the original GEOSGeometry does not get touched at all, right?
Also:
prepared geometries may already exist when you call prepare. Presumably we can just check the pointer and skip those that already have prepared geometries, since there aren't operations that would mutate the geometry separately from the prepared geometry (set_coordinates would create new geometries with no prepared geometries).

Indeed, we should skip preparing geometries if they exist already. It never makes sense to "reprepare" as geometries are immutable.

Another feature of this 'in place' preparing is that we can selectively prepare geometries based on some crude estimate:

pygeos.prepare(arr[pygeos.area(arr) > 10])
pygeos.prepare(arr[pygeos.get_type_id(arr) == 4])

[caspervdw]

@jorisvandenbossche I resumed the work on this PR. We can iron out the strtree implementation later, but I gathered from above conversations that the stuff in this PR is something we agree on already.

Before going through the full list of predicates, this should be reviewed.

[brendan-ward]

Thanks for getting this rolling again @caspervdw

Overall I think we were agreed about the overall behavior here, and evaluating this for strtree should be a different PR / discussion.

What do you think about adding a function is_prepared? This would return True / False based on whether or not the _ptr_prepared is NULL, and would avoid having to use that attribute during tests or elsewhere in the python API where we want to know if it is prepared. (except in tests of is_prepared where we need to test against that attribute to know it is working properly).

[brendan-ward]

It would be good to add a note here that this will return NULL if the geometry has not yet been prepared.

Also - get_geom has extra processing to handle cases where obj is None. It seems like it is important to do the same here, right?

[caspervdw]

Instead, I implemented get_geom_with_prepared. This because get_geom_prepared would always go after get_geom, better to combine them in one function (to be sure the checks are done).

[brendan-ward]

Should there be more guidance here or in docs about what those specific circumstances are?

[caspervdw]

Thanks for the review @brendan-ward . I like the idea of an is_prepared function, but to vectorize it properly I like to wait for #207 to go in (the O_b type function is exactly what I need for this)

OK like this to start implementing the other prepared predicates?

[caspervdw]

@jorisvandenbossche @brendan-ward Could you take a look at this PR?

[brendan-ward]

@caspervdw thanks for the updates here.

A few suggested changes mostly around docstrings and params.

Added #235 to leverage this for stree.query_bulk.

[brendan-ward]

Thanks @caspervdw - looking forward to having this available!

[jorisvandenbossche]

@caspervdw BTW, thank a lot for reviving and finishing this up!