feat: road, kerbs, and sidewalks on the DEM
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.
This commit is contained in:
zwitschi committed
1 parent
ece48649f5
commit
d42e39faa9
6 files changed
+300
-1
No files matched your search
@@ -13,6 +13,8 @@ Repeatable import/export and processing helpers. Anything that was done twice by
|
||||
## Helpers
|
||||
|
||||
- `check-lfs.sh` — repository guard: fails if `git-lfs` or the `lfs` filter is missing, or if an LFS-tracked file was committed as a plain blob.
|
||||
- `street_mesh.py` — road, kerbs, and sidewalks on the DEM, in the project frame; kerbs flush at the three crossings.
|
||||
- `terrain_from_dem.py` — terrain grid from the committed DEM, with building footprints cut out.
|
||||
- `blosm_blockout.py` — headless Blosm import of the committed `.osm` extract into a `.blend`, with a per-building scale report. Exact invocation in `Art/blender/README.md`.
|
||||
- `commons_fetch.py` — Wikimedia Commons reference fetcher. Search or exact `File:` titles, 1600 px thumbnails, provenance `sidecar.json` per folder. Writes only into the gitignored `Art/reference/`; reference images are never committed.
|
||||
|
||||
|
||||
@@ -0,0 +1,269 @@
|
||||
"""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()
|
||||
Reference in new issue
Block a user