isaaccorley/geospatial-skills/parallel-geospatial

Parallel Geospatial

Reference for designing and debugging large-scale parallel geospatial processing on a single fat node. You do not write geospatial code at this scale single-threaded. The defaults below come from real continental-scale (15M polygons, 6T-pixel rasters) failures.

The shape of every problem here is the same: partition along an axis where work is naturally independent, then run an embarrassingly parallel pool over the partitions. The skill is identifying that axis cheaply and avoiding the operation-specific footguns that crash GEOS or DuckDB at scale.

The decision tree (start here)

For a new geospatial operation at >1M-feature or >1B-pixel scale, ask in this order:

1. Is there a natural disjoint partition? state, county, tile, hilbert-curve cell, raster window, time slice. If yes → ProcessPool over partitions. Stop, you're done thinking.
2. Is the work output-disjoint but input-overlapping? (e.g. spatial join: each output row needs a small bbox query against an input). Build an STRtree once on inputs, partition over outputs. Workers share a read-only tree via fork-inherited memory or rebuild small per worker.
3. Is it a global aggregation that can't be partitioned? Decompose into a tree-reduce: partition → reduce-per-partition → small final merge. Almost every real "global" op decomposes this way (union, dissolve, count-distinct, hilbert sort).
4. Is it truly serial? Rare. Consider whether you actually need the operation, or whether a different formulation (e.g. raster mask AND vs. polygon difference) sidesteps it.

If you find yourself running anything covering a continent in a single thread, stop and reconsider — single-threaded for >1M features is almost always a design bug, not a "we'll optimize later" one.

Footguns that cost whole batch jobs (real failures)

These are operations whose textbook one-liner segfaults or silently corrupts at continental scale. Bake the workaround in from the start; do not "try the simple version first."

When in doubt about whether a given DuckDB-spatial or GEOS operation will scale, assume it will segfault. Cost of bake-in workaround: 30 min. Cost of discovering the segfault on a real run: a 4h job, wasted allocation, plus the time to find the bug after the empty stderr.

Embarrassingly-parallel patterns

Pattern A — Per-tile / per-state ProcessPool (the default)

The 80% case. Used for: raster polygonize, postprocess spatial joins, tile-edge pair finding, multi-year raster combines.

import os
from concurrent.futures import ProcessPoolExecutor
def worker_count(frac=0.95):
    return max(1, int((os.cpu_count() or 1) * frac))
def parallel_map(fn, items, max_workers):
    with ProcessPoolExecutor(max_workers=max_workers) as pool:
        return list(pool.map(fn, items))
partitions = partition_by_state(national)  # or by tile, region, etc.
results = parallel_map(
    _worker_fn,
    partitions,
    max_workers=worker_count(0.95),  # leave headroom for the OS / monitoring
)
merged = combine(results)

The _worker_fn should:

• Take ONLY picklable args (Arrow buffers, numpy arrays, paths, primitives — not DuckDB connections, not open file handles).
• Open its own DuckDB connection / Rasterio handle / etc. inside the worker. Workers share forked memory but you cannot share an open driver state.
• Return small results. If the result is large (e.g. 1GB of polygons), write to disk in the worker and return only the path.

Pattern B — Tree-reduce (when "global" aggregation is needed)

Example — a national tile-edge dissolve:

15M polygons globally
   │
   ├── partition by tile_id      # (a) cheap: just a column already on the row
   │
   ├── per-tile envelope         # (b) trivial groupby in DuckDB
   │
   ├── on_edge filter            # (c) keep only polygons whose bbox touches its tile envelope
   │       ↓
   │   ~500k polygons remain     # ~3% of original — STRtree now actually fits
   │
   ├── partition by state        # (d) tile-edge dissolves never cross state lines
   │
   ├── ProcessPool: STRtree+touches per state    # (e) pair finding
   │
   ├── union-find on combined pairs              # (f) connected components, fast in pure Python
   │
   └── ProcessPool: shapely.union_all per group  # (g) most groups have 2-3 members

Lesson: filter aggressively before the expensive step. The on_edge filter (c) drops 97% of work. Without it, even the parallel STRtree is too slow.

Pattern C — Output-disjoint, input-overlapping spatial join

Example — assigning each polygon its largest-overlap administrative boundary. STRtree of boundaries (small, ~3000 counties), partition over polygons:

boundary_tree = STRtree(boundary_geoms)
# In each worker:
def assign_county(polygon_chunk):
    out = []
    for poly in polygon_chunk:
        candidates = boundary_tree.query(poly, predicate="intersects")
        # Pick largest-overlap candidate by ST_Intersection area.
        ...
    return out

Workers share boundary_tree via fork (Linux). On Windows / macOS spawn, rebuild per worker (cheap if input is small).

Always-applicable defaults

These are the "no thinking required" choices for a new geospatial pipeline. Deviate only with a reason.

• CRS: pick one and lock it (EPSG:5070 for CONUS, EPSG:6933 for global equal-area). Never reproject silently mid-pipeline.
• GeoParquet: 1.1 with full-PROJJSON CRS. The short {id: {authority, code}} form breaks pyproj 3.x, geopandas, pyogrio, and GDAL readers.
• Bbox columns: explicit xmin/ymin/xmax/ymax in addition to geometry. Lets DuckDB row-group statistics prune queries before any GEOS work runs.
• Hilbert sort before write for any parquet that will be range-queried by bbox: ORDER BY ST_Hilbert(geom, BOX_2D). Co-locates spatially-near features in the same row group.
• Row group size: 50k–100k features. Smaller groups = better pruning, larger groups = better compression. 50k is a good default for ~30M-feature national datasets.
• Compression: zstd always. Never default snappy for vector — zstd is 30-50% smaller for similar speed.
• Workers: worker_count(0.95) of physical cores. The 5% headroom matters on shared nodes.
• dtype: int64 for IDs and labels. Cost of int64 vs int32: nothing. Cost of int32 overflow: silent data loss.
• Validation gate at every boundary: geom.notna() & ~geom.is_empty & geom.is_valid before any write that builds a spatial index (FGB, STRtree, DuckDB spatial parquet).

DuckDB-spatial-specific guidance

DuckDB-spatial is the right tool for a lot of this work but it has sharp edges. Treat these as load-bearing.

• Heavy GEOS ops segfault on big tables. ST_Union_Agg, ST_Touches self-joins, ST_Buffer over transformed lines. Materialise to Arrow + use shapely / numpy when you have >1M rows.
• Bbox stats are how you get pruning. Read parquet with explicit xmin/xmax/ymin/ymax; queries with WHERE xmax >= ? AND xmin <= ? ... will skip whole row groups before decoding any geometry.
• `always_xy => true` on every `ST_Transform` from `EPSG:4326`. Otherwise lat/lon swap → Infinity coords.
• Don't trust DuckDB's binder on wide self-joins with spatial predicates. Even when it doesn't crash, it produces "Failed to bind column reference" errors. Workaround: split into geometry-only spatial join CTE then attribute-join.
• `GROUP BY ALL` is great alone, broken with explicit keys. Use GROUP BY ALL by itself, or list every key explicitly. Do not mix.
• Use `BIND` over the parquet file, not over a TABLE you SELECTed into. read_parquet('path') lets DuckDB use file statistics; an in-memory CREATE TABLE ... AS SELECT ... discards them.

Shapely-specific guidance

• Vectorise everything. shapely.from_wkb(list) and shapely.buffer(geoms, m) are 50–100× faster than per-geometry loops. Never write [g.buffer(m) for g in geoms] over more than 1000 geoms.
• Always `make_valid` before union/intersect. Cheap, prevents 100% of TopologyException crashes.
• `shapely.coverage_union_all` for non-overlapping polygons; `union_all` for arbitrary. Coverage is ~10× faster but only correct on coverages (no overlaps, no gaps).
• STRtree is the only correct way to do bbox prefiltering at scale. Don't write your own grid index — shapely.strtree.STRtree is GEOS-backed and bulk-vectorized.
• STRtree.query returns a (2, N) array of (input_idx, tree_idx) pairs. When self-querying, filter a < b to dedupe. Predicate values: 'intersects', 'touches', 'contains', 'within'.

Shared clusters and batch schedulers