Partitioning GeoParquet by H3 Spatial Index

Partitioning a GeoParquet archive by H3 cell turns a directory of opaque columnar files into a spatially prunable dataset, where a query engine reads only the partitions whose hexagons intersect a bounding box instead of scanning every object in cold storage. This guide is for data engineers and GIS archivists who need a deterministic, cold-storage-aware partition scheme: choosing an H3 resolution that balances predicate-pushdown selectivity against a cloud provider’s minimum-object-size billing, deriving stable cell keys from geometry centroids, and writing a Hive-style h3_cell= layout that DuckDB, GDAL, and pyarrow can prune without opening files. Default writers ignore this entirely — they emit a flat file set with no partition column, so every spatial filter degrades into a full-archive scan. It operates under the Format Conversion & Pipeline Automation framework and extends the batch GeoParquet Migration Workflows reference with a partition-key strategy tuned for archival economics.

Why H3 for Archive Partitioning #

H3 is a hierarchical hexagonal grid: every cell has exactly one parent at the next-coarser resolution, so a partition scheme built on it is reprojectable up the hierarchy without recomputing keys from geometry. Hexagons also have uniform neighbor distance — unlike a quadkey grid, adjacent cells are equidistant, which keeps a radius query from fetching a lopsided set of partitions. The single decision that governs everything downstream is resolution: it fixes cell area, and therefore the feature count and byte size of each partition file.

The resolution you pick is a storage-economics decision as much as a query one, so make it against the same cold-tier minimums that drive Hot/Warm/Cold Tier Design for Geospatial Data — pick the coarsest resolution whose median file still clears the tier’s minimum billable object size.

Sizing the Resolution Against Cold-Storage Minimums #

Cold classes bill a minimum object size and a minimum retention duration regardless of the bytes you actually store. AWS Glacier Deep Archive rounds every object up to a 128 KB billable floor and adds a per-object overhead, so a partition scheme that shatters a dataset into millions of sub-megabyte files pays overhead on every one and cripples restore throughput. The countervailing pressure is pushdown: coarser cells mean fewer, larger files, so a tight bounding-box query still drags in gigabytes of irrelevant features.

Compute the resolution empirically rather than guessing. Sample the feature density of the source, project a per-cell byte estimate at each candidate resolution, and pick the coarsest one whose median file lands in your target range (a common archival target is 200 MB–1 GB per object, large enough to amortize per-object overhead, small enough to parallelize restore).

import h3
import geopandas
import numpy as np

# Sample 1% of a large source to estimate per-cell feature density.
sample = geopandas.read_parquet("s3://spatial-archive/staging/parcels_sample.parquet")
centroids = sample.geometry.to_crs(4326).centroid
avg_bytes_per_feature = 320  # measured: WKB geometry + attribute columns, post-ZSTD

for res in (3, 4, 5, 6):
    cells = [h3.latlng_to_cell(pt.y, pt.x, res) for pt in centroids]
    counts = np.array(list({c: cells.count(c) for c in set(cells)}.values()))
    # scale the 1% sample up to the full population
    est_features = counts * 100
    est_mb = est_features * avg_bytes_per_feature / 1_048_576
    print(f"res={res}: partitions={len(counts):>6d}  "
          f"median_file={np.median(est_mb):8.1f} MB  p95={np.percentile(est_mb,95):8.1f} MB")

Read the output as a curve: resolution is right when the median file clears your cold-tier floor and the p95 file stays under the restore-parallelism ceiling you can tolerate. If the p95 is an order of magnitude above the median, the data is spatially skewed (dense urban cells, empty ocean cells) and a single global resolution will not fit — see the mitigation below.

Deriving Stable Cell Keys and Writing the Layout #

Derive the H3 index from the geometry centroid in EPSG:4326, because H3 is defined on the WGS84 sphere and latlng_to_cell expects geographic coordinates. Deriving from a projected CRS silently places features in the wrong cell. Pin the CRS explicitly — the same discipline enforced by CRS Synchronization in Pipelines — and store the resolution in the archive manifest so future writers reproduce identical keys.

import json
import geopandas
import pyarrow as pa
import pyarrow.parquet as pq

H3_RES = 4  # pinned in the archive manifest; never change without a full repartition

def write_h3_partitioned(src: str, dst_root: str, crs_epsg: int = 4326):
    gdf = geopandas.read_parquet(src)
    # Centroid in EPSG:4326 is the ONLY correct input to H3.
    centroids = gdf.geometry.to_crs(4326).centroid
    gdf["h3_cell"] = [h3.latlng_to_cell(p.y, p.x, H3_RES) for p in centroids]

    geo_meta = {
        "version": "1.1.0",
        "primary_column": "geometry",
        "columns": {"geometry": {"encoding": "WKB",
                                  "geometry_types": sorted(gdf.geom_type.unique().tolist()),
                                  "crs": f"EPSG:{crs_epsg}"}},
    }
    gdf["geometry"] = gdf.geometry.to_wkb()

    for cell, part in gdf.groupby("h3_cell"):
        table = pa.Table.from_pandas(part, preserve_index=False)
        table = table.replace_schema_metadata({b"geo": json.dumps(geo_meta).encode()})
        # Hive-style directory: the engine reads h3_cell from the path, not the file.
        out = f"{dst_root}/h3_cell={cell}/part-0000.parquet"
        pq.write_table(table, out, compression="zstd", compression_level=3,
                       row_group_size=100_000)

The written tree is Hive-partitioned: s3://spatial-archive/parcels/h3_cell=8428309ffffffff/part-0000.parquet. The h3_cell value lives in the path, not inside the file, so a partition-pruning engine skips whole directories before it opens a single Parquet footer. Note the h3_cell column is dropped from the file body by convention — it is redundant with the path — which also shrinks each object.

For skewed data where one resolution cannot serve both dense and sparse regions, partition at a coarse base resolution and let row-group sizing handle intra-partition selectivity instead of pushing to a finer H3 level; the Row Group Sizing Strategies reference covers picking a row-group size that keeps in-file pruning effective on multi-gigabyte partitions.

Verifying Predicate Pushdown on the Partition Column #

Writing the layout is only half the job; confirm that engines actually prune it. DuckDB with hive_partitioning=1 exposes h3_cell as a virtual column and skips non-matching directories.

duckdb -c "
  SET enable_object_cache=true;
  EXPLAIN ANALYZE
  SELECT count(*) FROM read_parquet(
      's3://spatial-archive/parcels/**/*.parquet', hive_partitioning=1)
  WHERE h3_cell IN ('8428309ffffffff','842830bffffffff');"

The EXPLAIN ANALYZE output must show that only the matching partitions were read. Look for the Files Read counter — it should equal the number of cells in your IN list, not the archive total:

┌─────────────────────────────────────┐
│           PARQUET_SCAN              │
│   Files Read: 2 / 18443             │  ← pruned to the two matching cells
│   Rows Scanned: 214,880            │
│   Total Time: 0.31s                │
└─────────────────────────────────────┘

If Files Read equals the archive total, pruning failed — the filter did not reference the partition column, or hive_partitioning was off. To resolve a bounding box to the set of cells to query, use h3.polygon_to_cells on the query envelope and pass the result as the IN list; this is how you translate a map viewport into a partition predicate.

Troubleshooting Partition Pruning #

Operational Execution Checklist #

• Estimate per-cell feature density from a 1% sample and pick the coarsest H3 resolution whose median file clears the cold-tier minimum object size.
• Pin the chosen H3_RES in the archive manifest; treat any change as a full repartition, not an incremental edit.
• Derive every cell key from the geometry centroid reprojected to EPSG:4326 — never from a projected CRS.
• Write a Hive-style h3_cell=<index>/ tree and drop the redundant h3_cell column from the file body.
• Preserve the geo metadata block (version, primary_column, WKB encoding, CRS) on every partition write.
• Verify pruning with EXPLAIN ANALYZE and confirm Files Read matches only the intersecting cells.
• Translate query bounding boxes to partition predicates with h3.polygon_to_cells, not geometry-level ST_Intersects on the whole archive.

Related #

• Up: GeoParquet Migration Workflows — the parent reference for archive-wide partitioning and metadata heuristics this scheme plugs into.
• Converting Legacy Shapefiles to GeoParquet at Scale — the batch conversion job that generates the H3 keys written here.
• Writing GeoParquet Bounding-Box Metadata for Predicate Pushdown — the intra-file complement to partition pruning, skipping row groups after the directory prune.
• Calculating Optimal Row Group Size for Spatial Queries — sizing row groups so pruning still works inside a skewed, oversized partition.
• Spatial Archive Cost Modeling — the minimum-object-size and per-request math that sets the floor on your H3 resolution.

Part of the Spatial Data Archival knowledge base.