PostGIS to partitioned GeoParquet

[jordanp-dluhc]

I've been following the discussion on GeoParquet partitioning strategies and came across @kylebarron's comment about this library which might be of use.

We're looking to convert versions of somewhat large datasets to spatially partitioned GeoParquets. Each vary in size, for example the point dataset is ~150GB in postgres, whereas the polygons are ~1-3TB.

The datasets are available in annual snapshots going back ~20 years. The way people query it is usually temporally then spatially (i.e. Find me data in 2018 for this area).

We've been looking at various approaches such as KD-Trees with DuckDB mentioned in the above discussion.

I'm wondering if there's a way to construct the partitions straight from PostGIS using the PostGIS and GeoParquet functions?

Any thoughts or pointers would be amazing, in the meantime I'll look to construct a reproducible example if anyone is interested.

[kylebarron]

Hello 👋

I've been exploring spatially partitioned GeoParquet for a while. And it's something I'm motivated to get integrated here, but in the short term the answer to your question is no, this isn't fully ready yet via this library.

I'm getting ready for a new Rust and Python release in January, and then focusing more on integrating spatial support for DataFusion, an extensible query engine built on Rust, Arrow, and Parquet. I'm hoping that this integration will enable spatial partitioning of large datasets easy to do out of the box.

I'd also point you to https://github.com/kylebarron/geo-index, which is a very memory efficient way to create a KD tree and which has Python bindings. So it may be possible to create the KDTree in memory, even if you can't put the full dataset in memory. In the near term I'd like to add APIs to that which let you extract the spatial partitioning information from the generated kdtree. (Also, I'm not sure why but in that linked discussion, they seemed to say that generating the KDTree partitioning information took a few minutes with the recursive SQL approach, but with geo-index it only took 13 seconds from Python).

I'd like to have better spatial partitioning functionality in place here (and ideally easy to use) in Q1-Q2 2025. Contributions are also welcome!

If your data is already in PostGIS, I'd start by trying to construct the output partition information using PostGIS and then grouping by each of those output partitions and writing those to disk directly from PostGIS.

[jordanp-dluhc]

Hi Kyle, thanks for the quick reply! Apologies for the lateness of mine.

As a follow up, I managed to play around with the geo-index python bindings and the partitioning APIs you released.

Was simple to read from PostGIS with the geoarrow-rs python bindings, build an R-Tree and write to partitioned files with ogr2ogr following the general GeoParquet distributing guidance docs.

Took ~5 minutes on a test case of ~30million polygons, the bulk of which was from ogr2ogr. I think the geo-index R-tree took <1 second which is awesome.

I couldn't figure out the syntax to calculate the bounding box extents using geoarrow.rust.compute.envelope, so if I figured that out I imagine it'd be a little quicker.

Very interested to see how spatial support in DataFusion comes along, seems really promising.

[kylebarron]

That's awesome! Is any of the code you wrote public?

[jordanp-dluhc]

It was actually a test case using Overture data so no reason why it can't be, give me maybe a week or so in between my day job to tidy a few bits up and I'll send over a gist.

[jordanp-dluhc]

Just to check, what I did previously when creating the geo-index R-Tree was:

• read the data as a geoarrow table
• convert to shapely
• calculate bounds
• pass bounds to R-Tree

As below:

import geoarrow.rust.core as geoarrow_core
import geoarrow.rust.io as geoarrow_io
import shapely
from geoindex_rs import rtree as rt

table = geoarrow_io.read_postgis(
    connection_url="postgres://spatial_user:my_password@localhost:5432/spatial_db",
    sql="select id, geom from buildings_eu",
)

geometries = geoarrow_core.to_shapely(table["geometry"])

bounds = shapely.bounds(geometries)

tree_size = 65_535
builder = rt.RTreeBuilder(len(geometries), tree_size)

builder.add(bounds)

tree = builder.finish()

This works OK for moderately sized data (~24 million) but anything bigger and it crashes.

Am I right in assuming geoarrow.rust.compute.envelope is meant to return the bounding box of each geometry?

If so, I'm using it like this:

import geoarrow.rust.compute as geoarrow_compute
import geoarrow.rust.io as geoarrow_io

table = geoarrow_io.read_postgis(
    connection_url="postgres://spatial_user:my_password@localhost:5432/spatial_db",
    sql="select id, geom from buildings_eu",
)

geoarrow_compute.envelope(table["geometry"])

Which returns the following errors:

thread '<unnamed>' panicked at [/Users/runner/.cargo/registry/src/index.crates.io-6f17d22bba15001f/arrow-array-53.3.0/src/array/struct_array.rs:91:46](https://file+.vscode-resource.vscode-cdn.net/Users/runner/.cargo/registry/src/index.crates.io-6f17d22bba15001f/arrow-array-53.3.0/src/array/struct_array.rs:91:46):
called `Result::unwrap()` on an `Err` value: InvalidArgumentError("Incorrect number of arrays for StructArray fields, expected 4 got 8")

Am I missing something completely obvious?