gig-line/Art/blender/street/street_mesh.py

"""Build the street corridor — road, kerbs, sidewalks — on the DEM.

    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/street_mesh.py -- \
      --osm Geo/raw/osm-api-map-2026-11-06.osm \
      --dem Geo/raw/dem-3dep-1m-26911-2026-11-06.tif \
      --out Art/blender/street/street-2026-11-06.blend

Everything lands in the project frame from `Geo/manifests/local-origin.md`:
X east, Y north, Z up, metres, origin at Hollywood × Highland.

Geometry follows the measured profile rather than a guess:

* carriageway between gutters at ±10.5 m from the centreline (cross-sections
  put them at −10 and +10…+12; the crossing length says 20.9 m);
* kerb upstand **0.15 m**, which is an assumption — the DEM does not resolve a
  kerb (see `Geo/verification/2026-11-06-cross-sections.md`);
* kerbs are flushed to 0.02 m within ±2.5 m of each signalised crossing, with a
  1.5 m blend, so the crossings are actually usable on foot;
* each signalised crossing gets continental zebra markings — 0.5 m white bars
  with 0.5 m gaps across a 3 m band, laid 1 cm above the road surface — so the
  mid-block crossing (OSM way 385648319, s ≈ 130 m) reads as a designated
  crossing in game;
* sidewalk widths 4.0 m south and 5.3 m north, from the DEM-derived setback
  arithmetic;
* the DEM supplies every elevation, so the 3.4 m fall and the ~2 % crown come
  through unchanged.

The road surface is lifted 2 cm above the DEM so it does not z-fight the terrain
mesh. When the terrain is cut to the corridor later, that lift can go.
"""

import argparse
import glob
import math
import os
import shutil
import subprocess
import sys

import bpy

HERE = os.path.dirname(os.path.abspath(__file__))
sys.path.insert(0, os.path.join(HERE, "..", "..", "..", "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

SOUTH_GUTTER = -10.5
NORTH_GUTTER = 10.5
SOUTH_SIDEWALK = -14.5
NORTH_SIDEWALK = 15.8
CROSSINGS = (36.6, 130.1, 266.1)


def parse_args():
    argv = sys.argv
    argv = argv[argv.index("--") + 1:] if "--" in argv else []
    ap = argparse.ArgumentParser()
    ap.add_argument("--osm", required=True)
    ap.add_argument("--dem", required=True)
    ap.add_argument("--out", required=True)
    ap.add_argument("--along", type=float, default=1.0)
    ap.add_argument("--cross", type=float, default=0.75)
    ap.add_argument("--kerb", type=float, default=0.15)
    ap.add_argument("--ramp-half-width", type=float, default=2.5)
    ap.add_argument("--ramp-blend", type=float, default=1.5)
    return ap.parse_args(argv)


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
    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()]


class Dem:
    def __init__(self, path):
        entries, w, h, values = dem_stats.read_tiff(path)
        self.w, self.h, self.v = w, h, values
        scale, tie = entries[33550], 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):
        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)):
            val = self.v[r * self.w + c]
            if self.nodata is not None and val == self.nodata:
                return None
            p.append(val)
        return (p[0] * (1 - fx) * (1 - fy) + p[1] * fx * (1 - fy)
                + p[2] * (1 - fx) * fy + p[3] * fx * fy)


def kerb_height(s, base, half_width, blend):
    """Full kerb, except flushed near a crossing so people can step off."""
    for crossing in CROSSINGS:
        d = abs(s - crossing)
        if d <= half_width:
            return 0.02
        if d <= half_width + blend:
            t = (d - half_width) / blend
            return 0.02 + (base - 0.02) * t
    return base


def strip(rows, offsets, z_at, name, material):
    """Build a quad grid from per-row sampled vertices."""
    verts, faces = [], []
    cols = len(offsets)
    for r, (_s, base, normal) in enumerate(rows):
        for off in offsets:
            x = base[0] + off * normal[0]
            y = base[1] + off * normal[1]
            verts.append((x, y, z_at(r, off, x, y)))
    for r in range(len(rows) - 1):
        for c in range(cols - 1):
            i = r * cols + c
            # Winding faces up: the natural (i, i+1, i+cols+1, i+cols) order
            # produces downward normals (every surface faced the ground and was
            # invisible in engine until the terrain corridor cut exposed it).
            faces.append((i, i + cols, i + cols + 1, i + 1))
    mesh = bpy.data.meshes.new(name)
    mesh.from_pydata(verts, [], faces)
    mesh.update()
    if material:
        mesh.materials.append(material)
    obj = bpy.data.objects.new(name, mesh)
    bpy.context.scene.collection.objects.link(obj)
    return obj


def main():
    args = parse_args()
    nodes, ways = oc.load(args.osm)
    line, _, _, _, block_len = oc.block_extent(
        nodes, ways, "Hollywood Boulevard", "North Orange Drive", "North Highland Avenue")
    utm = project_to_utm(line, find_proj())
    local = [(e - ORIGIN_E, n - ORIGIN_N) for e, n in utm]

    # arc length + unit normals in the local frame
    acc, lengths = 0.0, [0.0]
    for a, b in zip(local, local[1:]):
        acc += math.hypot(b[0] - a[0], b[1] - a[1])
        lengths.append(acc)
    total = lengths[-1]
    print(f"centreline {total:.1f} m")

    def at(s):
        for i in range(len(local) - 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 = local[i]
                bx, by = local[i + 1]
                dx, dy = bx - ax, by - ay
                norm = math.hypot(dx, dy) or 1.0
                return ((ax + t * dx, ay + t * dy), (-dy / norm, dx / norm), (dx / norm, dy / norm))
        ax, ay = local[-1]
        return ((ax, ay), (0.0, 1.0), (1.0, 0.0))

    dem = Dem(args.dem)

    rows = []
    s = 0.0
    while s <= total + 1e-9:
        base, normal, tangent = at(s)
        rows.append((s, base, normal))
        s += args.along

    # materials
    asphalt = bpy.data.materials.new("Asphalt")
    asphalt.diffuse_color = (0.09, 0.09, 0.095, 1.0)
    concrete = bpy.data.materials.new("Concrete")
    concrete.diffuse_color = (0.55, 0.53, 0.50, 1.0)

    def dem_z(x, y):
        z = dem.at(x + ORIGIN_E, y + ORIGIN_N)
        return None if z is None else z - ORIGIN_Z

    road_offsets = [SOUTH_GUTTER + i * (NORTH_GUTTER - SOUTH_GUTTER) / 28 for i in range(29)]
    south_offsets = [SOUTH_SIDEWALK + i * (SOUTH_GUTTER - SOUTH_SIDEWALK) / 8 for i in range(9)]
    north_offsets = [NORTH_GUTTER + i * (NORTH_SIDEWALK - NORTH_GUTTER) / 10 for i in range(11)]

    def road_z(r, off, x, y):
        z = dem_z(x, y)
        return (z if z is not None else 0.0) + 0.02

    def south_z(r, off, x, y):
        z = dem_z(x, y)
        base = z if z is not None else 0.0
        return base + kerb_height(rows[r][0], args.kerb, args.ramp_half_width, args.ramp_blend)

    def north_z(r, off, x, y):
        return south_z(r, off, x, y)

    road = strip(rows, road_offsets, road_z, "Road_Hollywood", asphalt)
    sw_s = strip(rows, south_offsets, south_z, "Sidewalk_South", concrete)
    sw_n = strip(rows, north_offsets, north_z, "Sidewalk_North", concrete)

    # kerb faces: vertical quads joining road edge to sidewalk edge
    def kerb_faces(offset, name):
        verts, faces = [], []
        for idx, (s_val, base, normal) in enumerate(rows):
            x = base[0] + offset * normal[0]
            y = base[1] + offset * normal[1]
            z_road = road_z(idx, offset, x, y)
            h = kerb_height(s_val, args.kerb, args.ramp_half_width, args.ramp_blend)
            z_walk = (dem_z(x, y) or 0.0) + h
            verts.extend([(x, y, z_road), (x, y, z_walk)])
        for idx in range(len(rows) - 1):
            i = idx * 2
            faces.append((i, i + 1, i + 3, i + 2))
        mesh = bpy.data.meshes.new(name)
        mesh.from_pydata(verts, [], faces)
        mesh.update()
        mesh.materials.append(concrete)
        obj = bpy.data.objects.new(name, mesh)
        bpy.context.scene.collection.objects.link(obj)
        return obj

    kerb_faces(SOUTH_GUTTER, "Kerb_South")
    kerb_faces(NORTH_GUTTER, "Kerb_North")

    # Continental zebra markings: one band per signalised crossing. Bars run
    # along the road (traffic direction), side by side across the carriageway;
    # a pedestrian crosses perpendicular to them. Lifted 1 cm above the road
    # surface (which is itself 2 cm above the DEM) to avoid z-fighting.
    paint = bpy.data.materials.new("Road_Paint_White")
    paint.diffuse_color = (0.85, 0.85, 0.85, 1.0)

    # All three bands live in ONE object. Creating one object per crossing used
    # to produce Crossing_Markings, .001 and .002; site assembly appends by the
    # exact name, so only the first crossing ever reached the game.
    band_verts, band_faces = [], []

    def add_crossing_band(s_c):
        # Stripes are emitted as two triangles each: the road crown makes the
        # quad non-planar and the FBX SDK refuses to triangulate it on import
        # ("unable to triangulate mesh" — the mesh is omitted entirely).
        zone = 3.0        # band length along the road, m
        w, g = 0.5, 0.5   # stripe width / gap, m
        span = (SOUTH_GUTTER + 0.2, NORTH_GUTTER - 0.2)
        n = int(zone // (w + g))
        covered = n * w + (n - 1) * g   # actual stripe coverage, centred on s_c
        for k in range(n):
            s0 = s_c - covered / 2 + k * (w + g)
            base_idx = len(band_verts)  # sampled positions, then vertices
            pts = []
            for s_edge in (s0, s0 + w):
                base, normal, _ = at(s_edge)
                for off in span:
                    x = base[0] + off * normal[0]
                    y = base[1] + off * normal[1]
                    # +6 cm above the DEM = 4 cm above the road surface: the road
                    # interpolates the DEM along 1 m rows while the markings sample
                    # only their own corners, so a 1 cm lift left parts of the
                    # stripes sunk under the road and invisible in game.
                    pts.append((x, y, (dem_z(x, y) or 0.0) + 0.06))
            # pts order: (s0,off_s), (s0,off_n), (s1,off_s), (s1,off_n)
            band_verts.extend(pts)
            band_faces.append((base_idx + 0, base_idx + 3, base_idx + 1))
            band_faces.append((base_idx + 0, base_idx + 2, base_idx + 3))

    for s_c in CROSSINGS:
        add_crossing_band(s_c)

    mesh = bpy.data.meshes.new("Crossing_Markings")
    mesh.from_pydata(band_verts, [], band_faces)
    mesh.update()
    mesh.materials.append(paint)
    markings = bpy.data.objects.new("Crossing_Markings", mesh)
    bpy.context.scene.collection.objects.link(markings)
    print(f"crossing markings: {len(CROSSINGS)} bands, "
          f"{len(markings.data.vertices)} verts, {len(markings.data.polygons)} faces")

    for obj in (road, sw_s, sw_n):
        print(f"{obj.name:16} verts {len(obj.data.vertices):5} faces {len(obj.data.polygons):5}")

    # crown check on the built mesh
    mid = total / 2
    for off in (SOUTH_GUTTER, -5.0, 0.0, 5.0, NORTH_GUTTER):
        base, normal, _ = at(mid)
        x = base[0] + off * normal[0]
        y = base[1] + off * normal[1]
        zr = road_z(0, off, x, y)
        print(f"  crown at {off:+5.1f} m: road Z {zr:+.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()