gig-line/Art/blender/terrain/terrain_from_dem.py

"""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", "tools"))
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")
    ap.add_argument("--corridor", type=float, default=10.7,
                    help="half-width in m of the Hollywood Blvd corridor to cut out so the "
                         "street meshes (road, kerbs, sidewalks) are not buried under the 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 centreline_polyline(osm_path):
    """Hollywood Blvd centreline as a local-frame polyline (Highland -> Orange)."""
    nodes, ways = oc.load(osm_path)
    line, _, _, _, _ = oc.block_extent(
        nodes, ways, "Hollywood Boulevard", "North Orange Drive", "North Highland Avenue")
    utm = project_to_utm(line, find_proj())
    return [(e - ORIGIN_E, n - ORIGIN_N) for e, n in utm]


def dist_to_polyline(x, y, pts):
    best = float("inf")
    for (ax, ay), (bx, by) in zip(pts, pts[1:]):
        dx, dy = bx - ax, by - ay
        l2 = dx * dx + dy * dy
        if l2 == 0.0:
            continue
        t = max(0.0, min(1.0, ((x - ax) * dx + (y - ay) * dy) / l2))
        px, py = ax + t * dx, ay + t * dy
        d = (x - px) ** 2 + (y - py) ** 2
        if d < best:
            best = d
    return best ** 0.5


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 = []
    corridor = None
    if args.osm:
        polys = building_polygons(args.osm)
        print(f"masking {len(polys)} building footprints")
        corridor = centreline_polyline(args.osm)
        print(f"cutting {args.corridor} m corridor around the Hollywood Blvd centreline")

    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
            if corridor is not None and dist_to_polyline(x, y, corridor) <= args.corridor:
                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)
    # *.blend is LFS-lockable, so a fresh checkout marks it read-only. Make it
    # writable again or the save fails and downstream steps build from a stale
    # file — which is worse than failing here.
    if os.path.exists(out) and not os.access(out, os.W_OK):
        os.chmod(out, 0o644)
    result = bpy.ops.wm.save_as_mainfile(filepath=out)
    if result != {"FINISHED"}:
        raise RuntimeError(f"could not save {out}: {result}")
    print(f"saved {out}")


if __name__ == "__main__":
    main()