SpatialPandas design and features

[jonmmease]

spatialpandas

Pandas and Dask extensions for vectorized spatial and geometric operations.

This proposal is a plan towards extracting the functionality of the spatial/geometric utilities developed in Datashader into this separate, general-purpose library. Some of the functionality here has now (as of mid-2020) been implemented as marked below, but much remains to be done!

Goals

This project has several goals:

1. Provide a set of pandas ExtensionArrays for storing columns of discrete geometric objects (Points, Polygons, etc.). Unlike the GeoPandas GeoSeries, there will be a separate extension array for each geometry type (PointsArray, PolygonsArray, etc.), and the underlying representation for the entire array will be stored in a single contiguous ragged array that is suitable for fast processing using numba.
2. (partially done; see below) Provide a collection of vectorized geometric operations on these extension arrays, including most of the same operations that are included in shapely/geopandas, but implemented in numba rather than relying on the GEOS and libspatialindex C++ libraries. These vectorized Numba implementations will support CPU thread-based parallelization, and will also provide the foundation for future GPU support.
3. Provide Dask DataFrame extensions for spatially partitioning DataFrames containing geometry columns. Also provide Dask DataFrame/Series implementations of the same geometric operations, but that also take advantage of spatial partitioning. This would effectively replace the Datashader SpatialPointsFrame and would support all geometry types, not only points.
4. Provide round-trip serialization of pandas/dask data structures containing geometry columns to/from parquet. Writing a Dask DataFrame to a partitioned parquet file would optionally include Hilbert R-tree packing to optimize the partitions for later spatial access on read.
5. Fast import/export to/from shapely and geopandas. This may rely on the pygeos library to interface directly with the GEOS objects, rather than calling shapely methods.

These features will make it very efficient for Datashader to process large geometry collections. They will also make it more efficient for HoloViews to perform linked selection on large spatially distributed datasets.

Non-goals

spatialpandas will be focused on geometry only, not geography. As such:

• No built-in support for loading data from geography-specific file formats
• No dependency on GDAL/fiona
• No coordinate reference frame logic

Features

Extension Arrays

The spatialpandas.geometry package will contain geometry classes and pandas extension arrays for each geometry type

• Point/PointArray: single point / array of points, one point per element. Maps to shapely Point class.
• MultiPoint/MultiPointArray: multiple points / array of points, multiple points per element. Maps to shapely Points class.
• MultiLines/MultiLinesArray: One or more lines / array of lines, multiple lines per element. Maps to shapely LineString and MultiLineString classes.
• Ring/RingArray: Single closed ring / array of rings, one per element. Maps to shapely LinearRing class.
• Polygon/PolygonArray: One or more polygons, each with zero or more holes / array of polygons with holes, multiple per element. Maps to shapely Polygon/MultiPolygon classes.

Spatial Index

• The spatialpandas.spatialindex module will contain a vectorized and parallel numba implementation of a Hilbert-RTree.

• Each extension array has an sindex property that holds a lazily generated spatial index.

Extension Array Geometric Operations

The extension arrays above will have methods/properties for shapely-compatible geometric operations. These are implemented as parallel vectorized numba functions. Supportable operations include:

• area
• length
• bounds
• boundary
• buffer
• centroid
• convex_hull
• covers
• contains
• crosses
• difference
• disjoint
• distance
• envelope
• exterior
• hausdorff_distance
• interpolate
• intersection
• intersects_bounds
• minimum_rotated_rectangle
• overlaps
• project
• simplify
• union
• unary_union
• affine_transform
• rotate
• scale
• skew
• translate
• sjoin (spatial join) (see tools/sjoin) (partial, see intersection and spatial join support for Point arrays #21)

Only a minimal subset of these will be implemented by the core developers, but others can be added relatively easily by users and other contributors by starting with one of the implemented methods as an example, then adding code from other published libraries (but Numba-ized and Dask-ified if possible!).

Pandas accessors

Custom pandas Series accessor is included to expose these geometric operations at the Series level. E.g.

• df.column.spatial.area returns a pandas Series with the same index as df, containing area values.
• df.column.spatial.cx is a geopandas-style spatial indexer for filtering a series spatially using a spatial index.

Custom pandas DataFrame accessor is included to track the current "active" geometry column for a DataFrame, and provide DataFrame level operations. E.g.

• df.spatial.cx will filter a dataframe spatially based on the current active geometry column.

Dask accessors

A custom Dask Series accessor is included to expose geometric operations on a Dask geometry Series. E.g.

• ddf.column.spatial.area returns a dask Series with the same index as ddf, containing area values.

• The accessor also holds a spatial index of bounding box of the shapes in each partition. This allows spatial operations (e.g. cx, spatial join) to skip processing entire partitions that will not contain relevant geometries.

• A Custom dask DataFrame accessor is included that is exactly analogous to the pandas version.

Conversions

• Fast conversions to and from geopandas will be provided, potentially relying on pygeos.

Parquet support

• read/to parquet for Pandas DataFrames will be able to rely on the standard pandas parquet functions with extension array support.

• Special read/to parquet functions will be provided for Dask DataFrames.

• to_parquet will add extra metadata to the parquet file to specify the geometric bounds of each partition. There will also be the option for to_parquet to use Hilbert R-tree packing to optimize the partitions for later spatial access on read.

• read_parquet will read this partition bounding-box metadata and use it to pre-populate the spatial accessor for each geometry column with a partition-level spatial index.

Compatibility

• We would aim to use consistent naming between geopandas and spatialpandas whenever possible. Since spatialpandas will rely on a DataFrame accessor rather than a DataFrame subclass, spatial properties and methods are under the accessor rather than on the DataFrame.

For example, geodf.cx becomes spdf.spatial.cx.

• Eventually, a subclass compatibility layer could likely be added that would simply dispatch certain class methods/properties to the spatial accessor.

Testing

• The test suite will rely heavily on property testing (potentially with hypothesis) to compare the spatialpandas results to geopandas/shapely results.

[jonmmease]

@jbednar @jorisvandenbossche @philippjfr @jsignell

[brendancol]

@jonmmease this sounds awesome, and while it's a serious amount of work, at least GEOS already exists to be used as reference while implementing.

[brendancol]

Also, one anecdotal comment...

Two or three years ago, I attempted to implement an RTree using Numba and had trouble with numba optimizing trees with varying numbers of nodes in each level. It kept falling back to python mode and I gave up after the RTree couldn't beat a "brute force" query of the rectangles using numba in nopython mode.

3 thoughts:

1. My issue could have been user error (i.e. my error).
2. Numba may now be able to handle varying number of nodes in each level without falling back to python mode.
3. Numba may still have this limitation and @sklam may be a good person to comment on this thread.

[sklam]

@brendancol, can you point me to the old RTree implementation, or any current code that needs to be numba-ified?

[brendancol]

@brendancol, can you point me to the old RTree implementation, or any current code that needs to be numba-ified?

So its been a long time since I looked at this...and certainly don't want to waste your time. The associated code is here: https://github.com/parietal-io/py-rbush

We were probably running tests, profiling, and then just stopped mid-stream...

To boil down the approach, the RTree was using numba classes with deferred types. A simple question you may be able to help answer is:

"Since 2017, have there been significant developments in numba JITed classes and deferred types?"

@sklam Maybe we could think up a simple example (like the BTree example) which uses a varying number of nodes per level to help reproduce this issue and / or demonstrate an approach.

@jonmmease sorry to take over this thread...we can move this numba discussion somewhere else...

[jonmmease]

@brendancol, no problem. I'm actually playing around with a numba RTree implementation right now.

For the use cases I have in mind, I don't think there's a need to iteratively update the RTree after creation. In this case, I think I could get away with an array-based implementation of the binary tree (something like https://www.geeksforgeeks.org/binary-tree-array-implementation/). At least that's what I'm trying at the moment. No performance to report yet 🙂

[brendancol]

@jonmmease Awesome! I think part of my problem was probably "thinking by analogy"...Myself and @chbrandt were doing a port of https://github.com/mourner/rbush

[sklam]

Numba has support for a typed version of list, dict and heapq now. It would make things easier.

[jbednar]

Jon, this sounds like a fabulous proposal to me! Some musings:

• I think it would be good to lead the proposal (and thus any docs, blogs, home pages, etc.) with the goals that distinguish this project from other existing libraries. Otherwise, people can be confused and ask "why not just use PyGEOS", etc., until they reach the non-goals section. To my mind, those distinguishing and unique features are, in order:
  1. This project is explicitly independent of any assumptions about what space these spatial data structures and operations cover, in contrast to Shapely (which is "deeply rooted in the conventions of the geographic information systems (GIS) world" (even though it is often used outside of it)). To be useful for Datashader and other non-GIS applications, this project would need to have no such ties.
  2. As such, this project does not have any dependencies on the traditional geographic stack (GDAL, fiona, etc.), being built solely on top of general-purpose SciPy and PyData libraries and tools, which means that it should be easily installable using Conda or pip without any tricky binary-compatibility issues.
  3. All code included is pure Python compiled to machine code with Numba and able to run distributed and/or out of core using Dask (or Numba where appropriate), to allow it to be fully scalable for even the largest problems.
• At the same time, we need to make it clear that the traditional geo stack will always support much more functionality than this library ever will; this should be very clearly stated and acknowledged.
• As such, it is crucial that this library support fast and convenient conversion into data structures that can be used with the traditional stack, because it doesn't matter how fast its own operations are if users are forever having to spend code-development time filling in even the most trivial functions already available elsewhere.
• We need these basic capabilities for Datashader, and are willing to create the initial library and maintain it for Datashading purposes, but it should be clear that the Datashader maintainers are responsible only for the most basic function, and that any expansion of the supported functionality beyond that would be for the broader community to contribute. E.g. it seems to me that the initial library should support (and thus validate and test) only a small fraction of the geometric options listed, unless e.g. we could validate our tool against the existing tools' existing test suites. Most of the utilities involved are quite straightforward to implement and don't add up to much code, but they are well outside the purview of Datashader's maintainers, and so I'm reluctant for us to own the responsibility for documenting them, testing them, and tracking down any issues found in the tests.

[jonmmease]

@brendancol, here's the numba R-tree implementation I came up with today https://github.com/jonmmease/spatialpandas/blob/master/spatialpandas/spatialindex/rtree.py

I'll make up some benchmarks at some point, but what I'm seeing is that it matches brute force up to around 1000 elements, and then does significantly better the larger the inputs size gets. Several hundred times faster at a million elements. It also matches or is a bit faster than the rtree package, and its much faster to create the initial index than rtree. All of that is for a single lookup. What's really exciting is that these intersection calculations can be called from within numba functions. So I think it should be possible to write a really rast numba implementation of a spatial join on top of this.

[brendancol]

@brendancol, here's the numba R-tree implementation I came up with today https://github.com/jonmmease/spatialpandas/blob/master/spatialpandas/spatialindex/rtree.py

I'll make up some benchmarks at some point, but what I'm seeing is that it matches brute force up to around 1000 elements, and then does significantly better the larger the inputs size gets. Several hundred times faster at a million elements. It also matches or is a bit faster than the rtree package, and its much faster to create the initial index than rtree. All of that is for a single lookup. What's really exciting is that these intersection calculations can be called from within numba functions. So I think it should be possible to write a really rast numba implementation of a spatial join on top of this.

@jonmmease Sick!!! So excited to use this! Also another advantage to consider...libspatialindex isn't threadsafe and I've seen a bunch of seg faults related to rtree and multiprocessing. This will be super helpful for the Python GIS community!

Question: This has value outside of pandas. New standalone RTree library?

[jonmmease]

Question: This has value outside of pandas. New standalone RTree library?

Maybe eventually, sure. I'll probably want to wait to make sure it actually handles the usecases here before splitting it off though.

[jorisvandenbossche]

@jonmmease thanks a lot for writing up this proposal!

To repeat what I said in the datashader issue: I think there is certainly value in having a good, reusable implementation of geometry data support with a ragged array representation.
Below already some comments on the goals (from a geopandas developer, so somewhat critical :-)), will try write down more comments on the rest of the discussion later.

About the goals:

Goals 1 (ragged-array based extension arrays) and 2 (numba-based algorithms) are IMO distinct goals for a project like this.
But goal 3 (dask extension to have spatially partitioned DataFrames) and 4 (round-trip serialization with parquet making use of spatial index) are perfectly valid (and realistic) goals for GeoPandas as well (it's also what was written in the SBIR project, and more or less what the POC in https://github.com/mrocklin/dask-geopandas/ is doing (partially)).
And moreover, I think implementations of those, could be rather ignorant of the actual representation of the geometries under the hood.

So I would love to see or explore if we can build common infrastructure for this, instead of each of the projects doing this themselves.

One possible way (but need to think this through more) is to make the dask extensions, IO functionality, etc to work with different kinds of GeometryArrays (if they can rely on the public methods of those, that might be possible). But will think about this aspect some more.

[jbednar]

Just to clarify, goal 4 is only in service of the other goals -- anything we work on needs to have some sort of persistence/serialization format, so that it's practical to work with, and the goal is not to provide one (if at all possible), just for there to be one. So it would be ideal if there is such a format available that both GeoPandas and this non-GeoPandas code can both use.

In any case, I'm very excited about the potential for having algorithms, persistence formats, etc. that can work with both GEOS and a Python-native, Numba-izable object, which would be amazing if we can achieve that.