zwitschi
255c1a9e86
New Geo/scripts/street_profile.py measures the 279.4 m block centreline, the three signalised crossings (20.9-25.9 m across), sidewalk centrelines at 12.6 m and 13.6 m, and per-building setbacks of 13.0-15.8 m. Derived sidewalk widths of ~4.0 m south and ~5.3 m north cross-check against the mapped centrelines. Geo/manifests/local-origin.md fixes EPSG:26911 and the Highland intersection node as origin, defines the axes and the x100 metre-to-centimetre rule, and leaves the projected pair blank because PROJ is not installed here. Records the lane-count discrepancy (4 tagged lanes over a 20.9 m carriageway).
239 lines
9.0 KiB
Python
Executable File
239 lines
9.0 KiB
Python
Executable File
#!/usr/bin/env python3
|
|
"""Street profile for a block: frontage setbacks, carriageway width, crossings.
|
|
|
|
Usage:
|
|
Geo/scripts/street_profile.py Geo/raw/osm-api-map-2026-11-06.osm \
|
|
--blvd "Hollywood Boulevard" --west "North Orange Drive" --east "North Highland Avenue"
|
|
|
|
Measures, in metres, from a plain .osm extract (standard library only):
|
|
|
|
* the boulevard centreline, stitched into one ordered line between the two
|
|
boundary intersections;
|
|
* each building's frontage setback: the distance from its nearest vertex to the
|
|
centreline, and which side of the road it is on;
|
|
* each mapped crossing (`footway=crossing`) across the boulevard: its length in
|
|
metres, which is the carriageway width at that point, and its position along
|
|
the block;
|
|
* the offset of the mapped sidewalk centrelines from the road centreline.
|
|
|
|
What it cannot measure: carriageway width where no crossing or kerb is mapped,
|
|
and sidewalk width, because OSM models sidewalks here as centre lines with no
|
|
width tag. Treat those as recorded assumptions, not measurements.
|
|
"""
|
|
|
|
import argparse
|
|
import math
|
|
import xml.etree.ElementTree as ET
|
|
|
|
EARTH_R = 6378137.0
|
|
|
|
|
|
def to_xy(lat, lon, lat0, lon0):
|
|
return (
|
|
math.radians(lon - lon0) * EARTH_R * math.cos(math.radians(lat0)),
|
|
math.radians(lat - lat0) * EARTH_R,
|
|
)
|
|
|
|
|
|
def load(path):
|
|
root = ET.parse(path).getroot()
|
|
nodes = {n.get("id"): (float(n.get("lat")), float(n.get("lon"))) for n in root.findall("node")}
|
|
ways = []
|
|
for w in root.findall("way"):
|
|
tags = {c.get("k"): c.get("v") for c in w.findall("tag")}
|
|
ways.append((w.get("id"), tags, [nd.get("ref") for nd in w.findall("nd")]))
|
|
return nodes, ways
|
|
|
|
|
|
def way_points(nodes, refs):
|
|
return [nodes[r] for r in refs if r in nodes]
|
|
|
|
|
|
def stitch(nodes, ways, name, lon_min, lon_max):
|
|
"""Ordered polyline for a named street between two longitudes."""
|
|
segments = []
|
|
for _, tags, refs in ways:
|
|
if tags.get("name") != name:
|
|
continue
|
|
pts = way_points(nodes, refs)
|
|
if not pts:
|
|
continue
|
|
if max(p[1] for p in pts) < lon_min - 1e-6 or min(p[1] for p in pts) > lon_max + 1e-6:
|
|
continue
|
|
segments.append(pts)
|
|
if not segments:
|
|
raise SystemExit(f"no geometry for {name!r} inside the block")
|
|
|
|
merged = []
|
|
for pts in segments:
|
|
if not merged:
|
|
merged = list(pts)
|
|
continue
|
|
if pts[0] == merged[-1]:
|
|
merged.extend(pts[1:])
|
|
elif pts[-1] == merged[-1]:
|
|
merged.extend(reversed(pts[:-1]))
|
|
elif pts[0] == merged[0]:
|
|
merged = list(reversed(pts)) + merged[1:]
|
|
elif pts[-1] == merged[0]:
|
|
merged = pts + merged[1:]
|
|
return merged
|
|
|
|
|
|
def cumulative_lengths(points, lat0, lon0):
|
|
xs = []
|
|
total = 0.0
|
|
prev = None
|
|
for lat, lon in points:
|
|
p = to_xy(lat, lon, lat0, lon0)
|
|
if prev is not None:
|
|
total += math.dist(prev, p)
|
|
xs.append(total)
|
|
prev = p
|
|
return xs
|
|
|
|
|
|
def project(point, line, lat0, lon0, lengths):
|
|
"""Nearest point on the polyline: returns (distance_m, s_along_m, side)."""
|
|
px, py = to_xy(point[0], point[1], lat0, lon0)
|
|
best = (float("inf"), 0.0, 0)
|
|
for i in range(len(line) - 1):
|
|
ax, ay = to_xy(line[i][0], line[i][1], lat0, lon0)
|
|
bx, by = to_xy(line[i + 1][0], line[i + 1][1], lat0, lon0)
|
|
dx, dy = bx - ax, by - ay
|
|
seg_len2 = dx * dx + dy * dy
|
|
if seg_len2 == 0:
|
|
continue
|
|
t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / seg_len2))
|
|
cx, cy = ax + t * dx, ay + t * dy
|
|
dist = math.hypot(px - cx, py - cy)
|
|
if dist < best[0]:
|
|
cross = dx * (py - ay) - dy * (px - ax)
|
|
side = 1 if cross > 0 else -1
|
|
best = (dist, lengths[i] + t * (lengths[i + 1] - lengths[i]), side)
|
|
return best
|
|
|
|
|
|
def main():
|
|
ap = argparse.ArgumentParser()
|
|
ap.add_argument("osm")
|
|
ap.add_argument("--blvd", default="Hollywood Boulevard")
|
|
ap.add_argument("--west", default="North Orange Drive")
|
|
ap.add_argument("--east", default="North Highland Avenue")
|
|
args = ap.parse_args()
|
|
|
|
nodes, ways = load(args.osm)
|
|
|
|
inter = {}
|
|
for label, street in (("west", args.west), ("east", args.east)):
|
|
refs = {r for _, tags, rs in ways if tags.get("name") == street for r in rs}
|
|
pts = way_points(nodes, refs)
|
|
blvd = {r for _, tags, rs in ways if tags.get("name") == args.blvd for r in rs}
|
|
shared = [nodes[r] for r in refs & blvd if r in nodes]
|
|
if not shared:
|
|
raise SystemExit(f"no intersection with {street!r}")
|
|
inter[label] = min(p[1] for p in shared) if label == "west" else max(p[1] for p in shared)
|
|
|
|
line = stitch(nodes, ways, args.blvd, inter["west"], inter["east"])
|
|
|
|
# orient west -> east, then clip to the block between the two intersections
|
|
if line[0][1] > line[-1][1]:
|
|
line = list(reversed(line))
|
|
west_pt = min(line, key=lambda p: abs(p[1] - inter["west"]))
|
|
east_pt = min(line, key=lambda p: abs(p[1] - inter["east"]))
|
|
i0, i1 = line.index(west_pt), line.index(east_pt)
|
|
if i0 > i1:
|
|
i0, i1 = i1, i0
|
|
line = line[i0 : i1 + 1]
|
|
|
|
lat0 = sum(p[0] for p in line) / len(line)
|
|
lon0 = sum(p[1] for p in line) / len(line)
|
|
lengths = cumulative_lengths(line, lat0, lon0)
|
|
block_len = lengths[-1]
|
|
print(f"centreline: {len(line)} vertices, {block_len:.1f} m between the boundary intersections\n")
|
|
|
|
print("crossings across the boulevard (length = carriageway width there):")
|
|
for wid, tags, refs in ways:
|
|
if tags.get("footway") != "crossing":
|
|
continue
|
|
pts = way_points(nodes, refs)
|
|
if not pts:
|
|
continue
|
|
lats = [p[0] for p in pts]
|
|
lons = [p[1] for p in pts]
|
|
# must span the roadway: the along-road extent small, the across extent large
|
|
a = to_xy(pts[0][0], pts[0][1], lat0, lon0)
|
|
b = to_xy(pts[-1][0], pts[-1][1], lat0, lon0)
|
|
length = math.dist(a, b)
|
|
mid = (sum(lats) / len(lats), sum(lons) / len(lons))
|
|
dist, s, side = project(mid, line, lat0, lon0, lengths)
|
|
if min(lons) > inter["east"] or max(lons) < inter["west"]:
|
|
continue
|
|
if not (0.0 <= s <= block_len) or dist > 4.0:
|
|
continue
|
|
marks = tags.get("crossing:markings", "?")
|
|
print(
|
|
f" way {wid:>10} {length:6.1f} m across, at {s:6.1f} m along the block "
|
|
f"(offset {dist:4.1f} m) markings={marks:6} signals={tags.get('crossing')}"
|
|
)
|
|
|
|
print("\nsidewalk centrelines near the block (offset from road centreline):")
|
|
south, north = [], []
|
|
for wid, tags, refs in ways:
|
|
if tags.get("footway") != "sidewalk":
|
|
continue
|
|
pts = way_points(nodes, refs)
|
|
if not pts:
|
|
continue
|
|
if min(p[1] for p in pts) > inter["east"] or max(p[1] for p in pts) < inter["west"]:
|
|
continue
|
|
offsets = [project(p, line, lat0, lon0, lengths) for p in pts]
|
|
mean = sum(o[0] for o in offsets) / len(offsets)
|
|
side = offsets[0][2]
|
|
along = max(o[1] for o in offsets) - min(o[1] for o in offsets)
|
|
# keep only long lines that run roughly parallel within the block
|
|
if not (6.0 <= mean <= 20.0) or along < 30.0:
|
|
continue
|
|
target = south if side < 0 else north
|
|
target.append((mean, along, wid, tags.get("name", "")))
|
|
|
|
for label, rows in (("south side", south), ("north side", north)):
|
|
if not rows:
|
|
print(f" {label}: none mapped")
|
|
continue
|
|
rows.sort()
|
|
mean = sum(r[0] for r in rows) / len(rows)
|
|
print(f" {label}: {len(rows)} line(s), mean offset {mean:.1f} m, range {rows[0][0]:.1f}..{rows[-1][0]:.1f} m")
|
|
|
|
print("\nbuilding frontage setbacks (nearest vertex to centreline):")
|
|
rows = []
|
|
for wid, tags, refs in ways:
|
|
if "building" not in tags:
|
|
continue
|
|
pts = way_points(nodes, refs)
|
|
if not pts:
|
|
continue
|
|
if min(p[1] for p in pts) > inter["east"] or max(p[1] for p in pts) < inter["west"]:
|
|
continue
|
|
proj = [project(p, line, lat0, lon0, lengths) for p in pts]
|
|
nearest = min(proj, key=lambda o: o[0])
|
|
along = [o[1] for o in proj]
|
|
frontage = max(along) - min(along)
|
|
rows.append((nearest[0], nearest[2], wid, tags, frontage))
|
|
rows.sort()
|
|
for dist, side, wid, tags, frontage in rows:
|
|
side_name = "S" if side < 0 else "N"
|
|
print(
|
|
f" {side_name} way {wid:>10} {tags.get('addr:housenumber', '-'):>11} "
|
|
f"setback {dist:5.1f} m, extent {frontage:5.1f} m {tags.get('name', '')[:34]}"
|
|
)
|
|
|
|
print("\nassumptions to record, not measurements:")
|
|
print(" - lane width is not tagged; the crossing lengths above are the best carriageway")
|
|
print(" proxy, and crossing endpoints sit within about 1 m of the kerb line")
|
|
print(" - sidewalk width is not tagged; it is setback minus half carriageway")
|
|
|
|
|
|
if __name__ == "__main__":
|
|
main()
|