gig-line/Tools/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", "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()