From ece48649f5e8e4ef21cb86ce2fa79a32269f3ad2 Mon Sep 17 00:00:00 2001 From: Zwitschi Date: Fri, 6 Nov 2026 23:40:50 +0100 Subject: [PATCH] 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. --- Art/blender/README.md | 44 +++- .../blockout/blockout-2026-11-06.blend | 4 +- Art/blender/terrain/terrain-2026-11-06.blend | 3 + Geo/verification/2026-11-06-cross-sections.md | 4 + ROADMAP.md | 2 +- Tools/blosm_blockout.py | 123 ++++++++++ Tools/terrain_from_dem.py | 220 ++++++++++++++++++ 7 files changed, 391 insertions(+), 9 deletions(-) create mode 100644 Art/blender/terrain/terrain-2026-11-06.blend create mode 100644 Tools/terrain_from_dem.py diff --git a/Art/blender/README.md b/Art/blender/README.md index e54bfaa..096c982 100644 --- a/Art/blender/README.md +++ b/Art/blender/README.md @@ -29,18 +29,50 @@ Run time is about 36 s for this extent. Two headless-only workarounds are baked 1. Blosm refuses to import until an asset directory is set; the paid asset pack path is passed explicitly and is **not** committed. 2. In background mode the GUI never fills Blosm's default level table, and the add-on then takes a code path with a broken relative import. The script seeds one entry. +## Frame correction + +Blosm projects with a **sphere** of radius 6 378 137 m about the origin, not an ellipsoid. Against EPSG:26911 that is conformal but about 0.4 % oversized and rotated by the meridian convergence, which put the theatre ~2 m off. `Tools/blosm_blockout.py --origin` therefore fits a six-parameter affine from Blosm's frame into UTM on five spread control points and applies it to the imported roots. Fit residual: **0.0 cm**. After it, the blockout and the terrain share one metric frame. + +Two traps, both already handled in the script: transform only **root** objects (children inherit, and doing both applies the transform twice — that looked like a 2 m placement error), and fit an affine, not a similarity (a similarity leaves ~30 cm). + +## Terrain + +`terrain/terrain-2026-11-06.blend` — regular grid from the committed 3DEP DEM, built by `Tools/terrain_from_dem.py`. + +```bash +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 --osm Geo/raw/osm-api-map-2026-11-06.osm --out Art/blender/terrain/terrain-2026-11-06.blend +``` + +- 191 × 106 grid at 2 m over X −330…+50, Y −100…+110: 20 246 vertices, 13 640 faces. +- Z is NAVD88 minus 115.911 m, the DEM value at the origin, so the origin ground is Z 0. +- Building footprints from the OSM extract are **cut out** (6 310 cells), so rooftops never become terrain. 72 footprints masked. +- Z range −2.80…+8.19 m. The high values are the ground rising north of the boulevard, not rooftops: the maximum sits at the northern edge, Y +110. + +Spot checks against the measured cross-sections: + +| Point | Terrain Z | +| --- | ---: | +| Origin, Hollywood × Highland | 0.00 | +| Orange Drive end of the block | +3.39 | +| North sidewalk mid-block | +1.43 | +| South gutter, mid-block | +0.93 | +| Road centre, mid-block | +1.38 | + +That reproduces both the 3.4 m eastward fall and the ~2 % crown, so the DEM, the cross-sections, and this mesh agree. + ## Scale check (blockout, 2026-11-06) Blender works in metres, and the imported meshes match the OSM geometry, so the blockout is 1:1 before any Unreal export: -| Building | Measured in Blender | OSM geometry | Height source | +| Building | Blender dimensions | OSM way outline | Height source | | --- | --- | --- | --- | -| El Capitan Entertainment Centre | 25.68 × 54.88 × 14.50 m | 25.75 × 54.88 m | `height` tag 14.5 | -| El Capitan Theatre | 38.59 × 54.87 × 20.62 m | 38.59 × 54.87 m | `building:levels=6` × 3.2 m, plus roof | -| Madame Tussaud's | 42.41 × 57.20 × 15.90 m | 42.0 m frontage | `height` tag 15.9 | -| TCL Chinese Theatre | 43.75 × 76.19 × 18.10 m | 44.0 m frontage | `height` tag 18.1 | +| El Capitan Entertainment Centre | 25.72 × 54.47 × 14.50 m | 25.75 × 54.88 m | `height` tag 14.5 | +| El Capitan Theatre | 38.66 × 54.46 × 20.96 m | 38.59 × 54.87 m | `building:levels=6` × 3.2 m, plus roof | +| Madame Tussaud's | 42.48 × 56.77 × 15.90 m | 42.0 m frontage | `height` tag 15.9 | +| TCL Chinese Theatre | 43.78 × 75.90 × 18.10 m | 44.0 m frontage | `height` tag 18.1 | -The 0.07 m width difference on the hero is Blosm's own projection of the outline, not a units problem. Anything larger than a few centimetres would be a flag. +Widths agree to a few centimetres; depths differ by up to 0.4 m because Blosm assembles **multipolygon relations**, not the raw way I measured. Placement after the frame correction agrees with the UTM control points to 0.1–0.5 m per building — the residual is relation assembly and Blosm's edge straightening, not projection. Buildings **without** a `height` tag get levels × 3.2 m, which is a placeholder, not evidence. Every such height is a modelling decision to be replaced from the reference photography. diff --git a/Art/blender/blockout/blockout-2026-11-06.blend b/Art/blender/blockout/blockout-2026-11-06.blend index a87d37d..b9c1496 100644 --- a/Art/blender/blockout/blockout-2026-11-06.blend +++ b/Art/blender/blockout/blockout-2026-11-06.blend @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:6e50ca1a5c9e5c8513719d1dfed6ba92987221db115924dbd38e68ed19201d59 -size 250112 +oid sha256:2de637944c86dc656593ace2af0d3275f69174cc63133b7e80bbf81a0f55cc0c +size 250220 diff --git a/Art/blender/terrain/terrain-2026-11-06.blend b/Art/blender/terrain/terrain-2026-11-06.blend new file mode 100644 index 0000000..ef78caf --- /dev/null +++ b/Art/blender/terrain/terrain-2026-11-06.blend @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6f3522d20fd25df5e533cb821117b0237da9a38338b4dfb00e1ad391870100aa +size 470245 diff --git a/Geo/verification/2026-11-06-cross-sections.md b/Geo/verification/2026-11-06-cross-sections.md index 5171734..6d62f79 100644 --- a/Geo/verification/2026-11-06-cross-sections.md +++ b/Geo/verification/2026-11-06-cross-sections.md @@ -40,6 +40,10 @@ Consequence: kerb height and shape are an **assumption** until measured or obser - The north sidewalk sits higher than the road crown and keeps rising toward the building line, so the two sidewalks are not mirror images. - Crossing ramps must be modelled deliberately at the three crossings; nothing in the DEM supplies the local ramp geometry. +## Reproduced in the terrain mesh + +`Art/blender/terrain/terrain-2026-11-06.blend` samples the DEM on a 2 m grid and reproduces these numbers: origin Z 0.00, Orange Drive end +3.39, south gutter +0.93, road centre +1.38, north sidewalk +1.43 relative to the origin ground. The crown therefore survives from DEM to mesh. See `Art/blender/README.md`. + ## Open - Kerb height, kerb radius, and ramp slope: assumption until surveyed or measured from photographs. diff --git a/ROADMAP.md b/ROADMAP.md index adf656b..16ff4f0 100644 --- a/ROADMAP.md +++ b/ROADMAP.md @@ -42,7 +42,7 @@ Status detail in [Docs/setup.md](Docs/setup.md); workstation audit in [Docs/envi - [ ] Model hero façade for the El Capitan Entertainment Centre site. - [ ] Model surrounding block frontage. - [ ] Model simplified background buildings directly across Hollywood Boulevard. -- [ ] Build sidewalks, roads, and collision meshes. Ground must follow the DEM: 3.94 m fall along the block, ~2 % crown between the gutters at ±10–11 m from the centreline, north sidewalk higher than south. Kerb height is an unresolved assumption — the DEM does not resolve it (`Geo/verification/2026-11-06-cross-sections.md`). +- [-] Build sidewalks, roads, and collision meshes. Terrain mesh from the DEM exists (`Art/blender/terrain/terrain-2026-11-06.blend`); sidewalks, roads, kerbs and collision are not built. Ground must follow the DEM: 3.94 m fall along the block, ~2 % crown between the gutters at ±10–11 m from the centreline, north sidewalk higher than south. Kerb height is an unresolved assumption — the DEM does not resolve it (`Geo/verification/2026-11-06-cross-sections.md`). - [ ] Materials and look development for all exported assets. ### Gameplay diff --git a/Tools/blosm_blockout.py b/Tools/blosm_blockout.py index df8a251..bf01ab8 100644 --- a/Tools/blosm_blockout.py +++ b/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) diff --git a/Tools/terrain_from_dem.py b/Tools/terrain_from_dem.py new file mode 100644 index 0000000..750776e --- /dev/null +++ b/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()