geo: measure block elevation and pick the DEM source
Probe 22 points with the USGS 3DEP point service: walkable surfaces span 115.65-119.59 m, a 3.94 m relief, so elevation data is required by the AGENTS rule. Grade falls ~3.3 m eastward over the 279 m block (1.2%) and the north sidewalk sits 0.4-0.7 m above the south. 3DEP 1 m is public domain and becomes the chosen source; the tile itself is still to be downloaded. Factored the shared OSM helpers into Geo/scripts/osm_common.py so the street profile and the elevation probe read the same geometry.
This commit is contained in:
zwitschi committed
1 parent
255c1a9e86
commit
65fb3759c9
9 files changed
+621
-152
No files matched your search
+28
-147
@@ -7,111 +7,23 @@ Usage:
|
||||
|
||||
Measures, in metres, from a plain .osm extract (standard library only):
|
||||
|
||||
* the boulevard centreline, stitched into one ordered line between the two
|
||||
boundary intersections;
|
||||
* the boulevard centreline, stitched and clipped 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;
|
||||
* each mapped crossing (`footway=crossing`) across the boulevard: its length,
|
||||
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.
|
||||
width tag. Derive those from the numbers printed here and record them as
|
||||
assumptions.
|
||||
"""
|
||||
|
||||
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
|
||||
import osm_common as oc
|
||||
|
||||
|
||||
def main():
|
||||
@@ -122,59 +34,32 @@ def main():
|
||||
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]
|
||||
nodes, ways = oc.load(args.osm)
|
||||
line, lat0, lon0, lengths, block_len = oc.block_extent(
|
||||
nodes, ways, args.blvd, args.west, args.east
|
||||
)
|
||||
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)
|
||||
pts = oc.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"]:
|
||||
if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
|
||||
continue
|
||||
a = oc.to_xy(pts[0][0], pts[0][1], lat0, lon0)
|
||||
b = oc.to_xy(pts[-1][0], pts[-1][1], lat0, lon0)
|
||||
length = ((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) ** 0.5
|
||||
mid = (sum(p[0] for p in pts) / len(pts), sum(p[1] for p in pts) / len(pts))
|
||||
dist, s, _ = oc.project(mid, line, lat0, lon0, lengths)
|
||||
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')}"
|
||||
f"(offset {dist:4.1f} m) markings={tags.get('crossing:markings', '-'):6} "
|
||||
f"signals={tags.get('crossing')}"
|
||||
)
|
||||
|
||||
print("\nsidewalk centrelines near the block (offset from road centreline):")
|
||||
@@ -182,20 +67,17 @@ def main():
|
||||
for wid, tags, refs in ways:
|
||||
if tags.get("footway") != "sidewalk":
|
||||
continue
|
||||
pts = way_points(nodes, refs)
|
||||
pts = oc.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"]:
|
||||
if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
|
||||
continue
|
||||
offsets = [project(p, line, lat0, lon0, lengths) for p in pts]
|
||||
offsets = [oc.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", "")))
|
||||
(south if offsets[0][2] < 0 else north).append((mean, along, wid, tags.get("name", "")))
|
||||
|
||||
for label, rows in (("south side", south), ("north side", north)):
|
||||
if not rows:
|
||||
@@ -210,16 +92,15 @@ def main():
|
||||
for wid, tags, refs in ways:
|
||||
if "building" not in tags:
|
||||
continue
|
||||
pts = way_points(nodes, refs)
|
||||
pts = oc.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"]:
|
||||
if min(p[1] for p in pts) > line[-1][1] or max(p[1] for p in pts) < line[0][1]:
|
||||
continue
|
||||
proj = [project(p, line, lat0, lon0, lengths) for p in pts]
|
||||
proj = [oc.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.append((nearest[0], nearest[2], wid, tags, max(along) - min(along)))
|
||||
rows.sort()
|
||||
for dist, side, wid, tags, frontage in rows:
|
||||
side_name = "S" if side < 0 else "N"
|
||||
|
||||
Reference in new issue
Block a user