zwitschi
d42e39faa9
Tools/street_mesh.py builds the corridor in the project frame from the measured profile: carriageway -10.5..+10.5 m, kerbs 0.15 m (0.13 m as built), sidewalks 4.0 m south and 5.3 m north, all elevations sampled from the DEM. Kerbs flush at the three signalised crossings. Crown at mid-block reproduces the cross-sections to the centimetre.
270 lines
9.7 KiB
Python
270 lines
9.7 KiB
Python
"""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()
|