zwitschi
386d4bee78
Geo/tools/ is now every geographic script, with Geo/README.md documenting each one; the Blender-side builders moved beside their outputs in Art/blender/*/ and import their helpers from Geo/tools/osm_common.py. That module gains to_unreal_local() and to_blender_local(): the only place in the repository that multiplies metres by 100, so a second silent scale cannot be introduced. block_height_report.py adds height provenance for the block — 3 of 15 buildings have no height data and render at a placeholder 6.4 m. asset_metadata.py writes and audits <asset>.geo.json sidecars; four now exist (hero-facade, terrain, street, blockout) and the audit recomputes each centre from the origin, failing anything that disagrees.
221 lines
7.3 KiB
Python
221 lines
7.3 KiB
Python
"""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")
|
||
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()
|