Is there a pure DuckDB way to order GeoParquet output spatially?

[cholmes]

So with DuckDB 1.1 there's GeoParquet output, with bbox struct (awesome), and r-tree indexes (also awesome). But I loaded some geospatial data that wasn't ordered and then wrote out the geoparquet and it's clearly not indexed / ordered at all. Is there a way to tell DuckDB to order the data by the r-tree? I tried creating an r-tree index and then doing 'order by geometry', but it didn't seem to do anything.

My default way would be to create a quadkey and order by that - but it seems like the native r-tree should be able to help me and hopefully remove a couple steps.

Happy to contribute to the docs / share broadly if there's a way to do it.

[Maxxen]

Hello! Unfortunately, DuckDB does not support overriding the comparison operators (used to sort) for custom types, meaning ORDER BY <geometry> will just order based on the lexical order of the bytes in the internal binary representation. The r-tree index is completely separate from ordering. I know that other databases (like Postgres) sometimes can perform an index scan instead of explicitly sorting when executing an order-by, but this requires the ordering criteria to map to the physical order of the index. In the case of R-Trees, I can't think of any way that would work.

Perhaps there could be a "dummy" scalar function the optimizer could detect to do this if the column is r-tree indexed, or otherwise error out, e.g. SELECT geom FROM tbl ORDER BY rtree_order(geom). But again, the r-tree traversal order might not actually be the most efficient to use for spatial clustering, compared to e.g. a Hilbert curve. Additionally, doing index lookups for every row in the table is going to be significantly more expensive than a full table scan. DuckDB indexes don't store any actual data within themselves (unlike how indexes sometimes work in other databases), but just the row_id of the row in the indexed table. This means that every time an index scan is performed, DuckDB first has to traverse the index, collect the row ids, and then perform (what is most likely) random-access point-lookups from the underlying table row-by-row.

There are also some other issues related to the preservation (or lack thereof) of row-order when scanning from a table in parallel or writing to parquet. I think the only way to guarantee the output order of the rows within parquet files is to execute a ORDER BY, which will emit chunks with a "batch index", that the parquet writer can use to reorganize incoming chunks in flight and ensure they are written in sorted order, even when executing with multiple threads.

So uh, yeah, we could look into adding some sort of dummy marker function + optimizer rule + index-batch-scan operator, but I don't think its going to be much faster or produce much better results than just doing a full table scan + explicit sort, e.g. COPY (SELECT geom FROM tbl ORDER BY ST_Hilbert(geom) ) TO 'gpq.parquet

[cholmes]

Cool. That's all great to know. I think all I was really looking for was your last line - I didn't realize there was ST_Hilbert.

[Maxxen]

By default it maps the entire space of two floating point coordinates into a single uint32, which produces kind of messy Hilbert curves, but it is overloaded so you can pass an optional "bounds" parameter which will normalize the centroid of the geometry to 0-1 based on how close it maps to the boundary, which gives the resulting Hilbert curve better precision.