Partitioning Cold Spatial Archives by Administrative Boundary

Partitioning a cold spatial archive by administrative boundary — country, then state or province, then county — aligns the physical file layout with how governance, retention law, and analyst queries actually slice the data, at the cost of the wild size skew that uneven feature density inflicts on grid-free layouts. This guide is for the data engineers, GIS archivists, and compliance teams deciding whether boundary-aligned partitions beat a uniform spatial grid, and shows exactly how to build a hybrid admin-plus-grid tree that keeps sparse rural jurisdictions and dense metros both within a workable file-size band. It operates under the Compression Tuning & Storage Optimization framework, where partition layout drives predicate pushdown, per-request billing, and the blast radius of a retention deletion. Naive “one file per county” partitioning fails the moment a single urban county holds a hundred times the features of its neighbors.

When Boundary Partitions Beat a Grid #

A uniform grid — the H3, S2, or quadkey schemes — optimizes for spatial-proximity queries and even cell sizes. Administrative partitioning optimizes for a different axis: it wins when your queries, your retention obligations, and your access controls all follow political boundaries rather than geographic windows. Three signals say “partition by boundary.” First, queries name jurisdictions (“all parcels in Harris County”) far more often than they draw bounding boxes. Second, retention and legal-hold rules differ by jurisdiction — a state-mandated seven-year hold applies to one state’s records and not another’s, and you want that boundary to be a clean prefix your retention policy framework can target. Third, data ownership and deletion requests arrive per jurisdiction, so a boundary-aligned prefix makes a lawful deletion a single-prefix operation instead of a full-archive scan.

The cost of that alignment is size skew: administrative units are not equal-area or equal-density. A grid gives you controllable cell sizes; boundaries give you meaningful prefixes but leave one county holding gigabytes and another holding kilobytes. The fix is a hybrid tree that partitions by boundary down to a level that stays meaningful, then falls back to a spatial grid inside any partition that exceeds a size threshold.

Skew is not a nuisance you can average away, because it is bimodal in the way that hurts most. A national parcel or address archive typically has thousands of low-population counties that each compress to a few megabytes — well below any cold tier’s minimum object size, so each one wastes the minimum-billable footprint and adds a LIST/GET line item — alongside a handful of metropolitan counties that each hold enough features to make a single-file read pull hundreds of megabytes out of cold storage under a restore SLA. A pure boundary layout is therefore wrong at both ends of the distribution simultaneously: too many tiny objects and a few oversized ones. The hybrid approach fixes the oversized end by subdividing, and you address the tiny end by choosing how deep the administrative nesting goes — stopping at state rather than county for sparse regions, or coalescing small adjacent counties under a single state-level file. The threshold-driven procedure below encodes exactly that logic.

Step 1: Join Features to Administrative Units #

Assign every feature its jurisdiction codes with a spatial join against an authoritative boundary set — census TIGER, GADM, or a national cadastre. Use a stable code (FIPS, ISO 3166-2) rather than a display name, because names change spelling and codes do not. Perform the join in a projected CRS appropriate to the region so the point-in-polygon test is geometrically sound, then carry the codes as plain columns.

import geopandas as gpd

parcels = gpd.read_file("s3://spatial-archive/vector/2024/parcels_raw.fgb")
counties = gpd.read_file("s3://spatial-archive/reference/tiger_2023_counties.fgb")

# Project both to an equal-area CRS for a correct point-in-polygon assignment.
parcels_m = parcels.to_crs("EPSG:5070")      # CONUS Albers equal-area
counties_m = counties.to_crs("EPSG:5070")

joined = gpd.sjoin(
    parcels_m,
    counties_m[["STATEFP", "COUNTYFP", "geometry"]],
    how="left",
    predicate="within",
)

# Features that miss every polygon (offshore, digitizing slivers) must not be dropped.
unmatched = joined["COUNTYFP"].isna().sum()
print(f"unmatched features routed to _unassigned/: {unmatched}")
joined["COUNTYFP"] = joined["COUNTYFP"].fillna("_unassigned")
joined["STATEFP"] = joined["STATEFP"].fillna("_unassigned")

Step 2: Measure Skew and Set a Split Threshold #

Before writing the tree, size each candidate leaf partition. The goal is a target file-size band — large enough to clear the cold tier’s minimum object size, small enough that a single-jurisdiction query does not drag a multi-gigabyte file out of cold storage. Compute per-county byte estimates and flag the ones that need grid subdivision.

duckdb -c "
SELECT STATEFP, COUNTYFP,
       count(*)                       AS features,
       round(count(*) * 512.0 / 1e6, 1) AS est_mb   -- ~512 bytes/feature after ZSTD
FROM read_parquet('joined/*.parquet')
GROUP BY STATEFP, COUNTYFP
ORDER BY features DESC
LIMIT 5"

Any county whose est_mb exceeds your upper band (for example 512 MB) gets a second-level H3 grid; everything else is written as a single file per county.

Step 3: Write the Hybrid Tree #

Partition by STATEFP and COUNTYFP for normal jurisdictions, and add an h3_cell sub-level only inside the oversized ones. The administrative prefix survives either way, so a jurisdiction-scoped lifecycle rule or deletion still targets one clean path.

import h3
import pyarrow as pa
import pyarrow.parquet as pq

SPLIT_MB, H3_RES = 512, 7

for (statefp, countyfp), grp in joined.groupby(["STATEFP", "COUNTYFP"]):
    base = f"s3://spatial-archive/parcels/country=US/state={statefp}/county={countyfp}"
    est_mb = len(grp) * 512 / 1e6
    if est_mb > SPLIT_MB:
        pts = grp.geometry.to_crs(4326).representative_point()
        grp = grp.assign(h3_cell=[h3.latlng_to_cell(p.y, p.x, H3_RES) for p in pts])
        for cell, sub in grp.groupby("h3_cell"):
            _write(sub, f"{base}/h3_cell={cell}/part-0000.parquet")
    else:
        _write(grp, f"{base}/part-0000.parquet")

def _write(frame, path):
    table = pa.Table.from_pandas(frame.drop(columns="geometry").assign(
        geometry=frame.geometry.to_wkb()), preserve_index=False)
    pq.write_table(table, path, compression="zstd", compression_level=9,
                   row_group_size=100_000)

Compression level here leans aggressive because these are cold, rarely-read partitions; the ratio-versus-CPU trade is covered under tuning ZSTD compression for GeoParquet archives.

Validation: Confirm Balanced, Prunable Partitions #

Verify two properties: no partition file blew past the band, and a jurisdiction query prunes to a single prefix. Check the size distribution and run an EXPLAIN to confirm partition pruning fires.

# 1. No leaf file should exceed the upper band after the hybrid split.
duckdb -c "
SELECT max(file_size_mb), avg(file_size_mb), count(*)
FROM (SELECT filename, sum(bytes)/1e6 AS file_size_mb
      FROM parquet_file_metadata('s3://spatial-archive/parcels/**/*.parquet')
      GROUP BY filename)"

Expected output — the maximum stays within the band and the mean sits comfortably above the tier minimum:

┌───────────────────┬───────────────────┬──────────────┐
│ max(file_size_mb) │ avg(file_size_mb) │ count_star() │
│      double       │      double       │    int64     │
├───────────────────┼───────────────────┼──────────────┤
│       498.2       │       143.6       │     3187     │
└───────────────────┴───────────────────┴──────────────┘

A query filtered on state='06' AND county='037' should read only that county’s files. If EXPLAIN shows a full-tree scan, the partition columns are not in the path and pruning cannot fire — the single most common boundary-partitioning defect.

Boundary Partitioning Failure Modes #

Operational Execution Checklist #

• Spatially join features to an authoritative boundary set in an equal-area CRS, keying on stable codes (FIPS, ISO 3166-2) not names.
• Route unmatched features to an explicit _unassigned/ prefix so nothing is silently dropped during the join.
• Estimate per-jurisdiction partition sizes and set a split threshold that clears the cold tier minimum object size.
• Subdivide only oversized jurisdictions with a second-level H3 grid, preserving the administrative prefix on every leaf.
• Write state=/county= (and h3_cell= where split) as real Hive path segments so jurisdiction queries prune to one prefix.
• Verify the max leaf-file size stays within the target band and that EXPLAIN confirms partition pruning on a jurisdiction filter.
• Pin the boundary-set vintage in the dataset contract and version any re-vintage as a deliberate rewrite.

Related #

• Up: Spatial Partitioning Techniques — the parent reference for choosing between boundary, grid, and hybrid partition layouts.
• Choosing H3 vs S2 vs Quadkey for Archive Partitioning — the sibling grid-based approach and the H3 grid used for subdividing dense jurisdictions here.
• Implementing Lifecycle Rules for Shapefile Archives — jurisdiction-scoped retention rules that map cleanly onto the administrative prefix.
• GeoParquet Migration Workflows — the conversion pipeline that produces the columnar files this tree organizes.
• Spatial Archive Cost Modeling — model how partition count and per-jurisdiction file sizes drive storage and request costs.

Part of the Spatial Data Archival knowledge base.