gig-line/Geo/tools/osm_common.py

#!/usr/bin/env python3
"""Shared geo helpers: parse .osm, stitch a street, project points, convert to Unreal.

Standard library only. Imported by every script in `Geo/tools/`.

The unit boundary is deliberate and lives here, once:

    WGS84 (lat, lon)  ->  EPSG:26911 (E, N metres)  ->  project local metres
                      ->  Unreal centimetres

`to_unreal` is the only place that multiplies by 100. Nothing else in the
repository may convert metres to centimetres, so a second silent scale cannot be
introduced by accident.
"""

from __future__ import annotations

import math
import os
import xml.etree.ElementTree as ET

EARTH_R = 6378137.0

# Project origin, EPSG:26911. See Geo/manifests/local-origin.md.
ORIGIN_E = 376515.389
ORIGIN_N = 3774224.083
ORIGIN_Z_NAVD88 = 115.911
METRES_TO_UNREAL = 100.0


def to_unreal_local(easting, northing, elevation=None):
    """EPSG:26911 metres -> Unreal centimetres in the project frame.

    X = east of the origin, Y = south of the origin (the FBX preset maps Blender
    +Y/north to Unreal -Y), Z = height above the origin ground. This is the single
    meters->centimetres conversion for the repository.
    """
    x_cm = (easting - ORIGIN_E) * METRES_TO_UNREAL
    y_cm = -(northing - ORIGIN_N) * METRES_TO_UNREAL
    if elevation is None:
        return x_cm, y_cm, None
    return x_cm, y_cm, (elevation - ORIGIN_Z_NAVD88) * METRES_TO_UNREAL


def to_blender_local(easting, northing, elevation=None):
    """EPSG:26911 metres -> Blender metres in the project frame (X east, Y north, Z up)."""
    x = easting - ORIGIN_E
    y = northing - ORIGIN_N
    z = None if elevation is None else elevation - ORIGIN_Z_NAVD88
    return x, y, z


def find_proj():
    """The PROJ binary: environment override, PATH, then the QGIS Flatpak."""
    import glob
    import os
    import shutil

    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 found; set PROJ_BIN or install pyproj")
    return matches[0]


def project_wgs84_to_utm(points, proj_bin=None):
    """WGS84 (lat, lon) -> EPSG:26911 (E, N) for a list, in one PROJ call."""
    import subprocess

    proj_bin = proj_bin or find_proj()
    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", "%.3f"],
        input=payload, capture_output=True, text=True, env=env, check=True,
    ).stdout.strip().splitlines()
    result = [tuple(float(v) for v in line.split()) for line in out]
    if len(result) != len(points):
        raise SystemExit(f"PROJ returned {len(result)} points for {len(points)} inputs")
    return result


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]