geo: sample seven cross-sections from the DEM

Gutters sit at -10 m south and +10..12 m north, a ~2% crown between them, and
the envelope runs 0.6-1.0 m per section. Gutter spacing independently confirms
the 20.9 m carriageway measured from the crossing. The sharpest breaks are
0.09-0.12 m/m, so 3DEP does not resolve a kerb upstand: kerb height stays a
recorded assumption rather than a derived value.
This commit is contained in:
zwitschi committed 2026-11-06 23:02:50 +01:00
1 parent d97fcc7ad4
commit 6cded45e2a
5 files changed
+237 -3

No files matched your search

+1 -1
View File
@@ -13,7 +13,7 @@ OSM is layout data only. It is not authoritative for façade detail, height, or
| `Geo/manifests/` | One manifest per saved extract (see below) |
| `Geo/scripts/` | Repeatable fetch and conversion helpers |
Current state: the extract, both El Capitan footprints, the demo bounds, the building sets, the street profile, and the designated crossing are recorded in [`verification/`](verification/). The local CRS and origin are fixed in [`manifests/local-origin.md`](manifests/local-origin.md), with the projected easting/northing still to be computed. Elevation is measured and the 3DEP 1 m DEM window is committed (EPSG:26911, NAVD88, 984 KB). The hero façade width remains open.
Current state: the extract, both El Capitan footprints, the demo bounds, the building sets, the street profile, and the designated crossing are recorded in [`verification/`](verification/). The local CRS and origin are fixed in [`manifests/local-origin.md`](manifests/local-origin.md), with the projected easting/northing still to be computed. Elevation is measured and the 3DEP 1 m DEM window is committed (EPSG:26911, NAVD88, 984 KB). The street cross-sections (crown ~2 %, gutters at ±10–11 m) are in `verification/`. The hero façade width remains open.
## Required manifest
+187
View File
@@ -0,0 +1,187 @@
#!/usr/bin/env python3
"""Cross-sections of the street from the local 3DEP DEM.
Usage:
Geo/scripts/street_section.py Geo/raw/osm-api-map-2026-11-06.osm \
Geo/raw/dem-3dep-1m-26911-2026-11-06.tif --step 0.5
For several positions along the block it samples the DEM on a line perpendicular
to the boulevard centreline, from 18 m south to 18 m north, and reports:
* the height profile and total cross-fall;
* where the surface breaks most sharply, which is the best available candidate
for the kerb line;
* whether a kerb step is visible in the DEM at all.
The OSM geometry is WGS84, the DEM is EPSG:26911, so points are batch-projected
through the PROJ binary that ships with the QGIS flatpak (the host has no PROJ).
Caveat: 3DEP is bare-earth lidar at 1 m. Kerbs are short vertical steps and are
usually smoothed or removed by the DEM production process. If no break shows at
the expected offset, the DEM simply does not resolve kerbs and kerb height has to
come from a stated assumption instead.
"""
import argparse
import glob
import os
import shutil
import subprocess
import sys
import dem_stats
import osm_common as oc
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
pattern = os.path.expanduser(
"~/.local/share/flatpak/app/org.qgis.qgis/x86_64/stable/*/files/bin/proj"
)
matches = sorted(glob.glob(pattern))
if matches:
return matches[0]
raise SystemExit(
"no PROJ binary found. Set PROJ_BIN, or install pyproj, or use the QGIS flatpak."
)
def project_batch(points, proj_bin):
"""WGS84 (lat, lon) -> EPSG:26911 (E, N) for a list of points, one PROJ call."""
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:.8f} {lat:.8f}" for lat, lon in points) + "\n"
result = subprocess.run(
[proj_bin, "+proj=utm", "+zone=11", "+datum=NAD83", "-f", "%.4f"],
input=payload,
capture_output=True,
text=True,
env=env,
check=True,
)
out = []
for line in result.stdout.strip().splitlines():
parts = line.split()
out.append((float(parts[0]), float(parts[1])))
if len(out) != len(points):
raise SystemExit(f"PROJ returned {len(out)} points for {len(points)} inputs")
return out
class Dem:
def __init__(self, path):
entries, width, height, values = dem_stats.read_tiff(path)
self.w, self.h, self.v = width, height, values
scale = entries[33550]
tie = 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):
"""Bilinear sample; returns None outside the raster or on nodata."""
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)):
value = self.v[r * self.w + c]
if self.nodata is not None and value == self.nodata:
return None
p.append(value)
return (
p[0] * (1 - fx) * (1 - fy)
+ p[1] * fx * (1 - fy)
+ p[2] * (1 - fx) * fy
+ p[3] * fx * fy
)
def main():
ap = argparse.ArgumentParser()
ap.add_argument("osm")
ap.add_argument("dem")
ap.add_argument("--blvd", default="Hollywood Boulevard")
ap.add_argument("--west", default="North Orange Drive")
ap.add_argument("--east", default="North Highland Avenue")
ap.add_argument("--step", type=float, default=0.5, help="sampling step in metres")
ap.add_argument("--reach", type=float, default=18.0, help="sample +/- this many metres")
ap.add_argument("--section", type=float, action="append", default=[],
help="position along the block in metres; repeatable. Default: every 40 m")
args = ap.parse_args()
nodes, ways = oc.load(args.osm)
line, lat0, lon0, lengths, block_len = oc.block_extent(
nodes, ways, args.blvd, args.west, args.east
)
dem = Dem(args.dem)
proj_bin = find_proj()
sections = args.section or [40.0 * i for i in range(8) if 40.0 * i <= block_len]
offsets = []
x = -args.reach
while x <= args.reach + 1e-9:
offsets.append(round(x, 2))
x += args.step
# build every sample point, project them in one batch
grid = []
for s in sections:
row = []
for off in offsets:
row.append(oc.point_at(line, lengths, s, off, lat0, lon0))
grid.append(row)
flat = [p for row in grid for p in row]
projected = project_batch(flat, proj_bin)
print(f"block {block_len:.1f} m, {len(sections)} sections, "
f"offsets {offsets[0]:+.1f}..{offsets[-1]:+.1f} m at {args.step} m\n")
index = 0
for s, row in zip(sections, grid):
heights = []
for _ in row:
e, n = projected[index]
index += 1
heights.append(dem.at(e, n))
valid = [(o, h) for o, h in zip(offsets, heights) if h is not None]
if not valid:
print(f"s={s:6.1f} m: no DEM data")
continue
hs = [h for _, h in valid]
high_off, high_h = max(valid, key=lambda t: t[1])
low_off, low_h = min(valid, key=lambda t: t[1])
cross_fall = high_h - low_h
# sharpest break in the profile = candidate kerb
breaks = []
for (o1, h1), (o2, h2) in zip(valid, valid[1:]):
slope = abs(h2 - h1) / max(1e-9, abs(o2 - o1))
breaks.append((slope, (o1 + o2) / 2, h2 - h1))
breaks.sort(reverse=True)
top = breaks[:3]
print(f"s={s:6.1f} m crown {high_h:.2f} m at {high_off:+.1f} "
f"low {low_h:.2f} m at {low_off:+.1f} cross-fall {cross_fall:.2f} m")
print(" profile (2 m): " + " ".join(
f"{o:+.0f}:{h:.2f}" for o, h in valid if abs(o % 2) < 1e-6))
print(" sharpest breaks: " + ", ".join(
f"{slope:.2f} m/m at {off:+.1f} m ({dh:+.2f} m)" for slope, off, dh in top))
print()
if __name__ == "__main__":
sys.exit(main())
@@ -0,0 +1,47 @@
# Cross-sections — crowned road, gutters, sidewalk shelves
Sampled 2026-11-06 from the local 3DEP DEM (`Geo/raw/dem-3dep-1m-26911-2026-11-06.tif`) with `Geo/scripts/street_section.py`. Seven sections at 40 m spacing, offsets −18 m to +18 m from the boulevard centreline (negative = south), 1 m sampling step.
## What the profile shows
Every section has the same shape:
- a local low at about **−10 m** (south gutter),
- a rise to a flat-ish crown around **+2 to +8 m**,
- a second low at about **+10 to +12 m** (north gutter),
- then the north side climbing to the building line.
Those gutter positions bracket a carriageway of roughly **20–22 m**, which independently matches the 20.9 m measured from the mid-block crossing in [`2026-11-06-street-profile.md`](2026-11-06-street-profile.md). Two unrelated measurements agreeing is the strongest evidence in the geo work so far.
## Section table
| Along block | Lowest point | Road crown | Cross-fall across the road | Gutter candidates |
| ---: | ---: | ---: | ---: | --- |
| 0.0 m | 118.78 @ −10 | ~119.41 @ 0 | ~0.6 m / 2 % | −10, +10 |
| 40.0 m | 118.43 @ −11 | ~118.94 @ 0 | ~0.5 m / 2 % | −11, +10.5 |
| 80.0 m | 117.84 @ −10 | ~118.33 @ 0 | ~0.5 m / 2 % | −10, +10.5–11.5 |
| 120.0 m | 117.22 @ −10 | ~117.73 @ 0 | ~0.5 m / 2 % | −10, +11.5 |
| 160.0 m | 116.55 @ −10 | ~117.04 @ 0 | ~0.5 m / 2 % | −10, +11.5 |
| 200.0 m | 115.96 @ −10 | ~116.37 @ 0 | ~0.4 m / 2 % | −10, +10.5–11.5 |
| 240.0 m | 115.63 @ −10 | ~116.02 @ 0 | ~0.4 m / 2 % | −10, +10.5–11.5 |
Total envelope from the lowest south gutter to the highest north point is 0.6–1.0 m per section, depending on how far the north sidewalk has risen.
## Kerbs are not in the DEM
The sharpest breaks found at the expected kerb offsets are **0.09–0.12 m/m** slopes — that is a gutter transition spread over a metre, not a kerb upstand. 3DEP is bare-earth and smooths a 0.15 m vertical step at 1 m sampling, so the DEM cannot give kerb height.
Consequence: kerb height and shape are an **assumption** until measured or observed in the reference photography. Write the assumption down where it is used rather than burying it in the mesh. Do not infer a kerb from the DEM's gradient.
## What the demo ground mesh must do
- Crown the road at roughly 2 % from each gutter to the centre, not a flat plane and not a constant tilt.
- Step down into the gutters at about ±10.5 m from the centreline, then run the sidewalks at their own grade.
- 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.
## Open
- Kerb height, kerb radius, and ramp slope: assumption until surveyed or measured from photographs.
- The 200 m section peaks at +14 m rather than +18 m, which may be a tree canopy or a building edge in the DEM. Do not treat it as a road feature.
- Drainage inlets, driveway aprons, and pavement joints are not in this data.
+1 -1
View File
@@ -52,5 +52,5 @@ Window statistics: min 110.46 m, max 127.73 m, spread 17.28 m — building tops
## Next
1. Sample the kerb line and the crossing ramps specifically — the grid above is a street-level average, not a kerb profile.
1. ~~Sample the kerb line and the crossing ramps~~ — done: seven cross-sections in [`2026-11-06-cross-sections.md`](2026-11-06-cross-sections.md) locate the gutters at ±10–11 m and show a ~2 % crown. They also show the DEM does **not** resolve a kerb upstand, so kerb height stays an assumption.
2. Build the ground mesh from the DEM along the street, not from a flat plane, and assert the 3.94 m fall survives the Blender→Unreal path.