GeoParquet Spatial Sorting/Indexing/Partitioning

[jwass]

Now that Overture is publishing GeoParquet, we should start thinking about how to sort or order the data within Parquet datasets in order to accelerate spatial queries.

There's a lot of recent interest in this work and it feels somewhat new.

Some applicable background:

• Chris Holmes explorations using admin partitioned geoparquet datasets: https://medium.com/radiant-earth-insights/the-admin-partitioned-geoparquet-distribution-59f0ca1c6d96
• Chris also did some work to build a quadtree by recursively partitioning buildings data until each partition reached some max. size: https://github.com/opengeos/open-buildings/blob/f1703b286ca54f3a22a53d2e96ac4d229962f65a/open_buildings/overture/partition.py#L132
• My proposal for GeoParquet to define a standard to specify geometry bounding boxes: Defining Bounding Boxes for Optimizing Spatial Queries opengeospatial/geoparquet#188

Let's use this discussion to demonstrate some progress, brainstorm ideas/implementations, etc.

[mtravis]

@jwass thanks for kicking off this discussion. I'm really keen to see partitioning by country_iso for buildings and places. I think it would really help with developer advocacy and the general experience for Overture users who just want to get their hands on the data. I think what Chris has done so far with releasing the buildings is a fine example of what could be doing.

I also like the idea of using bbox per row and hope that does get adopted as part of the geoparquet standard.

Maybe we could start with buildings and see if we can add a coutry_iso per row there to use as the partition? In terms of the outputted geoparquet data from Overture, is that currently created using Apache Sedona?

[FengJiang2018]

Yes. Currently the geoparquet data was created by using Apache Sedona.

[jwass]

I'm personally less-inclined to split it all up by country. This does create a more nested directory structure with one directory per country. And we have to do it this way across the entire dataset. We also need to think through features that might not have a country.

My current thinking is that we can possibly just add a column for the country. If we get the spatial indexing right, it might not be so bad to filter on a country. But we should measure.

[jwass]

Also came across https://github.com/kylebarron/spatially-partitioned-geoparquet as an additional reference.

[jwass]

See @marklit's detailed comment #113 (comment)

[jwass]

I made some significant improvements to the spatial partitioning in our release data.

Here are two images showing the buildings data for January vs. what we released this morning. On the left of each image is the bounding box of each file. On the right is the bounding box of each row group

January buildings data:

March buildings data:

This is a basic range partitioning of the geohash. Each file is then sorted by geohash 15. I think there's still plenty of room for improvement but this is a basic approach that should buy significant improvements in query performance.

[saosebastiao]

This is incredible work, thank you so much. Will this be used in an upcoming release?

[mtravis]

I believe its in the newest release Overture 2024-03-12-alpha.0

I notice some significant speed up times on spatial queries when using the bbox

[marklit]

Fantastic work!

[jwass]

@saosebastiao Yup - this is in the March release that was out yesterday.

[Spothedog1]

Howdy, let me know if I should post this somewhere else. I was trying to download all the buildings with a bbox that represents Washington D.C. I ran a DuckDB query within a Python script to execute it, but after 10 hours it's still not done. I've tried multiple times for a few hours each time but it's not working.

I'm wondering if something in my query could be better, or if its possible my local desktop just can't fit all of it into memory. With the improved partitioning I assumed it'd be sort of fast. I'm running this on a 4 year old mac with 16gb memory.

Any help would be appreciated.

INSTALL spatial;
INSTALL httpfs;
LOAD spatial;
LOAD httpfs;

SET s3_region='us-west-2';

COPY (
    SELECT
        *
    FROM read_parquet('s3://overturemaps-us-west-2/release/2024-03-12-alpha.0/theme=buildings/type=building/*',
                        filename=true, hive_partitioning=1)
    WHERE 
        bbox.minX BETWEEN -77.119759 AND -76.909393 AND
        bbox.minY BETWEEN 38.791645 AND 38.995548
        
)
TO '/Volumes/external/overture/2024-03-12-alpha.0/buildings_building.parquet'
WITH (FORMAT 'parquet');

[marklit]

Was this with DuckDB 0.10 which has the parquet data struct predicate patches or was it with an older version?

The 45th PQ file is the only one with DC geometry.

SELECT
          filename, COUNT(*)
      FROM read_parquet('s3://overturemaps-us-west-2/release/2024-03-12-alpha.0/theme=buildings/type=building/*',
                          filename=true, hive_partitioning=1)
      WHERE
          bbox.minX BETWEEN -77.119759 AND -76.909393 AND
          bbox.minY BETWEEN 38.791645 AND 38.995548
      GROUP BY 1
      ORDER BY 1
  ;

┌────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┬──────────────┐
│                                                      filename                                                      │ count_star() │
│                                                      varchar                                                       │    int64     │
├────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┼──────────────┤
│ s3://overturemaps-us-west-2/release/2024-03-12-alpha.0/theme=buildings/type=building/part-00045-4dfc75cd-2680-4d…  │       331395 │
└────────────────────────────────────────────────────────────────────────────────────────────────────────────────────┴──────────────┘

This finished in 16 seconds for me:

COPY (
    SELECT
        *
    FROM read_parquet('s3://overturemaps-us-west-2/release/2024-03-12-alpha.0/theme=buildings/type=building/*00045*',
                        filename=true, hive_partitioning=1)
    WHERE 
        bbox.minX BETWEEN -77.119759 AND -76.909393 AND
        bbox.minY BETWEEN 38.791645 AND 38.995548
)
TO 'dc.parquet'
WITH (FORMAT 'parquet');

select count(*) from read_parquet('dc.parquet');

┌──────────────┐
│ count_star() │
│    int64     │
├──────────────┤
│       331395 │
└──────────────┘

.mode line
 select * from read_parquet('dc.parquet') limit 1;

              id = 08b2aa8793266fff0200cca5e44ad092
        geometry = \x00\x00\x00\x00\x03\x00\x00\x00\x01\x00\x00\x00\x09\xC0SG\xA10f\xD0\x19@CeYu\xBD\x1B\x7F\xC0SG\xA2\x19\x9A\xDB\xB9@CeYA\xBA\xB5s\xC0SG\xA2@1Hw@Ce[\x14\xF9\x8CH\xC0SG\xA3\xEC\xD9\xB2\x8E@CeZ\xBA`\xB9~\xC0SG\xA3v\xFD\x8D_@CeUd\xB6b\xFE\xC0SG\xA1\xB0\xBFP\x0D@CeU\xC1\xD3t\x87\xC0SG\xA1\xC5L\xA5\xCB@CeV\xC41\xF3\x98\xC0SG\xA0\xF9\x09k\xBA@CeV\xEB\x9F\x1F\xEB\xC0SG\xA10f\xD0\x19@CeYu\xBD\x1B\x7F
            bbox = {'minx': -77.1193802, 'maxx': -77.1192, 'miny': 38.7916685, 'maxy': 38.7918421}
         version = 0
     update_time = 2022-09-08T18:53:18.000Z
         sources = [{'property': , 'dataset': OpenStreetMap, 'record_id': w835174811@2, 'confidence': NULL}]
           names =
           class =
           level =
       has_parts = false
          height = 3.8
      num_floors =
    facade_color =
 facade_material =
   roof_material =
      roof_shape =
  roof_direction =
roof_orientation =
      roof_color =
     eave_height =
        filename = s3://overturemaps-us-west-2/release/2024-03-12-alpha.0/theme=buildings/type=building/part-00045-4dfc75cd-2680-4d52-b5e0-f4cc9f36b267-c000.zstd.parquet
           theme = buildings
            type = building

[Spothedog1]

@marklit Ah thanks, my DuckDB CLI version was 0.10 but my Python DuckDB version was 0.9. I just upgraded to 0.10, hopefully that fixes it. Thanks so much.

EDIT: Wow yea that was it, finished in a few seconds. Fixed

[ivanthewebber]

sorted by geohash 15

It seems like sorting by a Hilbert curve would be ideal. Then you should get row groups that are pretty contained without any special logic.

I assume your geohash is something like Bing's quadkeys, Metropolis Grid, or Google's S2. S2 is already a Hilbert curve, and the others can be mapped onto it.