Support pickling

[caspervdw]

Closes #46

After a failed approach with __getstate__ and __setstate__, I chose for the __reduce__ route. __setstate__ does not play nice with the fact that our geometry objects are immutable. Further reading: https://docs.python.org/3/library/pickle.html#pickling-class-instances

__reduce__ should return a tuple (callable, args). A convenient way for us would be the following:

• (pygeos.lib.Geometry, ("<WKT>", ))

However, WKT serialization is less efficient than WKB serialization. We could opt for:

• (pygeos.lib.from_wkb, (b"<WKB>", ))

When testing these approaches, a noticed that both of them have drawbacks.

• WKT does not support SRID serialization and it has an issue with empty points in multipoints ([Done] Raise ValueError on to_wkt of multipoint with empty point #171).
• WKB does not support linearrings

I am still in favour of going the WKB route. Do we accept that linearrings gets transformed into linestrings when pickled-unpickled?

>>> import pickle, pygeos
>>> g = pygeos.linearrings([[(0, 0), (0, 1), (1, 1), (0, 0)]])
>>> pickle.loads(pickle.dumps(g))
array([<pygeos.Geometry LINESTRING (0 0, 0 1, 1 1, 0 0)>], dtype=object)

[brendan-ward]

WKB seems reasonable, given that other similar efforts (e.g., feather / parquet support in geopandas) use WKB as the serializable representation.

linearrings only really make sense in the context of polygons, right? Thus I think of linearrings as an input toward making polygons or extracted from polygons and used as a linestring, but not something that should be serialized as a discrete entity. I think noting this as a known limitation (imposed by WKB spec) is sufficient, and perhaps a one-liner for getting back to a linearring if absolutely necessary?

ring = pg.get_exterior_ring(pg.Geometry("POLYGON ((1 0, 1 1, 0 1, 0 0, 1 0))"))
pg.linearrings(pg.get_coordinates(pickle.loads(pickle.dumps(ring))))

[brendan-ward]

At first glance, it looks like this adds a sizeable amount of complexity with some potential overlaps with the from / to WKB functions in ufuncs. Would there be a benefit to trying to refactor the common bits into functions that are shared between this and the ufuncs?

It looks like the C method gains a significant time advantage when there are multiple geometries:

>>> geom = np.repeat(pg.linestrings(np.arange(1000), np.arange(1000)), 100)
>>> timeit('pickle.loads(pickle.dumps(geom))', number=100, globals={'geom': geom, 'pg': pg, 'pickle': pickle})
0.020391266999922664
>>> timeit('pg.from_wkb(pickle.loads(pickle.dumps(pg.to_wkb(geom))))', number=100, globals={'geom': geom, 'pg': pg, 'pickle': pickle})
1.2352546539998457

[caspervdw]

Thanks for the feedback.

About the linearrings: let's indeed warn the user about this drawback of pickling.

The code duplication bothers me as well. Since #182 , we could implement the Geometry features in pure python, this would enable us to remove a lot of complexity. But this should be accompanied by some benchmarks.

The performance benchmark you sent is very surprising to me. I expected that the vectorization in pygeos.to_wkb would give you a performance benifit. But I don't know how the pickling of an object-dtype ndarray works precisely.

[jorisvandenbossche]

Agreed that WKB should do, the fact that linearrings are not preserved on roundtrip can be a documented limitation.

I am also quite surprised to see that performance characteristic, as I actually did that (converting to WKB vectorized, and then pickling the result of that) in geopandas (https://github.com/geopandas/geopandas/blob/55b99c66081b0281c3ad3ee3ea0fcf14502ba473/geopandas/array.py#L390-L397), assuming this would be faster than letting python pickle each object of the numpy object dtype array individually ...

Now, trying out some other cases, it seems to depend on the lenth of the array.

Original example with a small array (100 elements) of relatively large linestrings:

In [1]: import pygeos

In [2]: import pickle

In [3]: arr = np.repeat(pygeos.linestrings(np.arange(1000), np.arange(1000)), 100)

In [4]: arr.shape
Out[4]: (100,)

In [5]: %timeit pickle.loads(pickle.dumps(arr))
148 µs ± 8.52 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [6]: %timeit pygeos.from_wkb(pickle.loads(pickle.dumps(pygeos.to_wkb(arr))))
13.9 ms ± 1.37 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

Large array of tiny linestrings:

In [7]: arr2 = pygeos.linestrings(np.random.randn(100000, 3, 2))

In [8]: %timeit pickle.loads(pickle.dumps(arr2))
275 ms ± 6.31 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

In [9]: %timeit pygeos.from_wkb(pickle.loads(pickle.dumps(pygeos.to_wkb(arr2))))
157 ms ± 4.74 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

and now the other method is actually faster.

Now, for a comparison relevant to this PR, we should actually compare against having it in Python vs in C, so that's still something else

[jorisvandenbossche]

Ah, I think I understand what is happening. The array created with np.repeat actually contains all the same objects (since it is object dtype array, I suppose it can point to the same python object), and I suppose that pickling has some optimization to detect identical objects being pickled.
Anyway, when creating the same array but with unique objects (the more typical case in practice), the first pickling method is not significantly faster anymore:

In [12]: arr = np.array([pygeos.linestrings(np.arange(1000), np.arange(1000)) for _ in range(100)], dtype=object)

In [13]: %timeit pickle.loads(pickle.dumps(arr))
12.3 ms ± 213 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [14]: %timeit pygeos.from_wkb(pickle.loads(pickle.dumps(pygeos.to_wkb(arr))))
13.1 ms ± 1.26 ms per loop (mean ± std. dev. of 7 runs, 100 loops each)

(and I suppose the vectorization benefit will mainly show for larger arrays than 100 elements)

[jorisvandenbossche]

Comparing with my shapely branch that subclasses from Geometry in Python, and adding a simple __reduce__ there:

    def __reduce__(self):
        return (pygeos.from_wkb, (self.wkb, ))

and repeating the same examples as above:

In [9]: arr = np.array([pygeos.linestrings(np.arange(1000), np.arange(1000)) for _ in range(100)], dtype=object)

In [10]: %timeit pickle.loads(pickle.dumps(arr))
15.9 ms ± 747 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

In [11]: arr2 = pygeos.linestrings(np.random.randn(100000, 3, 2))

In [12]: %timeit pickle.loads(pickle.dumps(arr2))
1.4 s ± 77.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

So for the second case (array with many elements), there seems to be quite some overhead from this being implemented in python (275ms -> 1.4s).

Now, in practice in applications like GeoPandas, we can override this by ensuring we have a custom pickle method (using the vectorized to_wkb/from_wkb method) and so for GeoPandas such a slowdown wouldn't have any impact. But when working with raw numpy arrays of geometries, then it's clearly beneficial to have it in C

[caspervdw]

Thanks for all these benchmarks! And I am happy to read that this PR is already optimal.

We get the extra performance at the expense of some code duplication.

Is this PR good to go?

[jorisvandenbossche]

Working on the shapely subclasses branch (shapely/shapely#983) this morning, I was thinking about the linestring/linearring a bit more: in principle we could "relatively easy" pickle the geometry type id as well, together with the WKB? In which case we could ensure type stability?

Now, we would then need to have a custom "deserialize" function, which is not just the class constructor, so not sure how complex this part would be.

[brendan-ward]

pickle the geometry type id as well, together with the WKB

That feels like it is duplicating the geometry type already stored into the WKB for all those types already well-typed there.

Another idea would be to come up with a relatively safe 32 bit identifier for linear rings (as in, unlikely to collide with future geometry type ids) and store that into the WKB instead of the ID for linear string in bytes 2-5 (first serialize the WKB then update these bytes). Then on decoding, just decode those bytes back to an int, check against this well known code to set a flag that it is a linear ring instead of line string, rewrite the WKB (to flip back to linear string id) and decode as now, and convert to a linear ring afterward.

Since the pickled WKB is not intended to be portable outside this context, perhaps that would be a reasonable compromise? It seems like this would be relatively easy to implement since it mostly involves modifying bytes in a well known location... (not sure about actual complexity or overhead)

[caspervdw]

So we're back on supporting linearrings? We also have the GeometryCollections, these may contain LinearRings as well (recursively).

I am no big fan on going down this rabbit hole again.

[caspervdw]

Would serializing linearrings to WKT (instead of WKB) be acceptable?

[brendan-ward]

It is unfortunate that neither serialization spec (WKT or WKB) perfectly maps to the geometry structures we're able to implement here, and that they are not even consistent with each other. It seems like the decision not to acknowledge linear rings as a valid top-level geometry type was made in the WKB spec, and not really for us to override here. My suggestions above add complexity with unclear benefit.

Personally, I like the route of WKB with linear rings => linestrings as a documented issue. There are mutiple user-side workarounds (convert back to linear ring after unpickling, make into a polygon before pickling) that we can suggest. We can either warn on pickling linear rings, or raise an exception and force users to convert before pickling.

What are we thinking as a version stability policy for pickled objects? Perhaps for now, we warn that we make no stability promise, and that objects can only be unpickled with exactly the same version of pygeos as they were pickled?

[brendan-ward]

Sorry for not reviewing sooner! Thank you for your work on this. We can revisit consolidating code in #165 at a later time.

In addition to the work around pickling, supporting WKB in GeometryObject_new may open up future opportunities (not sure what yet).

A few suggested changes here, the biggest being around tests:

• we should expect an output of linearstrings from linear rings
• it would be good to add tests for 3D geometries and geometries with SRIDs, or create an issue to come back to those later.

Otherwise things looked pretty solid to me. Note: I'm not familiar with pickling, so I may well have missed something important here.

[jorisvandenbossche]

So we're back on supporting linearrings?

Not necessarily, we can still decide not to. The main reason I raised it again is because I noticed in my shapely branch that in current shapely LinearRing is preserved in pickle roundtrip. So it would be a regression for shapely to no longer do that. But so maybe something to discuss with the other shapely devs if we are fine with such a regression.

We also have the GeometryCollections, these may contain LinearRings as well

Indeed .. and here shapely is also returning a LineString after roundtrip:

In [1]: from shapely.geometry import LinearRing, GeometryCollection  

In [4]: g = GeometryCollection([LinearRing([(0,0), (1,1), (0,1), (0,0)])])   

In [5]: print(g)  
GEOMETRYCOLLECTION (LINEARRING (0 0, 1 1, 0 1, 0 0))

In [6]: import pickle 

In [8]: print(pickle.loads(pickle.dumps(g))) 
GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1, 0 1, 0 0))

[jorisvandenbossche]

But so maybe something to discuss with the other shapely devs if we are fine with such a regression.

I opened shapely/shapely#988 to discuss it there as well.

Now, while writing it, I realized we could perfectly override the default pickling behaviour in Python specifically for the LinearRing subclass. So that way we keep it fast / simple in C for all geometry types except LinearRings, but since that is a special case anyway, it's probably also not a problem that pickling many linearrings is a bit slower.

Now, for eg GeoPandas, that would still give an inconsistency, if it uses the vectorized to_wkb on the array first before pickling .. (but that is then something for GeoPandas to decide)

[jorisvandenbossche]

In addition to the work around pickling, supporting WKB in GeometryObject_new may open up future opportunities (not sure what yet).

Personally, I am not a big fan of overloading the Geometry(..) constructor like that, certainly if we don't have a specific use case in mind. I personally find using from_wkb much more explicit.

[mwtoews]

Another limitation of WKB to note here is that it cannot differentiate between POINT(nan nan) and POINT EMPTY (#196 (comment)).

[caspervdw]

Thanks for the review @brendan-ward . I implemented your suggestions. Also I wrote a documentation on Geometry objects, it was missing. Good to merge now?

[brendan-ward]

Thanks for updating this @caspervdw .

The docs are a really nice addtion!

Please see a few outstanding prior suggestions in the PR (comments for a few things in C, initializing a couple variables to NULL); I've resolved all of those that you addressed in your most recent changes.

[brendan-ward]

Suggested change

	be changed inplace. Every PyGEOS operation will result in a new object being returned.
	be changed in place. Every PyGEOS operation will result in a new object being returned.

[caspervdw]

The docs are a really nice addtion!

Thanks!

Please see a few outstanding prior suggestions in the PR (comments for a few things in C, initializing a couple variables to NULL); I've resolved all of those that you addressed in your most recent changes.

I missed them because github folded them for me. I think I addressed them all now.

[brendan-ward]

Thanks @caspervdw !

[jorisvandenbossche]

Can you add a test for hashing of empty point (then this closes also #196, I suppose)

[jorisvandenbossche]

Why not for empty points? (or do you mean it might not preserve the dimension?)

[jorisvandenbossche]

What is the added value of supporting WKB here? IMO that only duplicates from_wkb, while being less explicit.