#!/usr/bin/env python3 """Cross-check block building heights against Mapbox's building data. Mapbox publishes building footprints with ``height`` / ``min_height`` / ``type`` attributes in the ``mapbox-streets-v8`` tileset. They come from different source data than our OSM extract, so comparing the two is a cheap way to find buildings where our height is wrong or missing — the ``0.1.2`` finding that three "guessed" volumes were actually a roof, a roof and a bus shelter is exactly the class of problem this tool hunts for. Usage ----- # token from the environment or the ignored .env file python3 Geo/tools/mapbox_buildings.py \ --osm Geo/raw/osm-api-map-2026-11-06.osm \ --out Geo/verification/2026-11-09-mapbox-buildings.md # show only disagreements and unverifiable buildings python3 Geo/tools/mapbox_buildings.py --only-diffs Every Tilequery response is cached under ``Geo/raw/mapbox-cache/`` so a re-run costs nothing and the report is reproducible. The cache is gitignored along with ``Geo/raw/``. Design notes ------------ * Standard library only (``urllib``), like the rest of ``Geo/tools``. * One request per building, spaced by ``--delay`` (default 1.2 s) — Mapbox bills per request and rate-limits bursts. * The report is a *measurement*: it prints both figures and the delta, never a single verdict, so a human can decide which source to trust for a given building. * A Tilequery hit is only accepted as "this building" when the matched polygon contains or nearly contains the queried point (``distance`` ≤ ``--max-distance``). Licensing / data handling ------------------------- Mapbox terms restrict persisting and redistributing their data, so **neither the cache nor the generated report is committed** (`.gitignore` covers `Geo/raw/mapbox-cache/` and `Geo/verification/*mapbox*`). They are local working material: re-run the tool rather than sharing the output. The token is read from the environment or the ignored `.env` and is never written to the cache, the report or the JSON. """ from __future__ import annotations import argparse import json import math import os import sys import time import urllib.error import urllib.parse import urllib.request sys.path.insert(0, os.path.dirname(os.path.abspath(__file__))) import block_height_report as bhr # noqa: E402 import osm_common as oc # noqa: E402 ROOT = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..")) CACHE_DIR = os.path.join(ROOT, "Geo", "raw", "mapbox-cache") TILESET = "mapbox.mapbox-streets-v8" LAYER = "building" API = "https://api.mapbox.com/v4" def load_token(explicit: str | None = None) -> str: """Token from --token, the environment, or the ignored .env file.""" if explicit: return explicit.strip() token = os.environ.get("MAPBOX_ACCESS_TOKEN") if token: return token.strip() env_path = os.path.join(ROOT, ".env") if os.path.exists(env_path): with open(env_path, encoding="utf-8") as handle: for line in handle: line = line.strip() if line.startswith("MAPBOX_ACCESS_TOKEN="): return line.split("=", 1)[1].strip().strip('"').strip("'") raise SystemExit( "no Mapbox token: pass --token, set MAPBOX_ACCESS_TOKEN, or add it to .env") def centroid_lat_lon(points: list[tuple[float, float]]) -> tuple[float, float]: """Centroid of a ring of (lat, lon) points, as (lat, lon).""" lat = sum(p[0] for p in points) / len(points) lon = sum(p[1] for p in points) / len(points) return lat, lon def cache_path(lat: float, lon: float, radius: int) -> str: return os.path.join(CACHE_DIR, f"{lat:.6f}_{lon:.6f}_r{radius}.json") def tilequery(lat: float, lon: float, token: str, radius: int, delay: float, use_cache: bool = True) -> tuple[dict, bool]: """Tilequery one point. Returns (payload, from_cache).""" path = cache_path(lat, lon, radius) if use_cache and os.path.exists(path): with open(path, encoding="utf-8") as handle: return json.load(handle), True query = urllib.parse.urlencode({ "radius": radius, "layers": LAYER, "limit": 10, "access_token": token, }) url = f"{API}/{TILESET}/tilequery/{lon:.6f},{lat:.6f}.json?{query}" request = urllib.request.Request(url, headers={"User-Agent": "gig-line-geo/1.0"}) payload = None for attempt in range(4): try: with urllib.request.urlopen(request, timeout=30) as response: payload = json.load(response) break except urllib.error.HTTPError as exc: if exc.code == 429 or exc.code >= 500: wait = delay * (2 ** attempt) sys.stderr.write(f" HTTP {exc.code}; retrying in {wait:.1f}s\n") time.sleep(wait) continue raise except urllib.error.URLError as exc: sys.stderr.write(f" network error ({exc.reason}); retrying\n") time.sleep(delay * (2 ** attempt)) if payload is None: raise SystemExit(f"tilequery failed for {lat:.6f},{lon:.6f}") if token in json.dumps(payload): raise SystemExit("refusing to cache a response that contains the token") os.makedirs(CACHE_DIR, exist_ok=True) with open(path, "w", encoding="utf-8") as handle: json.dump(payload, handle) time.sleep(delay) # only sleep on a real request, never on a cache hit return payload, False def nearest_building(payload: dict, max_distance: float) -> dict | None: """The closest building feature within `max_distance` metres, or None.""" best = None for feature in payload.get("features", []): props = feature.get("properties", {}) if props.get("tilequery", {}).get("layer") != LAYER: continue distance = props.get("tilequery", {}).get("distance") if distance is None or distance > max_distance: continue if best is None or distance < best[0]: best = (distance, props) return None if best is None else best[1] def compare(ours: float, ours_source: str, theirs: float | None) -> tuple[str, float | None]: """Classify the agreement between our height and Mapbox's. Returns (verdict, delta_m). Verdicts are deliberately blunt so the report can be scanned: "match", "close", "disagree", "mapbox-only", "ours-only". """ ours_is_volume = not ours_source.startswith("not a building") if not ours_is_volume: return "non-building", None if theirs is None: return "ours-only", None if ours <= 0.0: return "mapbox-only", theirs delta = theirs - ours if abs(delta) <= 0.5: return "match", delta if abs(delta) <= 2.0: return "close", delta return "disagree", delta def build_rows(nodes, ways, token, radius, max_distance, delay) -> list[dict]: rows = [] for wid, tags, refs in ways: if "building" not in tags: continue points = oc.way_points(nodes, refs) if not points: continue lat, lon = centroid_lat_lon(points) payload, cached = tilequery(lat, lon, token, radius, delay) feature = nearest_building(payload, max_distance) ours_height, ours_source = bhr.classify_height(tags, way_id=wid) theirs = None if feature is None else feature.get("height") verdict, delta = compare(ours_height, ours_source, theirs) rows.append({ "way": wid, "name": tags.get("name") or tags.get("addr:housenumber") or "", "kind": tags.get("building"), "osm_height_m": round(ours_height, 2) if ours_height else 0.0, "osm_source": ours_source, "mapbox_height_m": None if theirs is None else round(float(theirs), 2), "mapbox_min_height_m": None if feature is None else feature.get("min_height"), "mapbox_type": None if feature is None else feature.get("type"), "mapbox_distance_m": None if feature is None else round(feature.get("tilequery", {}).get("distance", 0.0), 1), "verdict": verdict, "delta_m": None if delta is None else round(delta, 2), "cached": cached, "lat": round(lat, 6), "lon": round(lon, 6), }) return rows def render_markdown(rows: list[dict], args) -> str: order = {"disagree": 0, "mapbox-only": 1, "ours-only": 2, "non-building": 3, "close": 4, "match": 5} rows = sorted(rows, key=lambda r: (order.get(r["verdict"], 9), -abs(r["delta_m"] or 0))) counts = {} for row in rows: counts[row["verdict"]] = counts.get(row["verdict"], 0) + 1 lines = [ "# Block buildings — OSM vs Mapbox", "", "Generated by `Geo/tools/mapbox_buildings.py`; re-run it to reproduce.", "Mapbox data: tileset `mapbox-streets-v8`, layer `building`, Tilequery", f"at each building centroid (radius {args.radius} m, accept ≤ " f"{args.max_distance} m). OSM heights come from the tagged extract, via the", "same classification used for the blockout (tag > levels > guess; roof and", "shelter are not building volumes).", "", f"Verdicts: {', '.join(f'{k} {v}' for k, v in sorted(counts.items()))}.", "", "| Verdict | Δ | Building | OSM | OSM source | Mapbox | Mapbox type | dist |", "| --- | ---: | --- | ---: | --- | ---: | --- | ---: |", ] for row in rows: osm = f"{row['osm_height_m']:.1f} m" mb = "—" if row["mapbox_height_m"] is None else f"{row['mapbox_height_m']:.1f} m" delta = "—" if row["delta_m"] is None else f"{row['delta_m']:+.1f} m" dist = "—" if row["mapbox_distance_m"] is None else f"{row['mapbox_distance_m']:.0f} m" name = row["name"] or f"way {row['way']}" lines.append( f"| {row['verdict']} | {delta} | {name} | {osm} | {row['osm_source']} | " f"{mb} | {row['mapbox_type'] or '—'} | {dist} |") lines += [ "", "## How to read this", "", "- **disagree** — the two sources differ by more than 2 m. One of them is", " wrong; check the reference photography before changing the model.", "- **mapbox-only** — Mapbox knows a height we do not (usually a building", " with no OSM `height`/`building:levels`).", "- **ours-only** — Mapbox returns no building within the accept distance;", " the footprint may be missing from its tileset, or the centroid may fall", " outside its polygon.", "- **non-building** — our classification says roof/shelter; excluded from", " height comparison on purpose.", "", "Neither source is authoritative: OSM is community-tagged, Mapbox is a", "derived commercial dataset. Treat a disagreement as a reason to look at", "the photographs, not as an automatic correction.", "", ] if args.only_diffs: lines.insert(0, "> Filtered: only disagreements and unverifiable buildings.\n") return "\n".join(lines) def main() -> int: ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter) ap.add_argument("--osm", default="Geo/raw/osm-api-map-2026-11-06.osm") ap.add_argument("--out", help="Markdown report path; omit to print to stdout") ap.add_argument("--json", help="Also write the merged rows as JSON here") ap.add_argument("--token", help="Mapbox access token (default: env/.env)") ap.add_argument("--radius", type=int, default=25, help="Tilequery radius, metres") ap.add_argument("--max-distance", type=float, default=8.0, help="accept a hit only within this distance of the centroid") ap.add_argument("--delay", type=float, default=1.2, help="seconds between requests") ap.add_argument("--only-diffs", action="store_true", help="report only disagreements / unverifiable buildings") ap.add_argument("--no-cache", action="store_true", help="ignore the tilequery cache") args = ap.parse_args() token = load_token(args.token) nodes, ways = oc.load(args.osm) rows = build_rows(nodes, ways, token, args.radius, args.max_distance, args.delay) if args.no_cache: for row in rows: row["cached"] = False if args.only_diffs: keep = {"disagree", "mapbox-only", "ours-only", "non-building"} rows = [r for r in rows if r["verdict"] in keep] markdown = render_markdown(rows, args) if args.out: os.makedirs(os.path.dirname(os.path.abspath(args.out)), exist_ok=True) with open(args.out, "w", encoding="utf-8") as handle: handle.write(markdown + "\n") print(f"wrote {args.out} ({len(rows)} buildings)") else: print(markdown) if args.json: payload = { "generated_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()), "osm": args.osm, "mapbox": {"tileset": TILESET, "layer": LAYER, "radius_m": args.radius, "max_distance_m": args.max_distance}, "buildings": rows, } with open(args.json, "w", encoding="utf-8") as handle: json.dump(payload, handle, indent=2) print(f"wrote {args.json}") return 0 if __name__ == "__main__": raise SystemExit(main())