Pmtiles Pipeline

PMTiles Pipeline

Reference for building a single global PMTiles archive on one fat node, distilled from a real global change-layer build (~6T-px raster input, 176M polygons, 28 GB PMTiles output). The pipeline is:

Copy code
raster.zarr  →  geoparquet shards  →  FlatGeobuf batches  →  sharded tippecanoe  →  tile-join  →  pmtiles
   (1)              (2)                    (3)                     (4)                  (5)

Each stage's "why" matters more than the syntax. Use this skill whenever a step in this pipeline is being built, debugged, or optimized.

Stage 1 — Raster → polygons (as parquet shards)

• Use contourrs.shapes_arrow to vectorize a classified uint8 raster directly to Arrow (skip GeoJSON in/out).
• Pre-filter with scipy.ndimage.binary_opening(struct=NxN) per class to wipe salt-and-pepper noise. N=3 is a sane default at 10m resolution.
• Write GeoParquet via geopandas with compression="zstd". One parquet per super-tile (e.g., 8192² px) → tens of thousands of small shards. That granularity matters for downstream resume + parallelism.
• Track fields: class (uint8), area_deg2 (float64 — useful for downstream filtering), geometry.
• Connectivity 8 (diagonals count) and min-pixels=25 are fine starting defaults; higher min-pixels = aggressive smoothing.

Stage 2 — Field-size analysis before any pmtiles work

Always do this once. It informs the area filter at stage 3 and the simplification tolerance.

Sample ~100 random parquet shards, read class + area_deg2 only (no geometry — pq.read_table(p, columns=[...]).to_pandas(), geopandas refuses to read parquet without geom). Convert to m² via area_deg2 * 111320**2 (equator approximation, slight overestimate at high latitudes).

Look at percentiles + cumulative-below thresholds at: 0.25, 0.5, 1, 2.5, 5, 10 ha. The min observed area = polygonize floor + small expansion from simplify. Anything in the bottom decile is artifact noise; smallholder ag fields are typically 0.5–2 ha (median ~0.7 ha real-world).

Stage 3 — GeoParquet → FlatGeobuf (parallel, with filter + simplify)

Why FlatGeobuf: tippecanoe parallelizes natively on FGB input via the spatial index. No --read-parallel flag needed. GeoJSON-seq input forces single-thread JSON parsing — kills 96-core scaling. Tippecanoe does NOT read GeoParquet directly — FGB is the only fast format.

Implementation pattern:

• Round-robin shards into N batches (default N = physical core count) so every batch covers global geography.
• ProcessPool worker per batch. Each worker:
  1. Concat its parquet shards into one GeoDataFrame
  2. Filter by area_deg2 >= MIN (default 3.2e-7 deg² ≈ 0.4 ha drops ~15%)
  3. shapely.simplify(geoms, tolerance, preserve_topology=True) (default 2e-4 deg ≈ 22m at equator → smooths contour-step edges)
  4. Drop NULL/empty/invalid: geom.notna() & ~geom.is_empty & geom.is_valid. This is mandatory — pyogrio's FGB writer with spatial index errors out on NULL geometry mid-write. Without this filter you'll get FeatureError: NULL geometry not supported with spatial index.
  5. pyogrio.write_dataframe(gdf, path, driver="FlatGeobuf") — Arrow path, fast.
• Expected: 36 GB FGB from 26 GB parquet at this scale (5× expansion: spatial index + no row-level compression).
• Wall time: ~7 min on 96 cores for 124M filtered polygons.

Min-area threshold rule of thumb: at maxzoom=14 (9.6 m/px), a 50m feature = 5×5 px = the practical "easily visible" floor. Anything <2500 m² (25 px²) renders as 1–2 px even at maxzoom = noise. Default 0.4 ha is conservative; bump to 1 ha for cleaner output at the cost of dropping ~64% of polys (mostly smallholder fields).

Stage 4 — Build pmtiles via sharded tippecanoe

Single-process tippecanoe does NOT scale to 96 cores at 100M+ feature scale. Even with FGB native parallelism, tippecanoe's sort + final-write phases serialize to ~10–20 effective cores. A single-process build of a 124M-poly dataset timed out at 8h on a 96-core node. The sharded path finished in ~50 min.

Architecture

Copy code
1. Stage N FGBs from network storage (NFS) to local NVMe scratch   ← critical for I/O
2. Run N tippecanoes in parallel, TIPPECANOE_MAX_THREADS=1, one per FGB
   → N sub-pmtiles on local scratch
3. tile-join all N sub-pmtiles into one final pmtiles on local scratch
4. cp final pmtiles back to network storage

Step 1 is the biggest non-obvious win. N concurrent readers hit NFS hard when reading network storage directly (measured ~70 MB/s aggregate ceiling). NVMe local scratch gives 1–2 GB/s aggregate read. Always stage at this scale.

tippecanoe flags (zoom decisions are dataset-specific, not universal)

These flags are tuned for a 10 m agricultural change layer (polygons). Adjust for your dataset.

bashCopy code
-l <layer>                     # layer name in MVT
--minimum-zoom=3               # global view from this pmtiles. If you serve a separate raster overview for low zooms (a common split: pmtiles for z>=11 + COG for z<11), set higher.
--maximum-zoom=14              # 9.6 m/px ≈ 10m source native. Don't go past native unless you want client-side overzoom only.
--base-zoom=14                 # full detail held until maxzoom — no premature simplify. Use lower (e.g. 12) if you want smaller tiles at intermediate zooms.
--full-detail=14               # max vertex precision per tile (default 12). Bumping helps small-feature visibility at maxzoom.
--simplification=1             # near-zero DP simplify. Polygonize already simplified at stage 3; let high-zoom be detailed.
--buffer=8                     # pixel buffer on tile edges to prevent feature clipping mid-render. Default 5 sometimes shows seams.
--drop-rate=1                  # ~no rate-based dropping (default 2.5). Bigger pmtiles, more features at low zooms — set to 2–3 if you want a smaller archive.
--maximum-tile-bytes=2000000   # 2MB cap (default 500K). Web renderers handle 2MB fine; lower if mobile-perf matters.
--maximum-tile-features=1000000# 1M cap (default 200K). Same trade-off.
-y class                       # carry only the `class` attribute; drop other parquet columns to save bytes.
--coalesce-smallest-as-needed  # merge sub-pixel adjacent same-class polys at low zoom (preserves info vs dropping)
--drop-densest-as-needed       # safety net: tippecanoe ERRORS without this if a tile blows past the byte cap. Always include.
--detect-shared-borders        # geometry compression for adjacent same-class polys (5-10% size win, slight CPU cost). Drop in sharded mode — round-robin splits adjacent features across shards so the flag does little.

DO NOT use --no-feature-limit or --no-tile-size-limit. Both are present in many starter scripts; they produce tiles browsers can't render. Tippecanoe's defaults are correct.

Sharded caveat: tile-join (Felt v2.79) doesn't support --maximum-tile-bytes / --drop-densest-as-needed at merge time. Per-shard caps are the only enforcement. Round-robin sharding makes this OK — merged tile sizes recover to ~natural global density.

What didn't work + why (don't repeat these mistakes)

Stage 5 — Faster path for repeated runs: geographic-sharding + pmtiles merge

If you'll iterate on this pipeline regularly, the real speedup is here. The architecture worth investing in:

Copy code
1. Geographic shard at z=3 tile boundaries → 64 cells
2. Build z=0..2 once from all FGBs combined (small, ~1-2 min)
3. Build z=3..14 in 64 parallel tippecanoes, each over its own cell's FGB
   → 64 sub-pmtiles with DISJOINT (non-overlapping) tile coverage
4. go-pmtiles `pmtiles merge low.pmtiles cell_*.pmtiles → final.pmtiles`
   → byte-level concatenation, no MVT decode/encode, ~2-5 min

Why dramatically faster: pmtiles merge is O(bytes copied), tile-join is O(features × overlapping_inputs). For 124M features the difference is 5-10× on the merge step alone.

Caveats:

• Polygon clipping at cell boundaries needed (or accept seam artifacts at high zoom)
• z<3 must use feature-level merging (separate small build)
• All inputs must be clustered per pmtiles spec (tippecanoe writes clustered by default)

Decision rule: one-shot build → tile-join is fine. Annual / repeated builds → invest in geographic sharding.

Cluster and I/O notes

• At this scale, request a full fat node (96+ cores) and stage inputs to local NVMe scratch — never read 96-way concurrent from a network filesystem.
• The single-process path needs an 8h+ walltime budget (and timed out anyway); the sharded path fits comfortably in 1.5h.
• Make every stage resumable: polygonize skips existing parquet shards, the FGB step skips existing batch outputs, and wrapper scripts trap on EXIT and print the resume command.
• Never pass --force from a wrapper script. A wrapper that auto-passed --force once wiped a near-complete zarr store when re-run to fill in errored shards. Add a preflight that refuses to start if the existing store is missing expected arrays.

Quick reference: when each tool is right

Common errors + fixes

• FeatureError: NULL geometry not supported with spatial index during FGB write — drop geom.notna() & ~geom.is_empty & geom.is_valid rows before write.
• tile-join: unrecognized option '--maximum-tile-bytes' — Felt tile-join doesn't expose this; per-shard caps from sub-builds are the only enforcement.
• Tippecanoe progress hangs at "99.9% 14/x/y" — stuck on the very densest z=14 tiles in agricultural regions (Iowa/Punjab/Po Valley). Either wait, or drop --detect-shared-borders and rebuild.
• pmtiles file size grows then stops, but job is still running — single-threaded final-seal phase. Can take 30+ min on a multi-GB output written to network storage. Stage output to local scratch first to mitigate.