"""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; * 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", "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 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 faces.append((i, i + 1, i + cols + 1, i + cols)) 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 + 2, i + 3, i + 1)) 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") 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) bpy.ops.wm.save_as_mainfile(filepath=out) print(f"saved {out}") if __name__ == "__main__": main()