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

[Maxxen]

Yes! It's mostly a limitation of how the query engine is structured. You have to make a pass over every input geometry to calculate the bounds and then another pass to calculate the Hilbert code now that the bounds are known, but that's not possible with a normal scalar function. And even if it was, because DuckDB is a streaming engine you'd have to effectively stop the query pipeline and buffer the whole queryset (potentially to/from disk) which would be pretty expensive.

If extensions could provide custom statistics we could potentially keep the known bounds around for each column, but that's not possible yet. However, I have some ideas on how to manage/track SRID/projection info for each column id like to get to in the near future, and depending on how that goes we might be able to automatically infer the bounds as the bounds of the coordinate system.

[Maxxen]

Also, there isn't really a way in SQL to pass aggregates as input to scalar functions without separating it into multiple statements (first compute the aggregate result, e.g. st_extent_agg) and then pass it to the scalar (st_hilbert). So we can't create a macro or query-rewrite-level shortcut either.

[cholmes]

Ah, got it. That all makes a lot of sense. Could it be a mini-extension? Like is there a way with an extension to just have it run the two statements? Anyways, I'm happy to document it in the geoparquet publishing best practices I'm working on. And would be happy to contribute to the docs or a make a little tutorial. It's pretty simple, but it's not obvious at all from the info I could find.

[Maxxen]

The closes I can think of is a table macro, like this

CREATE OR REPLACE FUNCTION hilbert_auto(col_name, table_name) 
AS TABLE SELECT 
  st_hilbert(COLUMNS(col_name), (
    SELECT st_extent(ST_Extent_Agg(COLUMNS(col_name)))::BOX_2D FROM query_table(table_name))
  ) FROM query_table(table_name);

You'd use it like

SELECT * FROM hilbert_auto('geom', 't1');

This queries the table given by 'table_name' twice, once to compute the bounds and once to get the Hilbert values.

It's pretty simple, but it's not obvious at all from the info I could find.

Yeah this would for sure be good blog-post material! I want to add more guide/tutorial style content to docs at some point as well but I want to get some more stuff ready code-side before digging into that.

[cholmes]

The hilbert_auto didn't quite work for me, as it seemed to just return a single column. I asked chatgpt and got:

SELECT *
  FROM fsq
  ORDER BY st_hilbert(
      geom,
      (
          SELECT st_extent(ST_Extent_Agg(COLUMNS(geom)))::BOX_2D
          FROM fsq
      )
  )

Which seemed to work? But only tested on ~5 million rows. Not sure if there's something wrong with that approach? It also suggested:

WITH bounds_cte AS (
  SELECT st_extent(ST_Extent_Agg(COLUMNS(geom)))::BOX_2D AS bounds
  FROM fsq
)
SELECT f.*
FROM fsq AS f
ORDER BY st_hilbert(
    f.geom,
    (SELECT bounds FROM bounds_cte)
);

which I didn't try out, and it didn't seem more readable to me.

Yeah this would for sure be good blog-post material!

Cool, I'll try to throw something up. Let me know which of the above to do, or else another way to do it with two lines of SQL.

I want to add more guide/tutorial style content to docs at some point as well but I want to get some more stuff ready code-side before digging into that.

Whenever you start this I'd be happy to try to add a few sections for some of the things I do often.