feat: terrain mesh from the DEM, and correct the blockout frame

Tools/terrain_from_dem.py builds a 191x106 grid at 2 m over the block, Z being
NAVD88 minus the origin ground. Building footprints from the OSM extract are cut
out so rooftops never become terrain (72 footprints, 6310 cells). Spot checks
reproduce the measured profile: origin 0.00, Orange Dr +3.39, south gutter
+0.93, crown +1.38, north sidewalk +1.43.

Blosm projects with a sphere, not an ellipsoid, so its frame was ~0.4% oversized
and rotated by meridian convergence: about 2 m at the theatre. The blockout now
fits a six-parameter affine into EPSG:26911 on five control points (0.0 cm
residual) and transforms root objects only, since children inherit and double
transforming them read as a 2 m placement error.

Tools/blosm_blockout.py

@@ -41,6 +41,8 @@ def parse_args():
     ap.add_argument("--bounds", required=True, help="south,west,north,east in WGS84 degrees")
     ap.add_argument("--out", required=True, help="output .blend path")
     ap.add_argument("--level-height", type=float, default=3.2, help="metres per building level")
+    ap.add_argument("--origin", default=None,
+                    help="lat,lon to use as the Blender origin (project origin). Default: Blosm centres on the extract")
     ap.add_argument("--no-subdivide", action="store_true", help="skip OSM area subdivision")
     ap.add_argument("--load-missing", action="store_true",
                     help="let Blosm fetch relation members missing from the extract (needs network, breaks offline reproducibility)")
@@ -103,6 +105,13 @@ def import_osm(args):
     scene = bpy.context.scene
     osmp = scene.blosm
 
+    if args.origin:
+        lat, lon = (float(v) for v in args.origin.split(","))
+        scene["lat"] = lat
+        scene["lon"] = lon
+        osmp.relativeToInitialImport = True
+        print(f"origin set to {lat}, {lon} (project origin)")
+
     osmp.commandLineMode = True
     osmp.mode = args.mode
     assert osmp.mode == args.mode
@@ -134,6 +143,118 @@ def import_osm(args):
     print("import finished")
 
 
+def proj_batch(points, proj_bin):
+    """WGS84 (lat, lon) -> EPSG:26911 (E, N) in one PROJ call."""
+    import subprocess
+
+    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", "%.5f"],
+        input=payload, capture_output=True, text=True, env=env, check=True,
+    ).stdout.strip().splitlines()
+    return [tuple(float(v) for v in line.split()) for line in out]
+
+
+def find_proj():
+    import glob
+    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
+    pattern = os.path.expanduser(
+        "~/.local/share/flatpak/app/org.qgis.qgis/x86_64/stable/*/files/bin/proj"
+    )
+    matches = sorted(glob.glob(pattern))
+    if not matches:
+        raise SystemExit("no PROJ binary; set PROJ_BIN")
+    return matches[0]
+
+
+def correct_frame(args):
+    """Map Blosm's spherical Transverse Mercator onto true EPSG:26911.
+
+    Blosm projects with a sphere of radius 6378137 about the origin, which is a
+    conformal but ~0.3 % oversized frame relative to ellipsoidal UTM. The
+    difference between the two over a 500 m block is an affine map to well under
+    a centimetre, so fit a similarity transform on four spread control points and
+    apply it to the imported objects. After this, the blockout uses the same
+    metric frame as the DEM and the origin documented in Geo/manifests.
+    """
+    import math
+
+    import bpy
+    from mathutils import Matrix
+    from util.transverse_mercator import TransverseMercator
+
+    lat, lon = (float(v) for v in args.origin.split(","))
+    tm = TransverseMercator(lat=lat, lon=lon)
+
+    south, west, north, east = (float(v) for v in args.bounds.split(","))
+    corners = [(south, west), (south, east), (north, west), (north, east),
+               ((south + north) / 2, (west + east) / 2)]
+    blosm_xy = [tm.fromGeographic(la, lo)[:2] for la, lo in corners]
+    utm_xy = proj_batch(corners, find_proj())
+
+    # full affine fit (6 parameters): absorbs the small shear the spherical
+    # projection introduces across the block, which a similarity leaves as ~30 cm
+    n = len(blosm_xy)
+    bxy = [p for p in blosm_xy]
+    uxy = [q for q in utm_xy]
+    # solve [x, y, 1] * [a, b, c]^T = u  for both components, least squares
+    def solve(target):
+        ata = [[0.0] * 3 for _ in range(3)]
+        atb = [0.0] * 3
+        for (x, y), t in zip(bxy, target):
+            row = (x, y, 1.0)
+            for i in range(3):
+                for j in range(3):
+                    ata[i][j] += row[i] * row[j]
+                atb[i] += row[i] * t
+        # Gaussian elimination
+        m = [ata[i][:] + [atb[i]] for i in range(3)]
+        for i in range(3):
+            piv = max(range(i, 3), key=lambda r: abs(m[r][i]))
+            m[i], m[piv] = m[piv], m[i]
+            for r in range(3):
+                if r != i:
+                    f = m[r][i] / m[i][i]
+                    for c in range(i, 4):
+                        m[r][c] -= f * m[i][c]
+        return [m[i][3] / m[i][i] for i in range(3)]
+
+    a, b, tx = solve([q[0] for q in uxy])
+    d, e, ty = solve([q[1] for q in uxy])
+
+    residuals = [
+        math.hypot(a * p[0] + b * p[1] + tx - q[0], d * p[0] + e * p[1] + ty - q[1])
+        for p, q in zip(blosm_xy, utm_xy)
+    ]
+    print(f"frame correction affine: X = {a:.6f}*x {b:+.6f}*y {tx:+.2f}, "
+          f"Y = {d:+.6f}*x {e:.6f}*y {ty:+.2f}")
+    print(f"frame correction residual: max {max(residuals)*100:.1f} cm")
+
+    # local origin for the UTM frame: the project origin
+    origin_lat_lon = (lat, lon)
+    origin_xy = proj_batch([origin_lat_lon], find_proj())[0]
+
+    m = Matrix(((a, b, 0.0, tx - origin_xy[0]), (d, e, 0.0, ty - origin_xy[1]),
+                (0.0, 0.0, 1.0, 0.0), (0.0, 0.0, 0.0, 1.0)))
+    # only transform roots: children inherit their parent's transform, and
+    # transforming both would apply the correction twice
+    roots = [o for o in bpy.data.objects if o.parent is None]
+    for obj in roots:
+        obj.matrix_world = m @ obj.matrix_world
+    print(f"frame correction applied to {len(roots)} root objects "
+          f"({len(bpy.data.objects) - len(roots)} children inherit)")
+
+
 def report():
     """Per-object measurements in metres, so a scale bug is obvious in the log."""
     collections = {c.name: c for c in bpy.data.collections}
@@ -167,6 +288,8 @@ def main():
     enable_blosm()
     configure_preferences(args)
     import_osm(args)
+    if args.origin:
+        correct_frame(args)
     report()
 
     out = os.path.abspath(args.out)

Tools/terrain_from_dem.py

@@ -0,0 +1,220 @@
+"""Build a terrain mesh from the committed 3DEP DEM, in project-local metres.
+
+Run inside Blender (background or GUI):
+
+    B=/var/lib/flatpak/app/org.blender.Blender/x86_64/stable/active/files
+    LD_LIBRARY_PATH="$B/blender/lib" "$B/blender/blender" -b --factory-startup \
+      --python Tools/terrain_from_dem.py -- \
+      --dem Geo/raw/dem-3dep-1m-26911-2026-11-06.tif \
+      --out Art/blender/terrain/terrain-2026-11-06.blend
+
+Coordinates follow `Geo/manifests/local-origin.md`:
+
+* X = UTM easting − 376 515.389 (east positive)
+* Y = UTM northing − 3 774 224.083 (north positive)
+* Z = NAVD88 elevation − 115.911 (the DEM value at the origin)
+
+The DEM is already EPSG:26911, so no reprojection happens here: pixel centres are
+UTM, and only the origin subtraction is applied. Negative X is west, negative Y
+is south — the whole block sits at negative X.
+
+The mesh is a regular grid at `--step` metres (default 2, which still resolves
+the ~2 % road crown over the ~21 m carriageway). Nodata cells are skipped by
+dropping faces that touch them, so no fake geometry is invented at the edges.
+"""
+
+import argparse
+import os
+import sys
+
+import subprocess
+
+import bpy
+
+sys.path.insert(0, os.path.join(os.path.dirname(os.path.abspath(__file__)), "..", "Geo", "scripts"))
+import dem_stats  # noqa: E402
+import osm_common as oc  # noqa: E402
+
+ORIGIN_E = 376515.389
+ORIGIN_N = 3774224.083
+ORIGIN_Z = 115.911
+
+
+def parse_args():
+    argv = sys.argv
+    argv = argv[argv.index("--") + 1:] if "--" in argv else []
+    ap = argparse.ArgumentParser()
+    ap.add_argument("--dem", required=True)
+    ap.add_argument("--out", required=True)
+    ap.add_argument("--step", type=float, default=2.0, help="grid spacing in metres")
+    ap.add_argument("--xmin", type=float, default=-330.0)
+    ap.add_argument("--xmax", type=float, default=50.0)
+    ap.add_argument("--ymin", type=float, default=-100.0)
+    ap.add_argument("--ymax", type=float, default=110.0)
+    ap.add_argument("--osm", default=None,
+                    help="OSM extract; building outlines are cut out so rooftops never become terrain")
+    return ap.parse_args(argv)
+
+
+def find_proj():
+    import glob
+    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; set PROJ_BIN")
+    return matches[0]
+
+
+def project_to_utm(points, proj_bin):
+    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", "%.4f"],
+                         input=payload, capture_output=True, text=True, env=env, check=True)
+    return [tuple(float(v) for v in line.split()) for line in out.stdout.strip().splitlines()]
+
+
+def building_polygons(osm_path):
+    """Building footprints as local X/Y rings, via PROJ into EPSG:26911."""
+    nodes, ways = oc.load(osm_path)
+    rings = []
+    for _, tags, refs in ways:
+        if "building" not in tags:
+            continue
+        pts = [nodes[r] for r in refs if r in nodes]
+        if len(pts) < 3:
+            continue
+        rings.append(pts)
+    if not rings:
+        return []
+    flat = [p for ring in rings for p in pts_flat(ring)]
+    utm = project_to_utm(flat, find_proj())
+    out, i = [], 0
+    for ring in rings:
+        n = len(ring)
+        out.append([(utm[j][0] - ORIGIN_E, utm[j][1] - ORIGIN_N) for j in range(i, i + n)])
+        i += n
+    return out
+
+
+def pts_flat(ring):
+    return ring
+
+
+def inside(x, y, ring):
+    """Ray casting point-in-polygon."""
+    hit = False
+    n = len(ring)
+    for i in range(n):
+        x1, y1 = ring[i]
+        x2, y2 = ring[(i + 1) % n]
+        if (y1 > y) != (y2 > y):
+            xin = (x2 - x1) * (y - y1) / (y2 - y1 + 1e-12) + x1
+            if x < xin:
+                hit = not hit
+    return hit
+
+
+def main():
+    args = parse_args()
+    entries, width, height, values = dem_stats.read_tiff(args.dem)
+    scale = entries[33550]
+    tie = entries[33922]
+    e0 = tie[3] - tie[0] * scale[0]
+    n0 = tie[4] + tie[1] * scale[1]
+    sx, sy = scale[0], scale[1]
+    nodata = None
+    if 42113 in entries:
+        try:
+            nodata = float(entries[42113])
+        except (TypeError, ValueError):
+            nodata = None
+
+    step = args.step
+    cols = int((args.xmax - args.xmin) / step) + 1
+    rows = int((args.ymax - args.ymin) / step) + 1
+    print(f"grid {cols} x {rows} at {step} m over "
+          f"X {args.xmin}..{args.xmax}, Y {args.ymin}..{args.ymax}")
+
+    polys = []
+    if args.osm:
+        polys = building_polygons(args.osm)
+        print(f"masking {len(polys)} building footprints")
+
+    verts = []
+    ok = []
+    for r in range(rows):
+        y = args.ymin + r * step
+        for c in range(cols):
+            x = args.xmin + c * step
+            e, n = x + ORIGIN_E, y + ORIGIN_N
+            col = (e - e0) / sx
+            row = (n0 - n) / sy
+            ci, ri = int(col), int(row)
+            if ci < 0 or ri < 0 or ci >= width or ri >= height:
+                verts.append((x, y, 0.0))
+                ok.append(False)
+                continue
+            if polys and any(inside(x, y, ring) for ring in polys):
+                verts.append((x, y, 0.0))
+                ok.append(False)
+                continue
+            z = values[ri * width + ci]
+            if nodata is not None and z == nodata:
+                verts.append((x, y, 0.0))
+                ok.append(False)
+                continue
+            verts.append((x, y, z - ORIGIN_Z))
+            ok.append(True)
+
+    faces = []
+    skipped = 0
+    for r in range(rows - 1):
+        for c in range(cols - 1):
+            i = r * cols + c
+            quad = (i, i + 1, i + cols + 1, i + cols)
+            if not all(ok[k] for k in quad):
+                skipped += 1
+                continue
+            faces.append(quad)
+
+    mesh = bpy.data.meshes.new("Terrain_DEM")
+    mesh.from_pydata(verts, [], faces)
+    mesh.update()
+    obj = bpy.data.objects.new("Terrain_DEM", mesh)
+    bpy.context.scene.collection.objects.link(obj)
+
+    zs = [v[2] for v, good in zip(verts, ok) if good]
+    print(f"vertices {len(verts)}, faces {len(faces)}, skipped cells {skipped}")
+    print(f"Z range {min(zs):+.2f} .. {max(zs):+.2f} m (relative to origin ground {ORIGIN_Z} m)")
+
+    # spot checks against the cross-section findings: gutters near +/-10.5 m
+    for name, x, y in (("origin (Highland)", 0.0, 0.0),
+                       ("Orange Dr end", -279.6, 0.0),
+                       ("hero frontage", -133.0, -14.5),
+                       ("north sidewalk", -133.0, 13.6)):
+        if not (args.xmin <= x <= args.xmax and args.ymin <= y <= args.ymax):
+            continue
+        c = int((x - args.xmin) / step)
+        r = int((y - args.ymin) / step)
+        i = r * cols + c
+        if ok[i]:
+            print(f"  {name:22} ({x:7.1f}, {y:6.1f}) -> Z {verts[i][2]:+.2f} m")
+
+    out = os.path.abspath(args.out)
+    os.makedirs(os.path.dirname(out), exist_ok=True)
+    bpy.ops.wm.save_as_mainfile(filepath=out)
+    print(f"saved {out}")
+
+
+if __name__ == "__main__":
+    main()