gig-line/Geo/scripts/street_section.py

#!/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())