zwitschi
e75e5500a7
Merged PR #21: 0.1.1 collision patch (flush sidewalk edge, graded ramps) Co-authored-by: Zwitschi <zwitschi82@gmail.com>
381 lines
15 KiB
Python
381 lines
15 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;
|
||
* 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=None,
|
||
help="override the ramp run; default derives it from --max-ramp-grade")
|
||
ap.add_argument("--max-ramp-grade", type=float, default=0.08,
|
||
help="maximum crossing-ramp slope (0.08 = 8%%)")
|
||
args = ap.parse_args(argv)
|
||
args.ramp_blend = ramp_blend(args.kerb, args.max_ramp_grade, args.ramp_blend)
|
||
return args
|
||
|
||
|
||
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.
|
||
|
||
`blend` is the horizontal run over which the kerb drops from `base` to the
|
||
0.02 m flush height. It is derived from a target maximum grade (see
|
||
`ramp_blend`) so the crossing ramps stay within a real ramp's slope instead
|
||
of the ~8.7 % the old fixed 1.5 m run produced.
|
||
"""
|
||
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 ramp_blend(kerb, max_grade, override=None):
|
||
"""Horizontal run needed for the kerb ramp to stay at or under max_grade.
|
||
|
||
0.15 m kerb at an 8 % grade needs 1.63 m of run; the previous fixed 1.5 m
|
||
gave 8.7 %, and the walk pre-check measured the transition as a 10-13 cm
|
||
step over 1.5 m. When `override` is given it wins (kept for reproducing old
|
||
meshes), otherwise the run is computed from the live kerb height so it stays
|
||
correct if the kerb height assumption changes.
|
||
"""
|
||
if override is not None:
|
||
return max(0.0, override)
|
||
if max_grade <= 0.0:
|
||
return 0.0
|
||
return max(0.0, (kerb - 0.02)) / max_grade
|
||
|
||
|
||
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 sidewalk_z(r, off, inner, outer, x, y):
|
||
"""Sidewalk surface: full kerb upstand at the road-side (inner) edge,
|
||
tapering to zero at the building-side (outer) edge so the walk surface
|
||
meets the terrain flush.
|
||
|
||
Before this taper the outer lip stood the full kerb height above the DEM,
|
||
which the walkability pre-check measured as a 15.8 cm step where a sample
|
||
crosses from Terrain_DEM onto Sidewalk_South (with a matching lip on the
|
||
north side). A 0.15 m lip over a 4.0 m / 5.3 m strip is a ~3-4 % ramp:
|
||
walkable and visually negligible, and it removes the ledge entirely.
|
||
"""
|
||
z = dem_z(x, y)
|
||
base = z if z is not None else 0.0
|
||
h = kerb_height(rows[r][0], args.kerb, args.ramp_half_width, args.ramp_blend)
|
||
span = inner - outer
|
||
if abs(span) < 1e-9:
|
||
return base + h
|
||
t = min(1.0, max(0.0, (off - outer) / span))
|
||
return base + h * t
|
||
|
||
def south_z(r, off, x, y):
|
||
return sidewalk_z(r, off, SOUTH_GUTTER, SOUTH_SIDEWALK, x, y)
|
||
|
||
def north_z(r, off, x, y):
|
||
return sidewalk_z(r, off, NORTH_GUTTER, NORTH_SIDEWALK, 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()
|