Building a zoning vector tile pyramid with tippecanoe

# Measure the simplification error before choosing a value.
tippecanoe -z16 -Z10 -o /tmp/probe.pmtiles --simplification=1 parcels.geojson
python3 - <<'PY'
# Decode one tile and compare against the source geometry, per zoom.
from mvt_probe import max_symmetric_difference
for z in range(10, 17):
    err = max_symmetric_difference("/tmp/probe.pmtiles", "parcels.geojson", zoom=z)
    print(f"z{z}: worst-case symmetric difference {err:.2f} m²")
PY

NOISE_FLOOR_M2 = 2.0          # from the change detector; the tiles must stay under it

# Tolerance in metres, per zoom. Chosen so the resulting area error stays below the
# floor for a typical urban parcel perimeter (~120 m): error ≈ perimeter × tolerance.
TOLERANCE_M = {
    10: 8.0,    # regional: parcels are dissolved into districts anyway
    11: 4.0,
    12: 2.0,
    13: 1.0,
    14: 0.4,
    15: 0.15,
    16: 0.0,    # full geometry: this is the zoom people measure at
}


def tolerance_is_safe(zoom, typical_perimeter_m=120.0) -> bool:
    return TOLERANCE_M[zoom] * typical_perimeter_m <= NOISE_FLOOR_M2 or zoom <= 12

#!/usr/bin/env bash
set -euo pipefail

SRC=parcels.geojson          # EPSG:4326, one feature per parcel, id as a STRING
OUT=tiles
VERSION="$(date -u +%Y%m%dT%H%M%SZ)"

# Districts only, dissolved, for the regional zooms.
tippecanoe \
  --output="${OUT}/districts-${VERSION}.pmtiles" \
  --layer=districts \
  --minimum-zoom=8 --maximum-zoom=12 \
  --coalesce --reorder --detect-shared-borders \
  --simplification=4 \
  --no-tile-size-limit \
  districts.geojson

# Parcels, full detail at the top zoom, tapering below it.
tippecanoe \
  --output="${OUT}/parcels-${VERSION}.pmtiles" \
  --layer=parcels \
  --minimum-zoom=13 --maximum-zoom=16 \
  --use-attribute-for-id=feature_id \
  --preserve-input-order \
  --detect-shared-borders \
  --simplification=1 \
  --no-simplification-of-shared-nodes \
  --no-feature-limit --no-tile-size-limit \
  --no-tiny-polygon-reduction \
  "${SRC}"

import hashlib


def feature_id(parcel_id: str) -> int:
    """The MVT spec's id is an unsigned integer while parcel ids are strings, so a
    hash is required. 53 bits keeps it safe in JavaScript, and collisions must be
    checked across the WHOLE layer — a collision silently merges two parcels."""
    digest = hashlib.blake2b(parcel_id.encode(), digest_size=8).digest()
    return int.from_bytes(digest, "big") >> 11          # 53 bits


def prepare_source(parcels):
    """Emit GeoJSON with a numeric id AND the original string kept as an attribute."""
    seen = {}
    for p in parcels:
        fid = feature_id(p.parcel_id)
        if fid in seen and seen[fid] != p.parcel_id:
            raise ValueError(f"feature id collision: {p.parcel_id} and {seen[fid]} "
                             f"both hash to {fid}")
        seen[fid] = p.parcel_id
        yield {
            "type": "Feature",
            "id": fid,
            "properties": {
                "feature_id": fid,
                "parcel_id": p.parcel_id,        # the string, unmangled
                "zoning_class": p.canonical_class,
                "overlays": ",".join(sorted(p.overlay_codes)),
            },
            "geometry": p.geometry_geojson,
        }

import json


def write_manifest(out_dir, version, source_digest, effective_date, taxonomy_version):
    """Without this a screenshot of the map is an undated assertion."""
    manifest = {
        "pyramid_version": version,
        "generated_at_utc": version,             # the version IS the timestamp
        "source_snapshot_digest": source_digest,
        "as_of": effective_date.isoformat(),     # the zoning effective date, not the build time
        "taxonomy_version": taxonomy_version,    # what the style file must match
        "tolerance_m_by_zoom": TOLERANCE_M,
        "noise_floor_m2": NOISE_FLOOR_M2,
        "layers": {"districts": [8, 12], "parcels": [13, 16]},
    }
    (out_dir / f"manifest-{version}.json").write_text(json.dumps(manifest, indent=2))
    return manifest

def test_simplification_stays_under_the_noise_floor():
    for zoom in range(13, 17):                   # zooms where parcels are rendered
        err = max_symmetric_difference(TILES, SOURCE, zoom=zoom)
        assert err <= NOISE_FLOOR_M2, (
            f"z{zoom} simplification error {err:.2f} m² exceeds the {NOISE_FLOOR_M2} m² "
            f"change-detection floor — the map and the API disagree about geometry")


def test_no_parcel_is_dropped():
    for zoom in range(13, 17):
        assert count_features(TILES, "parcels", zoom) == SOURCE_FEATURE_COUNT


def test_feature_ids_are_stable_across_rebuilds():
    a = feature_ids(build(SOURCE, version="v1"))
    b = feature_ids(build(SOURCE, version="v2"))
    assert a == b


def test_parcel_id_keeps_its_leading_zero():
    props = feature_properties(TILES, feature_id("0714221"))
    assert props["parcel_id"] == "0714221"       # not 714221


def test_shared_borders_have_no_gaps():
    """Two adjacent parcels simplified independently drift apart. This is the visual
    artefact that makes a simplified parcel map look broken."""
    for zoom in (13, 14, 15):
        assert max_gap_between_neighbours(TILES, zoom) < 0.05     # 5 cm


def test_manifest_is_present_and_dates_the_data():
    m = read_manifest(TILES)
    assert m["as_of"] and m["source_snapshot_digest"] and m["taxonomy_version"]