diff --git a/CHANGELOG.md b/CHANGELOG.md index 886fc79..0a3a435 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/README.md b/README.md index 466f795..d31c839 100644 --- a/README.md +++ b/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) | diff --git a/src/ropeway/dem.py b/src/ropeway/dem.py index 2ff1ee2..0d05e31 100644 --- a/src/ropeway/dem.py +++ b/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, diff --git a/tests/test_curved_corridor.py b/tests/test_curved_corridor.py new file mode 100644 index 0000000..1eb4247 --- /dev/null +++ b/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)])