Spatially sort output features

[brawer]

As per GeoParquet best practices it would be good to spatially sort the output features, for two reasons:

• Typically, clients will restrict their parquet queries to a region of interest. Without spatial sorting, most parquet chunks contain data from locations all over the planet. The current lack of spatial correlation means that clients have to decompress almost all chunks in the parquet file, no matter their query bounding box — this is expensive. If the data was spatially sorted, queries would be (much) faster because they’d only have to decode those few parquet chunks that actually intersect the queried bounding box.
• Also, spatial sorting will likely reduce the output file size. Because nearby features tend to share tags like street and city names, there will be a higher chance of sharing tags between features of a single parquet chunk.

Currently, class GeoParquetWriter seems to emit features in the same order as they happen to be passed from libosmium. Consider extending the implementation of GeoParquetWriter to calculate the center lat/lon of each feature’s bounding box. Then, find the position of that point along a space-filling Hilbert curve, and use this number as a sort key for an external sort. There's several python libraries for hilbert curves, and likewise for external sorting.

To check the difference, try gt sort hilbert from GeoParquet tools. Perhaps you could simply call this tool in a post-processing step, before uploading the layercake output. But it seems a little heavy to bundle DuckDB; doing this yourself from python seems easy enough.

[ianthetechie]

I had the same thought. But I think there are a few things that contribute to this being tricky.

OSM data files are "special" ;) They are a collection of lists of different primitives with the nodes first, then the ways, then the relations. You have to make two passes through the file (or build up a huge cache) to make this efficient.

Spatial clustering of any sort is a non-goal of the OSM file format. And by extension, PyOsmium, which powers this tool, will not be processing anything in a spatially optimized order.

IMO delegating this to DuckDB is a perfectly reasonable option. I'm not really sure how "easy" this would be from Python, unless you're only looking at tiny extracts. Doing this in (pure) Python is not going to be very efficient.

[jake-low]

I agree this is a very good idea, and I'm looking into how to implement this. Layercake's runtime and operating costs are closely related to how much RAM is required for processing, so I'm trying a couple of approaches to see which one is the most memory efficient.

[ianthetechie]

Here are some anecdotal findings using the latest Asia extract and extracting boundaries...

An "unoptimized" parquet weighs in at 1.35GB. There is essentially no compression happening at the moment for a variety of reasons, including both sorting, tiny row groups, and some "compatibility" defaults at the writer level like snappy compression instead of zstd.

After the following SQL (DuckDB), I get a 674MB parquet ;)

INSTALL spatial;
LOAD spatial;

SET VARIABLE input_parquet = 'out/boundaries.parquet';

COPY (
    SELECT *
    FROM read_parquet(getvariable('input_parquet'))
    ORDER BY ST_Hilbert(
        geometry,
        (
          SELECT ST_Extent(ST_Extent_Agg(geometry))
          FROM   read_parquet(getvariable('input_parquet'))
        )
    )
)
TO 'out/boundaries-optimized.parquet' (
    FORMAT PARQUET,
    COMPRESSION ZSTD
    -- TBD: Partitioning; row group size by bytes?
);

[ianthetechie]

Layercake's runtime and operating costs are closely related to how much RAM is required for processing

Yeah, that's a really good point. I'll also post findings as they arise either here or in another ticket. One of my first counter-intuitive discoveries was that using the dense_mem_array index mode made it run approximately 3x faster on my MacBook, despite appearing to use roughly the same amount of RAM on my MacBook.

I haven't dug into a detailed trace or anything, but thought I'd share that haha.

[missinglink]

FWIW, GDAL has a SORT_BY_BBOX option to handle this:

https://gdal.org/en/latest/drivers/vector/parquet.html

SORT_BY_BBOX=[YES​/​NO]: (GDAL >= 3.9) Defaults to NO. Whether features should be sorted based on the bounding box of their geometries, before being written in the final file. Sorting them enables faster spatial filtering on reading, by grouping together spatially close features in the same group of rows.

Note however that enabling this option involves creating a temporary GeoPackage file (in the same directory as the final Parquet file), and thus requires temporary storage (possibly up to several times the size of the final Parquet file, depending on Parquet compression) and additional processing time.

The efficiency of spatial filtering depends on the ROW_GROUP_SIZE. If it is too large, too many features that are not spatially close will be grouped together. If it is too small, the file size will increase, and extra processing time will be necessary to browse through the row groups.

Note also that when this option is enabled, the Arrow writing API (which is for example triggered when using ogr2ogr to convert from Parquet to Parquet), fallbacks to the generic implementation, which does not support advanced Arrow types (lists, maps, etc.).

[ianthetechie]

Note also that when this option is enabled, the Arrow writing API (which is for example triggered when using ogr2ogr to convert from Parquet to Parquet), fallbacks to the generic implementation, which does not support advanced Arrow types (lists, maps, etc.).

This last bit would be a major bummer. I am pretty confident that we DO want to use the "advanced" (amusing wording) Arrow types. Specifically, I'm not sure of any better way to represent localized names than via a map. We certainly aren't going to explode dozens of columns, and it'd be a bummer of we end up "deciding" that a language just isn't important. Maps provide a really convenient solution to that, which would let the consumer decide which language(s) they want to ingest.

Fortunately, the DuckDB implementation seems to be more "arrow-native."

[missinglink]

Potentially stupid question but.. do you need to choose a single primary indexing method per file (ie. either bbox or primary ID / row ID)?

I agree that doing spatial filtering is a killer feature, but I guess we would lose the ability to do ID based lookups in O(~1)?

[missinglink]

Or is it that since it's a columnar format that you can do most lookups in O(log n) but the bbox index is 'special' because it's a composite index of 4 fields?

[ianthetechie]

Possibly the blind leading the blind here, but my understanding is that Parquet does not actually have a concept of an index per se; only row group stats and a file footer.

Thus, it's easy to effectively index over a single column, and there are ways to effectively spatial index (i.e. hilbert curve), but there is not a general way to index over more than one column. Thus, yeah, you effectively do have to choose your "index" method once and somewhat carefully based on how your data will be accessed. It's one of the tradeoffs of Parquet.

[missinglink]

Yeah, after thinking about it some more, maybe what's happening here is that, while it's fast enough to enumerate a list of Row IDs which match the query, the subsequent operation of loading the rest of the column data would be significantly faster if they were all next to each other in the data blocks.

[missinglink]

@brawer mate, do you understand how this works? 😆

[brawer]

How Parquet works? There's a nice description of the format here, have a look. Also, have a look at this pull request on the GeoParquet spec.

[jtbaker]

Just a datapoint on the DuckDB postprocessing approach:

Ran a Hilbert sort on a moderately spec'd laptop (M4 MB Air, 24GB RAM) with these configurations (added the memory_limit and preserve_insertion_order options after a few OOMs):

duckdb ./file.db
INSTALL spatial; LOAD spatial;

D SET memory_limit='14GB';
D SET preserve_insertion_order=false;
D COPY(SELECT
     * FROM 'highways.parquet' 
     ORDER BY st_hilbert( geometry, ( SELECT st_extent(ST_Extent_Agg(geometry))::BOX_2D FROM './highways.parquet' ) )
   ) TO 'highways_sorted.parquet'(FORMAT 'parquet', COMPRESSION 'zstd');

Took about 2 hours, and cut the file size by ~1/3. I'm not sure how much of that is due to the compaction of the row groups and being closer/sharing more properties, or the ZSTD compression (huge fan!). It's possible copying the parquet file to a real duckdb table first to operate on may speed that up some as well.