gig-line/Geo/scripts/osm_common.py

#!/usr/bin/env python3
"""Shared helpers for the Geo scripts: parse .osm, stitch a street, project points.

Standard library only. Imported by `block_report.py`, `street_profile.py`, and
`elevation_probe.py` when run from this directory.
"""

from __future__ import annotations

import math
import xml.etree.ElementTree as ET

EARTH_R = 6378137.0


def to_xy(lat: float, lon: float, lat0: float, lon0: float) -> tuple[float, float]:
    """Local planar metres, east/north, centred on (lat0, lon0)."""
    return (
        math.radians(lon - lon0) * EARTH_R * math.cos(math.radians(lat0)),
        math.radians(lat - lat0) * EARTH_R,
    )


def from_xy(x: float, y: float, lat0: float, lon0: float) -> tuple[float, float]:
    lat = lat0 + math.degrees(y / EARTH_R)
    lon = lon0 + math.degrees(x / (EARTH_R * math.cos(math.radians(lat0))))
    return lat, lon


def load(path: str):
    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 refs_by_name(ways, name):
    return {r for _, tags, refs in ways if tags.get("name") == name for r in refs}


def stitch(nodes, ways, name, lon_min, lon_max):
    """Ordered polyline for a named street, restricted to a longitude window."""
    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)

    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, prev = [], 0.0, 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: (distance_m, s_along_m, side) with side +1 = north."""
    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
        seg2 = dx * dx + dy * dy
        if seg2 == 0:
            continue
        t = max(0.0, min(1.0, ((px - ax) * dx + (py - ay) * dy) / seg2))
        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)
            best = (dist, lengths[i] + t * (lengths[i + 1] - lengths[i]), 1 if cross > 0 else -1)
    return best


def point_at(line, lengths, s, offset, lat0, lon0):
    """Point on the polyline at arc length s, pushed offset metres to the north side."""
    for i in range(len(line) - 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 = 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
            norm = math.hypot(dx, dy) or 1.0
            nx, ny = -dy / norm, dx / norm
            x = ax + t * dx + offset * nx
            y = ay + t * dy + offset * ny
            return from_xy(x, y, lat0, lon0)
    return line[-1]


def block_extent(nodes, ways, blvd, west, east):
    """Stitch and clip the boulevard between two cross streets.

    Returns (line, lat0, lon0, lengths, block_len) with `line` running west to east.
    """
    inter = {}
    blvd_refs = refs_by_name(ways, blvd)
    for label, street in (("west", west), ("east", east)):
        shared = [nodes[r] for r in refs_by_name(ways, street) & blvd_refs if r in nodes]
        if not shared:
            raise SystemExit(f"no intersection between {blvd!r} and {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, blvd, inter["west"], inter["east"])
    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)
    return line, lat0, lon0, lengths, lengths[-1]