geo: one tool directory, one unit boundary, per-asset provenance

Geo/tools/ is now every geographic script, with Geo/README.md documenting each
one; the Blender-side builders moved beside their outputs in Art/blender/*/ and
import their helpers from Geo/tools/osm_common.py.

That module gains to_unreal_local() and to_blender_local(): the only place in the
repository that multiplies metres by 100, so a second silent scale cannot be
introduced. block_height_report.py adds height provenance for the block — 3 of
15 buildings have no height data and render at a placeholder 6.4 m.
asset_metadata.py writes and audits <asset>.geo.json sidecars; four now exist
(hero-facade, terrain, street, blockout) and the audit recomputes each centre from
the origin, failing anything that disagrees.

Geo/tools/asset_metadata.py

@@ -0,0 +1,132 @@
+#!/usr/bin/env python3
+"""Read, write and audit per-asset geographic metadata — the provenance sidecar.
+
+Every exported asset gets `<asset>.geo.json` next to it carrying the boundary it
+was modelled against, its position in EPSG:26911 and in Unreal centimetres, the
+reasoning behind every non-measured dimension, and the unit conversion used. The
+file makes three things checkable after the fact:
+
+  * where the asset belongs in the world, without re-deriving it;
+  * which measurements are evidence and which are assumptions;
+  * that the metres-to-centimetres conversion was applied once, consistently.
+
+    Geo/tools/asset_metadata.py write --asset hero-facade --kind building \\
+        --bounds "OSM way 424024147" --centre-easting 376381.9 --centre-northing 3774184 \\
+        --out Art/blender/hero
+    Geo/tools/asset_metadata.py audit Art/blender/hero
+"""
+
+import argparse
+import datetime
+import glob
+import json
+import os
+import sys
+
+import osm_common as oc
+
+SCHEMA = 1
+
+
+def build(asset, kind, bounds, easting, northing, elevation, report_path, extra):
+    x, y, z = oc.to_unreal_local(easting, northing, elevation if elevation else 0.0)
+    doc = {
+        "schema": SCHEMA,
+        "asset": asset,
+        "kind": kind,
+        "generated_utc": datetime.datetime.now(datetime.timezone.utc).isoformat(timespec="seconds"),
+        "geographic": {
+            "crs": "EPSG:26911",
+            "source": bounds,
+            "centre_easting_m": easting,
+            "centre_northing_m": northing,
+            "elevation_navd88_m": elevation,
+            "origin": {
+                "easting_m": oc.ORIGIN_E,
+                "northing_m": oc.ORIGIN_N,
+                "reference": "Geo/manifests/local-origin.md",
+            },
+        },
+        "unreal": {
+            "centre_cm": [round(x, 1), round(y, 1), None if z is None else round(z, 1)],
+            "conversion": f"{oc.METRES_TO_UNREAL:g} uu per metre, applied once in osm_common.to_unreal_local",
+            "axes": "X east, Y south (FBX preset flips north), Z up",
+        },
+        "dims_report": report_path,
+        "notes": extra,
+    }
+    return doc
+
+
+def cmd_write(args):
+    doc = build(args.asset, args.kind, args.bounds, args.centre_easting, args.centre_northing,
+                args.elevation, args.report, args.note or [])
+    path = os.path.join(args.out, args.asset + ".geo.json")
+    os.makedirs(args.out, exist_ok=True)
+    with open(path, "w") as f:
+        json.dump(doc, f, indent=2)
+        f.write("\n")
+    print(f"wrote {path}")
+
+
+def cmd_audit(targets):
+    files = []
+    for t in targets:
+        if os.path.isdir(t):
+            # assets keep their sidecars either in the same folder as the .blend
+            # or one level down, so a directory target must recurse
+            files += glob.glob(os.path.join(t, "**", "*.geo.json"), recursive=True)
+        else:
+            files.append(t)
+    bad = 0
+    for path in sorted(files):
+        try:
+            doc = json.load(open(path))
+        except Exception as exc:  # noqa: BLE001
+            print(f"UNREADABLE {path}: {exc}")
+            bad += 1
+            continue
+        missing = [k for k in ("asset", "kind", "geographic", "unreal") if k not in doc]
+        if missing:
+            print(f"MISSING {path}: {missing}")
+            bad += 1
+            continue
+        u = doc["unreal"]["centre_cm"]
+        geo = doc["geographic"]
+        expect_x = (geo["centre_easting_m"] - oc.ORIGIN_E) * oc.METRES_TO_UNREAL
+        expect_y = -(geo["centre_northing_m"] - oc.ORIGIN_N) * oc.METRES_TO_UNREAL
+        if abs(u[0] - expect_x) > 1.0 or abs(u[1] - expect_y) > 1.0:
+            print(f"STALE {path}: centre_cm {u} disagrees with the origin conversion "
+                  f"({expect_x:.1f}, {expect_y:.1f})")
+            bad += 1
+            continue
+        print(f"ok {doc['asset']:24} {doc['kind']:10} centre_cm {u[0]:10.1f} {u[1]:10.1f}  "
+              f"src {geo['source'][:40]}")
+    print(f"\n{len(files)} sidecars, {bad} problem(s)")
+    return 1 if bad else 0
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    sub = ap.add_subparsers(dest="cmd", required=True)
+    w = sub.add_parser("write")
+    w.add_argument("--asset", required=True)
+    w.add_argument("--kind", required=True, choices=("building", "terrain", "street", "prop", "character"))
+    w.add_argument("--bounds", required=True, help="what the asset was modelled against")
+    w.add_argument("--centre-easting", type=float, required=True)
+    w.add_argument("--centre-northing", type=float, required=True)
+    w.add_argument("--elevation", type=float, default=None, help="NAVD88 metres")
+    w.add_argument("--report", default=None, help="dimension report produced by the build script")
+    w.add_argument("--note", action="append")
+    w.add_argument("--out", required=True)
+    a = sub.add_parser("audit")
+    a.add_argument("targets", nargs="+")
+    args = ap.parse_args()
+    if args.cmd == "write":
+        cmd_write(args)
+        return 0
+    return cmd_audit(args.targets)
+
+
+if __name__ == "__main__":
+    sys.exit(main())

Geo/tools/block_height_report.py

@@ -0,0 +1,75 @@
+#!/usr/bin/env python3
+"""Report which block buildings have real height data and which are guesses.
+
+Usage:
+    Geo/tools/block_height_report.py Geo/raw/osm-api-map-2026-11-06.osm
+
+`height` and `building:levels` are the only height evidence OSM carries. Where
+neither exists, Blosm falls back to `defaultLevels * levelHeight`, which is a
+placeholder: the demo renders a constant-height block instead of a real skyline.
+This tool lists every building in the demo bounds with its height source, so the
+substitutions are visible and measurable rather than buried in the model.
+
+Exit code is 0 either way; the count of guessed buildings is the deliverable.
+"""
+
+import argparse
+import math
+import sys
+
+import osm_common as oc
+
+LEVEL_HEIGHT = 3.2
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("osm")
+    args = ap.parse_args()
+    nodes, ways = oc.load(args.osm)
+
+    line, lat0, lon0, lengths, block_len = oc.block_extent(
+        nodes, ways, "Hollywood Boulevard", "North Orange Drive", "North Highland Avenue")
+
+    def in_block(pts):
+        if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
+            return False
+        return True
+
+    rows = []
+    for wid, tags, refs in ways:
+        if "building" not in tags:
+            continue
+        pts = oc.way_points(nodes, refs)
+        if not pts or not in_block(pts):
+            continue
+        tag_height = tags.get("height")
+        levels = tags.get("building:levels")
+        if tag_height:
+            height, source = float(tag_height), "height tag"
+        elif levels:
+            height, source = float(levels) * LEVEL_HEIGHT, f"levels x {LEVEL_HEIGHT} m"
+        else:
+            height, source = 2 * LEVEL_HEIGHT, f"GUESS: default 2 levels x {LEVEL_HEIGHT} m"
+        nearest = min(oc.project(p, line, lat0, lon0, lengths)[0] for p in pts)
+        side = "S" if oc.project(pts[0], line, lat0, lon0, lengths)[2] < 0 else "N"
+        rows.append((side, nearest, wid, tags.get("name", "") or tags.get("addr:housenumber", ""),
+                     height, source))
+
+    rows.sort(key=lambda r: (r[0], r[1]))
+    guessed = [r for r in rows if r[5].startswith("GUESS")]
+    print(f"{len(rows)} buildings in the block, {len(guessed)} with a guessed height\n")
+    for side, setback, wid, name, height, source in rows:
+        flag = "  <-- guessed" if source.startswith("GUESS") else ""
+        print(f"{side} way {wid:>10}  setback {setback:5.1f} m  height {height:5.1f} m  "
+              f"({source:26})  {name[:34]}{flag}")
+
+    if guessed:
+        print("\nGuessed heights to replace from the reference photography:")
+        for side, setback, wid, name, height, source in guessed:
+            print(f"  {side} way {wid} {name}")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

Geo/tools/block_report.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+"""Report streets, boundary intersections and the demo building set for an OSM extract.
+
+Usage:
+    Geo/tools/block_report.py Geo/raw/osm-api-map-2026-11-06.osm \
+        --blvd "Hollywood Boulevard" --west "North Orange Drive" --east "North Highland Avenue"
+
+Reads a plain .osm XML file (the output of the OSM API `map` endpoint or an
+Overpass export). No network access, no dependencies beyond the standard library.
+
+Output: the boulevard's block span in metres, the two boundary intersections, the
+buildings on the boulevard frontage split by side (even/odd is not trusted —
+latitude relative to the centreline is), and each building's distance from the
+centreline. Buildings behind the frontage are listed but marked as off-buffer.
+"""
+
+import argparse
+import math
+import xml.etree.ElementTree as ET
+
+EARTH_R = 6378137.0
+
+
+def local_metres(lat, lon, lat0, lon0):
+    x = math.radians(lon - lon0) * EARTH_R * math.cos(math.radians(lat0))
+    y = math.radians(lat - lat0) * EARTH_R
+    return x, y
+
+
+def load(path):
+    root = ET.parse(path).getroot()
+    nodes = {n.get("id"): (float(n.get("lat")), float(n.get("lon"))) for n in root.findall("node")}
+    ways = []
+    for w in root.findall("way"):
+        tags = {c.get("k"): c.get("v") for c in w.findall("tag")}
+        ways.append((w.get("id"), tags, [nd.get("ref") for nd in w.findall("nd")]))
+    return nodes, ways
+
+
+def refs_by_name(ways, name):
+    return {ref for _, tags, refs in ways if tags.get("name") == name for ref in refs}
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("osm")
+    ap.add_argument("--blvd", default="Hollywood Boulevard")
+    ap.add_argument("--west", default="North Orange Drive")
+    ap.add_argument("--east", default="North Highland Avenue")
+    ap.add_argument("--buffer", type=float, default=50.0, help="north-side buffer in metres, from the centreline")
+    args = ap.parse_args()
+
+    nodes, ways = load(args.osm)
+    blvd = refs_by_name(ways, args.blvd)
+    if not blvd:
+        raise SystemExit(f"no highway named {args.blvd!r} in {args.osm}")
+
+    ends = {}
+    for label, street in (("west", args.west), ("east", args.east)):
+        inter = refs_by_name(ways, street) & blvd
+        pts = sorted(nodes[i][1] for i in inter if i in nodes)
+        if not pts:
+            raise SystemExit(f"no intersection between {args.blvd!r} and {street!r}")
+        ends[label] = pts[0] if label == "west" else pts[-1]
+        print(f"{street} x {args.blvd}: {len(pts)} shared node(s), lon {pts[0]:.6f}..{pts[-1]:.6f}")
+
+    lo_min, lo_max = ends["west"], ends["east"]
+    lat0 = sum(nodes[i][0] for i in blvd if lo_min <= nodes[i][1] <= lo_max) / max(
+        1, len([i for i in blvd if lo_min <= nodes[i][1] <= lo_max])
+    )
+    span = math.radians(lo_max - lo_min) * EARTH_R * math.cos(math.radians(lat0))
+    print(f"block span: {span:.1f} m of frontage, centreline lat {lat0:.6f}\n")
+
+    rows = []
+    for wid, tags, refs in ways:
+        if "building" not in tags:
+            continue
+        pts = [nodes[r] for r in refs if r in nodes]
+        if not pts:
+            continue
+        lat = sum(p[0] for p in pts) / len(pts)
+        lon = sum(p[1] for p in pts) / len(pts)
+        if not (lo_min - 0.00005 <= lon <= lo_max + 0.00005):
+            continue
+        side = "S" if lat < lat0 else "N"
+        _, y = local_metres(lat, lon, lat0, lon)
+        rows.append((side, abs(y), lon, wid, tags))
+
+    rows.sort(key=lambda r: r[2])
+    for side, dist, lon, wid, tags in rows:
+        if side == "S":
+            flag = ""
+        else:
+            flag = "" if dist <= args.buffer else "  (outside buffer)"
+        print(
+            f"{side} way {wid:>10}  {tags.get('addr:housenumber', '-'):>11} "
+            f"{tags.get('addr:street', '')[:22]:22} h={tags.get('height', '-'):>5} "
+            f"lv={tags.get('building:levels', '-'):>3} d={dist:5.1f}m {tags.get('name', '')[:30]}{flag}"
+        )
+
+
+if __name__ == "__main__":
+    main()

Geo/tools/dem_stats.py

@@ -0,0 +1,173 @@
+#!/usr/bin/env python3
+"""Read a plain float32 GeoTIFF, report its georeferencing and elevation stats.
+
+Usage:
+    Geo/tools/dem_stats.py Geo/raw/dem-3dep-1m-26911-2026-11-06.tif
+    Geo/tools/dem_stats.py <tif> --sample 376515.389,3774224.083
+
+Why this exists: GDAL is not runnable on the host (its flatpak libraries do not
+resolve outside the sandbox), so this reads the strip-based GeoTIFF that the
+3DEP export service produces with the standard library only. It handles
+uncompressed, single-sample, strip-organised float32 images with the usual
+GeoTIFF tags; it is a small reader, not a general GeoTIFF implementation.
+
+Values are orthometric metres (NAVD88 for 3DEP CONUS). Nodata is excluded from
+the statistics.
+"""
+
+import argparse
+import struct
+import sys
+
+TYPE_SIZES = {1: 1, 2: 1, 3: 2, 4: 4, 5: 8, 7: 1, 11: 4, 12: 8, 16: 8}
+
+
+def read_tiff(path):
+    with open(path, "rb") as handle:
+        blob = handle.read()
+
+    if blob[:2] == b"II":
+        endian = "<"
+    elif blob[:2] == b"MM":
+        endian = ">"
+    else:
+        raise SystemExit("not a TIFF file")
+    if struct.unpack(endian + "H", blob[2:4])[0] != 42:
+        raise SystemExit("not a classic TIFF (BigTIFF is unsupported)")
+
+    ifd_offset = struct.unpack(endian + "I", blob[4:8])[0]
+    count = struct.unpack(endian + "H", blob[ifd_offset : ifd_offset + 2])[0]
+    entries = {}
+    for i in range(count):
+        base = ifd_offset + 2 + i * 12
+        tag, typ, n = struct.unpack(endian + "HHI", blob[base : base + 8])
+        size = TYPE_SIZES.get(typ, 1) * n
+        raw = blob[base + 8 : base + 8 + 4] if size <= 4 else blob[
+            struct.unpack(endian + "I", blob[base + 8 : base + 12])[0] :
+        ][:size]
+        fmt = {1: "B", 3: "H", 4: "I", 7: "B", 11: "f", 12: "d", 2: "c"}.get(typ, "B")
+        if typ == 2:
+            values = raw[:size].split(b"\0")[0].decode("latin-1")
+        elif fmt == "c":
+            values = raw
+        else:
+            values = list(struct.unpack(endian + fmt * n, raw[: struct.calcsize(endian + fmt * n)]))
+        entries[tag] = values
+
+    width = entries[256][0]
+    height = entries[257][0]
+    bits = entries.get(258, [8])[0]
+    compression = entries.get(259, [1])[0]
+    sample_format = entries.get(339, [1])[0]
+    if compression != 1 or bits != 32 or sample_format != 3:
+        raise SystemExit(
+            f"unsupported TIFF layout: compression={compression} bits={bits} sample_format={sample_format}"
+        )
+
+    if 322 in entries:
+        # tiled layout: tags 322/323 give tile size, 324/325 the offsets and byte counts
+        tile_w = entries[322][0]
+        tile_h = entries[323][0]
+        offsets = entries[324]
+        counts = entries[325]
+        across = (width + tile_w - 1) // tile_w
+        down = (height + tile_h - 1) // tile_h
+        if len(offsets) < across * down:
+            raise SystemExit(f"tiled TIFF with {len(offsets)} tiles, expected {across * down}")
+        grid = [[0.0] * width for _ in range(height)]
+        for index, (offset, nbytes) in enumerate(zip(offsets, counts)):
+            tile_row, tile_col = divmod(index, across)
+            tile = struct.unpack(endian + "f" * (tile_w * tile_h), blob[offset : offset + nbytes])
+            for r in range(tile_h):
+                target_row = tile_row * tile_h + r
+                if target_row >= height:
+                    break
+                base = r * tile_w
+                for c in range(tile_w):
+                    target_col = tile_col * tile_w + c
+                    if target_col >= width:
+                        break
+                    grid[target_row][target_col] = tile[base + c]
+        flat = [value for row in grid for value in row]
+        return entries, width, height, flat
+
+    # strip layout
+    offsets = entries[273]
+    counts = entries[279]
+    rows_per_strip = entries.get(278, [height])[0]
+    chunks = [blob[offset : offset + nbytes] for offset, nbytes in zip(offsets, counts)]
+    pixels = b"".join(chunks)
+    values = list(struct.unpack(endian + "f" * (width * height), pixels[: width * height * 4]))
+    return entries, width, height, values
+
+
+def geokeys(entries):
+    keys = {}
+    raw = entries.get(34735)
+    if not raw:
+        return keys
+    # GeoKeyDirectory: version, revision, minor, count, then 4 shorts per key
+    for i in range(4, len(raw), 4):
+        key, tag_location, count, value = raw[i : i + 4]
+        keys[key] = (tag_location, count, value)
+    return keys
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("tif")
+    ap.add_argument("--sample", action="append", default=[], help="E,N to sample, repeatable")
+    args = ap.parse_args()
+
+    entries, width, height, values = read_tiff(args.tif)
+    scale = entries.get(33550)
+    tiepoint = entries.get(33922)
+    keys = geokeys(entries)
+    nodata = None
+    if 42113 in entries:
+        try:
+            nodata = float(entries[42113])
+        except ValueError:
+            nodata = None
+
+    print(f"size: {width} x {height} px, float32, uncompressed")
+    if nodata is not None:
+        print(f"nodata: {nodata}")
+
+    valid = [v for v in values if nodata is None or v != nodata]
+    if valid:
+        print(
+            f"elevation: min {min(valid):.2f} m, max {max(valid):.2f} m, "
+            f"relief {max(valid) - min(valid):.2f} m, mean {sum(valid)/len(valid):.2f} m"
+        )
+        print(f"valid pixels: {len(valid)} of {width*height}")
+
+    if scale and tiepoint:
+        sx, sy = scale[0], scale[1]
+        px, py, pz, ex, ey = tiepoint[0], tiepoint[1], tiepoint[2], tiepoint[3], tiepoint[4]
+        x0 = ex - px * sx
+        y0 = ey + py * sy  # north-up: tiepoint row 0 is the top edge
+        print(f"pixel size: {sx} x {sy} m")
+        print(f"extent: E {x0:.1f} .. {x0 + width*sx:.1f}, N {y0 - height*sy:.1f} .. {y0:.1f}")
+        crs = keys.get(3072, (None, None, None))[2]
+        model = keys.get(1024, (None, None, None))[2]
+        print(f"model type: {model} (1 = projected), projected CRS: EPSG:{crs}")
+
+        for spec in args.sample:
+            try:
+                e, n = (float(x) for x in spec.split(","))
+            except ValueError:
+                raise SystemExit(f"bad --sample value: {spec!r}")
+            col = int((e - x0) / sx)
+            row = int((y0 - n) / sy)
+            if not (0 <= col < width and 0 <= row < height):
+                print(f"  sample E{e:.1f} N{n:.1f}: outside the raster")
+                continue
+            value = values[row * width + col]
+            print(f"  sample E{e:.3f} N{n:.3f} -> row {row}, col {col}: {value:.3f} m")
+    else:
+        print("no pixel scale or tiepoint found; georeferencing missing", file=sys.stderr)
+
+
+if __name__ == "__main__":
+    main()

Geo/tools/elevation_probe.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+"""Probe ground elevation across the demo block with the USGS 3DEP point service.
+
+Usage:
+    Geo/tools/elevation_probe.py Geo/raw/osm-api-map-2026-11-06.osm
+
+Samples a grid along the boulevard (centreline and both sidewalk centrelines),
+the hero footprint corners, and the crossing midpoints, then reports the relief
+across the demo and whether a separate elevation source is needed.
+
+The decision rule follows `AGENTS.md`: acquire separately licensed elevation data
+unless the DEM shows less than 1 m of height variation across the demo bounds.
+
+Writes the raw responses to `Geo/raw/epqs-<date>.json` so the numbers can be
+re-checked without re-querying. Service: USGS Elevation Point Query Service,
+3DEP 1 m seamless DEM (CONUS). Standard library only.
+"""
+
+import argparse
+import concurrent.futures
+import datetime
+import json
+import os
+import time
+import urllib.parse
+import urllib.request
+
+import osm_common as oc
+
+EPQS = "https://epqs.nationalmap.gov/v1/json"
+UA = "gig-line-geo/0.1 (local development, elevation check)"
+
+
+def query(lat: float, lon: float, retries: int = 5) -> dict:
+    params = urllib.parse.urlencode(
+        {"x": f"{lon:.7f}", "y": f"{lat:.7f}", "wkid": "4326", "units": "Meters", "includeDate": "false"}
+    )
+    url = f"{EPQS}?{params}"
+    for attempt in range(retries):
+        req = urllib.request.Request(url, headers={"User-Agent": UA})
+        try:
+            with urllib.request.urlopen(req, timeout=40) as resp:
+                payload = json.loads(resp.read().decode("utf-8"))
+            value = payload.get("value")
+            if value in (None, "", "No Data"):
+                raise ValueError(f"no elevation returned: {payload}")
+            return {
+                "lat": round(lat, 7),
+                "lon": round(lon, 7),
+                "elevation_m": float(value),
+                "raster_id": payload.get("rasterId"),
+                "resolution_m": payload.get("resolution"),
+            }
+        except Exception as exc:  # noqa: BLE001
+            if attempt == retries - 1:
+                # a single dead point must not kill the whole probe
+                print(f"    ! no elevation at {lat:.6f},{lon:.6f}: {exc}", flush=True)
+                return {"lat": round(lat, 7), "lon": round(lon, 7), "elevation_m": None,
+                        "raster_id": None, "resolution_m": None}
+            time.sleep(3 * (attempt + 1))
+    return {}
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("osm")
+    ap.add_argument("--blvd", default="Hollywood Boulevard")
+    ap.add_argument("--west", default="North Orange Drive")
+    ap.add_argument("--east", default="North Highland Avenue")
+    ap.add_argument("--hero", type=int, default=424024147, help="hero building OSM way id")
+    ap.add_argument("--workers", type=int, default=4, help="parallel requests; the service is slow")
+    ap.add_argument("--pause", type=float, default=0.0)
+    args = ap.parse_args()
+
+    nodes, ways = oc.load(args.osm)
+    line, lat0, lon0, lengths, block_len = oc.block_extent(nodes, ways, args.blvd, args.west, args.east)
+
+    samples = []
+    grid = [block_len * i / 4 for i in range(5)]
+    for s in grid:
+        for label, offset in (("south sidewalk", -12.6), ("centreline", 0.0), ("north sidewalk", 13.6)):
+            lat, lon = oc.point_at(line, lengths, s, offset, lat0, lon0)
+            samples.append({"label": f"s={s:5.1f} m {label}", "lat": lat, "lon": lon})
+
+    hero_refs = next((refs for wid, _, refs in ways if str(wid) == str(args.hero)), None)
+    if hero_refs:
+        hero_pts = oc.way_points(nodes, hero_refs)
+        lats = [p[0] for p in hero_pts]
+        lons = [p[1] for p in hero_pts]
+        corners = [
+            ("near NW", min(lats), min(lons)),
+            ("near NE", min(lats), max(lons)),
+            ("far SW", max(lats), min(lons)),
+            ("far SE", max(lats), max(lons)),
+        ]
+        for label, lat, lon in corners:
+            samples.append({"label": f"hero {label}", "lat": lat, "lon": lon})
+
+    for wid, tags, refs in ways:
+        if tags.get("footway") != "crossing":
+            continue
+        pts = oc.way_points(nodes, refs)
+        if not pts:
+            continue
+        mid = (sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts))
+        dist, s, _ = oc.project(mid, line, lat0, lon0, lengths)
+        # only crossings of THIS boulevard, not side-street crossings inside the window
+        if 0 <= s <= block_len and dist <= 4.0:
+            samples.append({"label": f"crossing way {wid} at s={s:.0f} m", "lat": mid[0], "lon": mid[1]})
+
+    def fetch(sample):
+        r = query(sample["lat"], sample["lon"])
+        r["label"] = sample["label"]
+        return r
+
+    started = time.time()
+    with concurrent.futures.ThreadPoolExecutor(max_workers=max(1, args.workers)) as pool:
+        futures = {pool.submit(fetch, s): s for s in samples}
+        results = []
+        for done, future in enumerate(concurrent.futures.as_completed(futures), start=1):
+            r = future.result()
+            results.append(r)
+            value = "no data" if r["elevation_m"] is None else f"{r['elevation_m']:7.2f} m"
+            print(f"  [{done:02d}/{len(samples)}] {r['label']:28} {value:>10}", flush=True)
+            if args.pause:
+                time.sleep(args.pause)
+    order = {s["label"]: i for i, s in enumerate(samples)}
+    results.sort(key=lambda r: order.get(r["label"], 999))
+    print(f"query time {time.time() - started:.0f} s for {len(samples)} samples")
+
+    walkable = [r for r in results if r["elevation_m"] is not None and not r["label"].startswith("hero")]
+    building = [r for r in results if r["elevation_m"] is not None and r["label"].startswith("hero")]
+    values = [r["elevation_m"] for r in walkable]
+    lo, hi = min(values), max(values)
+    relief = hi - lo
+
+    print("\nsamples by group (walkable surfaces only below; building interiors are diagnostics):")
+    for key in ("crossing", "hero"):
+        group = [r for r in results if r["elevation_m"] is not None and key in r["label"]]
+        if group:
+            g = [r["elevation_m"] for r in group]
+            print(f"  {key:12} n={len(group):2}  min {min(g):.2f}  max {max(g):.2f}  spread {max(g)-min(g):.2f} m")
+
+    print("\ncross-slope check (south sidewalk minus north sidewalk):")
+    for s in [block_len * i / 4 for i in range(5)]:
+        south = next(r for r in results if r["label"] == f"s={s:5.1f} m south sidewalk")
+        north = next(r for r in results if r["label"] == f"s={s:5.1f} m north sidewalk")
+        drop = south["elevation_m"] - north["elevation_m"]
+        print(f"  s={s:5.1f} m: south {south['elevation_m']:.2f} m, north {north['elevation_m']:.2f} m, "
+              f"difference {drop:+.2f} m over ~26 m")
+
+    if building:
+        print("\nwarning: points inside the hero footprint are not ground reference; 3DEP can return")
+        print("         building or void values there. They are excluded from the relief figure.")
+    print(f"\nelevation range on walkable surfaces: {lo:.2f} to {hi:.2f} m, relief {relief:.2f} m")
+    first = next(r for r in walkable)
+    print(f"source: USGS 3DEP point service, raster {first['raster_id']}, resolution {first['resolution_m']} m")
+    if relief < 1.0:
+        print("decision: relief under 1 m, a separate elevation source is NOT required (AGENTS.md rule)")
+    else:
+        print("decision: relief at or above 1 m, a separately licensed elevation source IS required")
+    print("datum: 3DEP CONUS heights are orthometric NAVD88; treat absolute values as such, not ellipsoidal")
+
+    today = datetime.date.today().isoformat()
+    out_path = os.path.join("Geo", "raw", f"epqs-{today}.json")
+    os.makedirs(os.path.dirname(out_path), exist_ok=True)
+    payload = {
+        "service": EPQS,
+        "processed_utc": datetime.datetime.now(datetime.timezone.utc).isoformat(timespec="seconds"),
+        "source_dem": "USGS 3DEP seamless 1 m (CONUS)",
+        "vertical_datum": "NAVD88 orthometric (3DEP CONUS)",
+        "query": {
+            "blvd": args.blvd,
+            "west": args.west,
+            "east": args.east,
+            "hero_way": args.hero,
+            "block_len_m": round(block_len, 2),
+        },
+        "summary": {"min_m": lo, "max_m": hi, "relief_m": round(relief, 3), "samples": len(results)},
+        "samples": results,
+    }
+    with open(out_path, "w", encoding="utf-8") as handle:
+        json.dump(payload, handle, indent=2)
+    print(f"wrote {out_path}")
+
+
+if __name__ == "__main__":
+    main()

Geo/tools/osm_common.py

@@ -0,0 +1,226 @@
+#!/usr/bin/env python3
+"""Shared geo helpers: parse .osm, stitch a street, project points, convert to Unreal.
+
+Standard library only. Imported by every script in `Geo/tools/`.
+
+The unit boundary is deliberate and lives here, once:
+
+    WGS84 (lat, lon)  ->  EPSG:26911 (E, N metres)  ->  project local metres
+                      ->  Unreal centimetres
+
+`to_unreal` is the only place that multiplies by 100. Nothing else in the
+repository may convert metres to centimetres, so a second silent scale cannot be
+introduced by accident.
+"""
+
+from __future__ import annotations
+
+import math
+import os
+import xml.etree.ElementTree as ET
+
+EARTH_R = 6378137.0
+
+# Project origin, EPSG:26911. See Geo/manifests/local-origin.md.
+ORIGIN_E = 376515.389
+ORIGIN_N = 3774224.083
+ORIGIN_Z_NAVD88 = 115.911
+METRES_TO_UNREAL = 100.0
+
+
+def to_unreal_local(easting, northing, elevation=None):
+    """EPSG:26911 metres -> Unreal centimetres in the project frame.
+
+    X = east of the origin, Y = south of the origin (the FBX preset maps Blender
+    +Y/north to Unreal -Y), Z = height above the origin ground. This is the single
+    meters->centimetres conversion for the repository.
+    """
+    x_cm = (easting - ORIGIN_E) * METRES_TO_UNREAL
+    y_cm = -(northing - ORIGIN_N) * METRES_TO_UNREAL
+    if elevation is None:
+        return x_cm, y_cm, None
+    return x_cm, y_cm, (elevation - ORIGIN_Z_NAVD88) * METRES_TO_UNREAL
+
+
+def to_blender_local(easting, northing, elevation=None):
+    """EPSG:26911 metres -> Blender metres in the project frame (X east, Y north, Z up)."""
+    x = easting - ORIGIN_E
+    y = northing - ORIGIN_N
+    z = None if elevation is None else elevation - ORIGIN_Z_NAVD88
+    return x, y, z
+
+
+def find_proj():
+    """The PROJ binary: environment override, PATH, then the QGIS Flatpak."""
+    import glob
+    import os
+    import shutil
+
+    env = os.environ.get("PROJ_BIN")
+    if env and os.path.exists(env):
+        return env
+    which = shutil.which("proj")
+    if which:
+        return which
+    matches = sorted(glob.glob(os.path.expanduser(
+        "~/.local/share/flatpak/app/org.qgis.qgis/x86_64/stable/*/files/bin/proj")))
+    if not matches:
+        raise SystemExit("no PROJ binary found; set PROJ_BIN or install pyproj")
+    return matches[0]
+
+
+def project_wgs84_to_utm(points, proj_bin=None):
+    """WGS84 (lat, lon) -> EPSG:26911 (E, N) for a list, in one PROJ call."""
+    import subprocess
+
+    proj_bin = proj_bin or find_proj()
+    app_files = os.path.dirname(os.path.dirname(proj_bin))
+    env = dict(os.environ)
+    env["LD_LIBRARY_PATH"] = f"{app_files}/lib:{env.get('LD_LIBRARY_PATH', '')}"
+    payload = "\n".join(f"{lon:.9f} {lat:.9f}" for lat, lon in points) + "\n"
+    out = subprocess.run(
+        [proj_bin, "+proj=utm", "+zone=11", "+datum=NAD83", "-f", "%.3f"],
+        input=payload, capture_output=True, text=True, env=env, check=True,
+    ).stdout.strip().splitlines()
+    result = [tuple(float(v) for v in line.split()) for line in out]
+    if len(result) != len(points):
+        raise SystemExit(f"PROJ returned {len(result)} points for {len(points)} inputs")
+    return result
+
+
+def to_xy(lat: float, lon: float, lat0: float, lon0: float) -> tuple[float, float]:
+    """Local planar metres, east/north, centred on (lat0, lon0)."""
+    return (
+        math.radians(lon - lon0) * EARTH_R * math.cos(math.radians(lat0)),
+        math.radians(lat - lat0) * EARTH_R,
+    )
+
+
+def from_xy(x: float, y: float, lat0: float, lon0: float) -> tuple[float, float]:
+    lat = lat0 + math.degrees(y / EARTH_R)
+    lon = lon0 + math.degrees(x / (EARTH_R * math.cos(math.radians(lat0))))
+    return lat, lon
+
+
+def load(path: str):
+    root = ET.parse(path).getroot()
+    nodes = {n.get("id"): (float(n.get("lat")), float(n.get("lon"))) for n in root.findall("node")}
+    ways = []
+    for w in root.findall("way"):
+        tags = {c.get("k"): c.get("v") for c in w.findall("tag")}
+        ways.append((w.get("id"), tags, [nd.get("ref") for nd in w.findall("nd")]))
+    return nodes, ways
+
+
+def way_points(nodes, refs):
+    return [nodes[r] for r in refs if r in nodes]
+
+
+def refs_by_name(ways, name):
+    return {r for _, tags, refs in ways if tags.get("name") == name for r in refs}
+
+
+def stitch(nodes, ways, name, lon_min, lon_max):
+    """Ordered polyline for a named street, restricted to a longitude window."""
+    segments = []
+    for _, tags, refs in ways:
+        if tags.get("name") != name:
+            continue
+        pts = way_points(nodes, refs)
+        if not pts:
+            continue
+        if max(p[1] for p in pts) < lon_min - 1e-6 or min(p[1] for p in pts) > lon_max + 1e-6:
+            continue
+        segments.append(pts)
+
+    merged = []
+    for pts in segments:
+        if not merged:
+            merged = list(pts)
+            continue
+        if pts[0] == merged[-1]:
+            merged.extend(pts[1:])
+        elif pts[-1] == merged[-1]:
+            merged.extend(reversed(pts[:-1]))
+        elif pts[0] == merged[0]:
+            merged = list(reversed(pts)) + merged[1:]
+        elif pts[-1] == merged[0]:
+            merged = pts + merged[1:]
+    return merged
+
+
+def cumulative_lengths(points, lat0, lon0):
+    xs, total, prev = [], 0.0, None
+    for lat, lon in points:
+        p = to_xy(lat, lon, lat0, lon0)
+        if prev is not None:
+            total += math.dist(prev, p)
+        xs.append(total)
+        prev = p
+    return xs
+
+
+def project(point, line, lat0, lon0, lengths):
+    """Nearest point on the polyline: (distance_m, s_along_m, side) with side +1 = north."""
+    px, py = to_xy(point[0], point[1], lat0, lon0)
+    best = (float("inf"), 0.0, 0)
+    for i in range(len(line) - 1):
+        ax, ay = to_xy(line[i][0], line[i][1], lat0, lon0)
+        bx, by = to_xy(line[i + 1][0], line[i + 1][1], lat0, lon0)
+        dx, dy = bx - ax, by - ay
+        seg2 = dx * dx + dy * dy
+        if seg2 == 0:
+            continue
+        t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / seg2))
+        cx, cy = ax + t * dx, ay + t * dy
+        dist = math.hypot(px - cx, py - cy)
+        if dist < best[0]:
+            cross = dx * (py - ay) - dy * (px - ax)
+            best = (dist, lengths[i] + t * (lengths[i + 1] - lengths[i]), 1 if cross > 0 else -1)
+    return best
+
+
+def point_at(line, lengths, s, offset, lat0, lon0):
+    """Point on the polyline at arc length s, pushed offset metres to the north side."""
+    for i in range(len(line) - 1):
+        if lengths[i] <= s <= lengths[i + 1]:
+            seg = lengths[i + 1] - lengths[i]
+            t = 0.0 if seg == 0 else (s - lengths[i]) / seg
+            ax, ay = to_xy(line[i][0], line[i][1], lat0, lon0)
+            bx, by = to_xy(line[i + 1][0], line[i + 1][1], lat0, lon0)
+            dx, dy = bx - ax, by - ay
+            norm = math.hypot(dx, dy) or 1.0
+            nx, ny = -dy / norm, dx / norm
+            x = ax + t * dx + offset * nx
+            y = ay + t * dy + offset * ny
+            return from_xy(x, y, lat0, lon0)
+    return line[-1]
+
+
+def block_extent(nodes, ways, blvd, west, east):
+    """Stitch and clip the boulevard between two cross streets.
+
+    Returns (line, lat0, lon0, lengths, block_len) with `line` running west to east.
+    """
+    inter = {}
+    blvd_refs = refs_by_name(ways, blvd)
+    for label, street in (("west", west), ("east", east)):
+        shared = [nodes[r] for r in refs_by_name(ways, street) & blvd_refs if r in nodes]
+        if not shared:
+            raise SystemExit(f"no intersection between {blvd!r} and {street!r}")
+        inter[label] = min(p[1] for p in shared) if label == "west" else max(p[1] for p in shared)
+
+    line = stitch(nodes, ways, blvd, inter["west"], inter["east"])
+    if line[0][1] > line[-1][1]:
+        line = list(reversed(line))
+    west_pt = min(line, key=lambda p: abs(p[1] - inter["west"]))
+    east_pt = min(line, key=lambda p: abs(p[1] - inter["east"]))
+    i0, i1 = line.index(west_pt), line.index(east_pt)
+    if i0 > i1:
+        i0, i1 = i1, i0
+    line = line[i0 : i1 + 1]
+
+    lat0 = sum(p[0] for p in line) / len(line)
+    lon0 = sum(p[1] for p in line) / len(line)
+    lengths = cumulative_lengths(line, lat0, lon0)
+    return line, lat0, lon0, lengths, lengths[-1]

Geo/tools/street_profile.py

@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+"""Street profile for a block: frontage setbacks, carriageway width, crossings.
+
+Usage:
+    Geo/tools/street_profile.py Geo/raw/osm-api-map-2026-11-06.osm \
+        --blvd "Hollywood Boulevard" --west "North Orange Drive" --east "North Highland Avenue"
+
+Measures, in metres, from a plain .osm extract (standard library only):
+
+* the boulevard centreline, stitched and clipped between the two boundary
+  intersections;
+* each building's frontage setback: the distance from its nearest vertex to the
+  centreline, and which side of the road it is on;
+* each mapped crossing (`footway=crossing`) across the boulevard: its length,
+  which is the carriageway width at that point, and its position along the block;
+* the offset of the mapped sidewalk centrelines from the road centreline.
+
+What it cannot measure: carriageway width where no crossing or kerb is mapped,
+and sidewalk width, because OSM models sidewalks here as centre lines with no
+width tag. Derive those from the numbers printed here and record them as
+assumptions.
+"""
+
+import argparse
+
+import osm_common as oc
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("osm")
+    ap.add_argument("--blvd", default="Hollywood Boulevard")
+    ap.add_argument("--west", default="North Orange Drive")
+    ap.add_argument("--east", default="North Highland Avenue")
+    args = ap.parse_args()
+
+    nodes, ways = oc.load(args.osm)
+    line, lat0, lon0, lengths, block_len = oc.block_extent(
+        nodes, ways, args.blvd, args.west, args.east
+    )
+    print(f"centreline: {len(line)} vertices, {block_len:.1f} m between the boundary intersections\n")
+
+    print("crossings across the boulevard (length = carriageway width there):")
+    for wid, tags, refs in ways:
+        if tags.get("footway") != "crossing":
+            continue
+        pts = oc.way_points(nodes, refs)
+        if not pts:
+            continue
+        if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
+            continue
+        a = oc.to_xy(pts[0][0], pts[0][1], lat0, lon0)
+        b = oc.to_xy(pts[-1][0], pts[-1][1], lat0, lon0)
+        length = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5
+        mid = (sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts))
+        dist, s, _ = oc.project(mid, line, lat0, lon0, lengths)
+        if not (0.0 <= s <= block_len) or dist > 4.0:
+            continue
+        print(
+            f"  way {wid:>10}  {length:6.1f} m across, at {s:6.1f} m along the block "
+            f"(offset {dist:4.1f} m)  markings={tags.get('crossing:markings', '-'):6} "
+            f"signals={tags.get('crossing')}"
+        )
+
+    print("\nsidewalk centrelines near the block (offset from road centreline):")
+    south, north = [], []
+    for wid, tags, refs in ways:
+        if tags.get("footway") != "sidewalk":
+            continue
+        pts = oc.way_points(nodes, refs)
+        if not pts:
+            continue
+        if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
+            continue
+        offsets = [oc.project(p, line, lat0, lon0, lengths) for p in pts]
+        mean = sum(o[0] for o in offsets) / len(offsets)
+        along = max(o[1] for o in offsets) - min(o[1] for o in offsets)
+        if not (6.0 <= mean <= 20.0) or along < 30.0:
+            continue
+        (south if offsets[0][2] < 0 else north).append((mean, along, wid, tags.get("name", "")))
+
+    for label, rows in (("south side", south), ("north side", north)):
+        if not rows:
+            print(f"  {label}: none mapped")
+            continue
+        rows.sort()
+        mean = sum(r[0] for r in rows) / len(rows)
+        print(f"  {label}: {len(rows)} line(s), mean offset {mean:.1f} m, range {rows[0][0]:.1f}..{rows[-1][0]:.1f} m")
+
+    print("\nbuilding frontage setbacks (nearest vertex to centreline):")
+    rows = []
+    for wid, tags, refs in ways:
+        if "building" not in tags:
+            continue
+        pts = oc.way_points(nodes, refs)
+        if not pts:
+            continue
+        if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
+            continue
+        proj = [oc.project(p, line, lat0, lon0, lengths) for p in pts]
+        nearest = min(proj, key=lambda o: o[0])
+        along = [o[1] for o in proj]
+        rows.append((nearest[0], nearest[2], wid, tags, max(along) - min(along)))
+    rows.sort()
+    for dist, side, wid, tags, frontage in rows:
+        side_name = "S" if side < 0 else "N"
+        print(
+            f"  {side_name} way {wid:>10}  {tags.get('addr:housenumber', '-'):>11} "
+            f"setback {dist:5.1f} m, extent {frontage:5.1f} m  {tags.get('name', '')[:34]}"
+        )
+
+    print("\nassumptions to record, not measurements:")
+    print("  - lane width is not tagged; the crossing lengths above are the best carriageway")
+    print("    proxy, and crossing endpoints sit within about 1 m of the kerb line")
+    print("  - sidewalk width is not tagged; it is setback minus half carriageway")
+
+
+if __name__ == "__main__":
+    main()

Geo/tools/street_section.py

@@ -0,0 +1,187 @@
+#!/usr/bin/env python3
+"""Cross-sections of the street from the local 3DEP DEM.
+
+Usage:
+    Geo/tools/street_section.py Geo/raw/osm-api-map-2026-11-06.osm \
+        Geo/raw/dem-3dep-1m-26911-2026-11-06.tif --step 0.5
+
+For several positions along the block it samples the DEM on a line perpendicular
+to the boulevard centreline, from 18 m south to 18 m north, and reports:
+
+* the height profile and total cross-fall;
+* where the surface breaks most sharply, which is the best available candidate
+  for the kerb line;
+* whether a kerb step is visible in the DEM at all.
+
+The OSM geometry is WGS84, the DEM is EPSG:26911, so points are batch-projected
+through the PROJ binary that ships with the QGIS flatpak (the host has no PROJ).
+
+Caveat: 3DEP is bare-earth lidar at 1 m. Kerbs are short vertical steps and are
+usually smoothed or removed by the DEM production process. If no break shows at
+the expected offset, the DEM simply does not resolve kerbs and kerb height has to
+come from a stated assumption instead.
+"""
+
+import argparse
+import glob
+import os
+import shutil
+import subprocess
+import sys
+
+import dem_stats
+import osm_common as oc
+
+
+def find_proj():
+    env = os.environ.get("PROJ_BIN")
+    if env and os.path.exists(env):
+        return env
+    which = shutil.which("proj")
+    if which:
+        return which
+    pattern = os.path.expanduser(
+        "~/.local/share/flatpak/app/org.qgis.qgis/x86_64/stable/*/files/bin/proj"
+    )
+    matches = sorted(glob.glob(pattern))
+    if matches:
+        return matches[0]
+    raise SystemExit(
+        "no PROJ binary found. Set PROJ_BIN, or install pyproj, or use the QGIS flatpak."
+    )
+
+
+def project_batch(points, proj_bin):
+    """WGS84 (lat, lon) -> EPSG:26911 (E, N) for a list of points, one PROJ call."""
+    app_files = os.path.dirname(os.path.dirname(proj_bin))
+    env = dict(os.environ)
+    env["LD_LIBRARY_PATH"] = f"{app_files}/lib:{env.get('LD_LIBRARY_PATH', '')}"
+    payload = "\n".join(f"{lon:.8f} {lat:.8f}" for lat, lon in points) + "\n"
+    result = subprocess.run(
+        [proj_bin, "+proj=utm", "+zone=11", "+datum=NAD83", "-f", "%.4f"],
+        input=payload,
+        capture_output=True,
+        text=True,
+        env=env,
+        check=True,
+    )
+    out = []
+    for line in result.stdout.strip().splitlines():
+        parts = line.split()
+        out.append((float(parts[0]), float(parts[1])))
+    if len(out) != len(points):
+        raise SystemExit(f"PROJ returned {len(out)} points for {len(points)} inputs")
+    return out
+
+
+class Dem:
+    def __init__(self, path):
+        entries, width, height, values = dem_stats.read_tiff(path)
+        self.w, self.h, self.v = width, height, values
+        scale = entries[33550]
+        tie = entries[33922]
+        self.x0 = tie[3] - tie[0] * scale[0]
+        self.y0 = tie[4] + tie[1] * scale[1]
+        self.sx, self.sy = scale[0], scale[1]
+        nodata = entries.get(42113)
+        self.nodata = float(nodata) if nodata not in (None, "") else None
+
+    def at(self, e, n):
+        """Bilinear sample; returns None outside the raster or on nodata."""
+        col = (e - self.x0) / self.sx - 0.5
+        row = (self.y0 - n) / self.sy - 0.5
+        c0, r0 = int(col // 1), int(row // 1)
+        if c0 < 0 or r0 < 0 or c0 + 1 >= self.w or r0 + 1 >= self.h:
+            return None
+        fx, fy = col - c0, row - r0
+        p = []
+        for r, c in ((r0, c0), (r0, c0 + 1), (r0 + 1, c0), (r0 + 1, c0 + 1)):
+            value = self.v[r * self.w + c]
+            if self.nodata is not None and value == self.nodata:
+                return None
+            p.append(value)
+        return (
+            p[0] * (1 - fx) * (1 - fy)
+            + p[1] * fx * (1 - fy)
+            + p[2] * (1 - fx) * fy
+            + p[3] * fx * fy
+        )
+
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("osm")
+    ap.add_argument("dem")
+    ap.add_argument("--blvd", default="Hollywood Boulevard")
+    ap.add_argument("--west", default="North Orange Drive")
+    ap.add_argument("--east", default="North Highland Avenue")
+    ap.add_argument("--step", type=float, default=0.5, help="sampling step in metres")
+    ap.add_argument("--reach", type=float, default=18.0, help="sample +/- this many metres")
+    ap.add_argument("--section", type=float, action="append", default=[],
+                    help="position along the block in metres; repeatable. Default: every 40 m")
+    args = ap.parse_args()
+
+    nodes, ways = oc.load(args.osm)
+    line, lat0, lon0, lengths, block_len = oc.block_extent(
+        nodes, ways, args.blvd, args.west, args.east
+    )
+    dem = Dem(args.dem)
+    proj_bin = find_proj()
+
+    sections = args.section or [40.0 * i for i in range(8) if 40.0 * i <= block_len]
+    offsets = []
+    x = -args.reach
+    while x <= args.reach + 1e-9:
+        offsets.append(round(x, 2))
+        x += args.step
+
+    # build every sample point, project them in one batch
+    grid = []
+    for s in sections:
+        row = []
+        for off in offsets:
+            row.append(oc.point_at(line, lengths, s, off, lat0, lon0))
+        grid.append(row)
+    flat = [p for row in grid for p in row]
+    projected = project_batch(flat, proj_bin)
+
+    print(f"block {block_len:.1f} m, {len(sections)} sections, "
+          f"offsets {offsets[0]:+.1f}..{offsets[-1]:+.1f} m at {args.step} m\n")
+
+    index = 0
+    for s, row in zip(sections, grid):
+        heights = []
+        for _ in row:
+            e, n = projected[index]
+            index += 1
+            heights.append(dem.at(e, n))
+
+        valid = [(o, h) for o, h in zip(offsets, heights) if h is not None]
+        if not valid:
+            print(f"s={s:6.1f} m: no DEM data")
+            continue
+
+        hs = [h for _, h in valid]
+        high_off, high_h = max(valid, key=lambda t: t[1])
+        low_off, low_h = min(valid, key=lambda t: t[1])
+        cross_fall = high_h - low_h
+
+        # sharpest break in the profile = candidate kerb
+        breaks = []
+        for (o1, h1), (o2, h2) in zip(valid, valid[1:]):
+            slope = abs(h2 - h1) / max(1e-9, abs(o2 - o1))
+            breaks.append((slope, (o1 + o2) / 2, h2 - h1))
+        breaks.sort(reverse=True)
+        top = breaks[:3]
+
+        print(f"s={s:6.1f} m  crown {high_h:.2f} m at {high_off:+.1f}  "
+              f"low {low_h:.2f} m at {low_off:+.1f}  cross-fall {cross_fall:.2f} m")
+        print("   profile (2 m): " + " ".join(
+            f"{o:+.0f}:{h:.2f}" for o, h in valid if abs(o % 2) < 1e-6))
+        print("   sharpest breaks: " + ", ".join(
+            f"{slope:.2f} m/m at {off:+.1f} m ({dh:+.2f} m)" for slope, off, dh in top))
+        print()
+
+
+if __name__ == "__main__":
+    sys.exit(main())