gig-line/Geo/tools/mapbox_buildings.py

#!/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())