feat(phase-12f): curved corridor DEM extraction

CHANGELOG.md

@@ -27,7 +27,23 @@ a documentation / housekeeping commit on `main`.
   `min_cable_elev_m`. `tests/test_forced_flyover.py` — 10 new tests
   (dataclass invariants, deficit math, evaluator penalty firing /
   silent / outside-corridor, GA-routes-up integration, coexistence
-  with `NoTowerZone`). Suite 184 → 199.
+  with `NoTowerZone`).
+
+- **Phase 12f — curved corridor DEM extraction**
+  (`extract_profile_along_route(dem, waypoints, kind="polyline"|"bezier")`,
+  `src/ropeway/dem.py`). The case studies showed urban routes always
+  exceed the straight-line geodesic (Línea Roja 1.5×, Medellín 1.3×,
+  Cablebús 2.2×). `routes.py` already built curved 2-D paths but they
+  were never wired into the DEM sampler. The new helper accepts a list
+  of waypoints (raw `(lon, lat)` tuples or `Waypoint` objects), builds
+  a polyline or Catmull-Rom Bezier route in local UTM, samples the
+  DEM along it, and returns a standard `TerrainProfile` keyed on
+  along-route arc length — drops straight into `optimize(profile_fn=…,
+  corridor_length=…)`. Two-waypoint calls collapse to the existing
+  `extract_profile_from_dem` behaviour. `tests/test_curved_corridor.py`
+  — 6 new tests (two-waypoint ≡ straight, intermediate waypoint
+  lengthens corridor, Bezier ≥ straight, returned profile drops into
+  optimizer, tuple/Waypoint interchangeable, fewer-than-2 raises).
 
 - **Phase 15j/k/l — three more case studies, twelve installations,
   five continents.** Reinforces the urban-transit story (MGD and

README.md

@@ -104,6 +104,7 @@ ropeway run \
 | 13  | ✅ | AI assistant (LLM natural-language → optimization) — `ropeway ask "..."` |
 | 14  | ✅ | htmx web client (shareable links) — `/` + `/htmx/run/{rid}` |
 | 12d | ✅ | `ForcedFlyOverZone` — bridge / freeway / shipping-channel corridor intervals where the cable must clear a minimum absolute elevation (dual of `NoTowerZone`) |
+| 12f | ✅ | Curved corridor DEM extraction — `extract_profile_along_route(dem, waypoints, kind="polyline"|"bezier")` closes urban route-length gaps |
 | 15b-redux / 15c-redux | ✅ | Re-ran Línea Roja with Cementerio pinned (3-station topology now matches) + Whistler with Fitzsimmons no-tower valley (2 806 m single span, 7 % vs as-built 3 024 m, down from 47 %) |
 | 15d | ✅ | Case study: Roosevelt Island Tramway (NYC urban jig-back) |
 | 15e | ✅ | Case study: Medellín Metrocable Línea K (MGD, 4-station, 2 pinned waypoints) |

src/ropeway/dem.py

@@ -139,6 +139,70 @@ def extract_profile_from_dem(
     )
 
 
+def extract_profile_along_route(
+    dem_path: str | Path,
+    waypoints: "list[tuple[float, float]] | list",
+    *,
+    sample_spacing_m: float = 10.0,
+    kind: str = "polyline",
+) -> TerrainProfile:
+    """Phase 12f: sample a DEM along a *curved* corridor defined by waypoints.
+
+    ``waypoints`` is a list of (lon, lat) pairs **or** ``routes.Waypoint``
+    objects. ``kind`` is ``"polyline"`` (straight legs between waypoints)
+    or ``"bezier"`` (smooth Catmull-Rom through them). The returned
+    ``TerrainProfile.distance`` is the **along-route arc length** in
+    metres — exactly what the optimizer expects — so the result drops
+    straight into ``optimize(profile_fn=..., corridor_length=...)``
+    without any other change to downstream code. When the waypoints
+    collapse to just the two endpoints this returns the same profile as
+    ``extract_profile_from_dem`` (the straight-line case).
+    """
+    from .routes import RouteKind, Waypoint, build_route
+
+    if not _HAVE_RASTERIO:
+        raise RuntimeError("rasterio not installed; install with `pip install rasterio`")
+
+    pts: list[Waypoint] = []
+    for w in waypoints:
+        if isinstance(w, Waypoint):
+            pts.append(w)
+        else:
+            lon, lat = w
+            pts.append(Waypoint(lon=float(lon), lat=float(lat)))
+    if len(pts) < 2:
+        raise ValueError("need >= 2 waypoints")
+
+    route_kind = RouteKind(kind) if isinstance(kind, str) else kind
+    route = build_route(pts, kind=route_kind, sample_spacing_m=sample_spacing_m)
+
+    # Sample the DEM at each (xs_utm[i], ys_utm[i]).
+    with rasterio.open(str(dem_path)) as ds:
+        dem_crs = ds.crs
+        xs_dem, ys_dem = rio_transform(
+            CRS.from_epsg(route.utm_epsg), dem_crs,
+            route.xs_utm.tolist(), route.ys_utm.tolist(),
+        )
+        coords = list(zip(xs_dem, ys_dem))
+        elev = np.array([v[0] for v in ds.sample(coords)], dtype=float)
+        nodata = ds.nodata
+        if nodata is not None:
+            elev = np.where(elev == nodata, np.nan, elev)
+
+    if np.any(np.isnan(elev)):
+        mask = np.isnan(elev)
+        elev[mask] = np.interp(np.flatnonzero(mask),
+                               np.flatnonzero(~mask), elev[~mask])
+
+    return TerrainProfile(
+        distance=route.along.copy(),
+        elevation=elev,
+        start_lonlat=(pts[0].lon, pts[0].lat),
+        end_lonlat=(pts[-1].lon, pts[-1].lat),
+        utm_epsg=route.utm_epsg,
+    )
+
+
 def synthetic_profile(
     length_m: float = 3000.0,
     spacing_m: float = 10.0,

tests/test_curved_corridor.py

@@ -0,0 +1,107 @@
+"""Phase 12f — curved-corridor DEM extraction tests."""
+
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import pytest
+
+from ropeway.dem import (
+    extract_profile_along_route,
+    extract_profile_from_dem,
+)
+from ropeway.routes import Waypoint
+
+
+# Pick the smallest committed DEM tile so the test stays fast.
+_DEM = Path("data/dem/Copernicus_DSM_N22_E113.tif")
+pytestmark = pytest.mark.skipif(
+    not _DEM.exists(),
+    reason="DEM tile not present; case-study runs download it on demand",
+)
+
+
+def test_two_waypoint_polyline_matches_straight_extract():
+    """With only the two terminals, the curved-route extractor must
+    reproduce the straight-line ``extract_profile_from_dem`` result —
+    the polyline degenerates to a single straight leg."""
+    start = (113.94150, 22.28950)
+    end = (113.90500, 22.25650)
+    straight = extract_profile_from_dem(_DEM, start, end, sample_spacing_m=20.0)
+    curved = extract_profile_along_route(
+        _DEM, [start, end], sample_spacing_m=20.0, kind="polyline",
+    )
+    # Length and endpoint elevations match within DEM sampling noise.
+    assert curved.total_length == pytest.approx(straight.total_length, rel=1e-3)
+    assert curved.elevation[0] == pytest.approx(straight.elevation[0], abs=1.0)
+    assert curved.elevation[-1] == pytest.approx(straight.elevation[-1], abs=1.0)
+
+
+def test_added_waypoint_lengthens_corridor_vs_straight_line():
+    """An intermediate waypoint off the geodesic must produce a longer
+    total along-route distance than the two-endpoint straight line."""
+    start = (113.94150, 22.28950)
+    via = (113.93000, 22.27000)   # deliberate detour
+    end = (113.90500, 22.25650)
+    straight = extract_profile_from_dem(_DEM, start, end, sample_spacing_m=20.0)
+    curved = extract_profile_along_route(
+        _DEM, [start, via, end], sample_spacing_m=20.0, kind="polyline",
+    )
+    assert curved.total_length > straight.total_length
+
+
+def test_bezier_route_smooths_dogleg():
+    """Bezier through three waypoints must produce a longer corridor
+    than the polyline with the same waypoints (curves overshoot the
+    chord) — and both must exceed the straight start->end line."""
+    start = (113.94150, 22.28950)
+    via = (113.93200, 22.27200)
+    end = (113.90500, 22.25650)
+    poly = extract_profile_along_route(
+        _DEM, [start, via, end], sample_spacing_m=20.0, kind="polyline",
+    )
+    bez = extract_profile_along_route(
+        _DEM, [start, via, end], sample_spacing_m=20.0, kind="bezier",
+    )
+    straight = extract_profile_from_dem(_DEM, start, end, sample_spacing_m=20.0)
+    assert poly.total_length > straight.total_length
+    assert bez.total_length > straight.total_length
+    # Bezier and polyline can differ by a few % around the same chord;
+    # the important contract is "both lengthen the corridor".
+
+
+def test_returned_profile_drops_into_optimizer():
+    """The curved extractor returns a standard TerrainProfile — confirm
+    its ``as_function()`` is callable and finite, as the optimizer
+    consumes it that way."""
+    profile = extract_profile_along_route(
+        _DEM,
+        [(113.94150, 22.28950), (113.93000, 22.27000), (113.90500, 22.25650)],
+        sample_spacing_m=25.0, kind="polyline",
+    )
+    fn = profile.as_function()
+    xs = np.linspace(0.0, profile.total_length, 50)
+    ys = fn(xs)
+    assert np.all(np.isfinite(ys))
+    assert ys.shape == xs.shape
+
+
+def test_waypoint_objects_and_tuples_interchangeable():
+    """``waypoints`` accepts both raw (lon, lat) tuples and Waypoint
+    instances — the extractor normalises them."""
+    p_tuple = extract_profile_along_route(
+        _DEM, [(113.94150, 22.28950), (113.90500, 22.25650)],
+        sample_spacing_m=30.0,
+    )
+    p_obj = extract_profile_along_route(
+        _DEM,
+        [Waypoint(113.94150, 22.28950), Waypoint(113.90500, 22.25650)],
+        sample_spacing_m=30.0,
+    )
+    assert p_tuple.total_length == pytest.approx(p_obj.total_length, rel=1e-9)
+
+
+def test_too_few_waypoints_raises():
+    with pytest.raises(ValueError, match=">= 2 waypoints"):
+        extract_profile_along_route(_DEM, [(113.94150, 22.28950)])