"""Elevation profile with a smoothed grade for a saved route (docs/profile_storms_v0_1_design.md section 2.1). The server does all of it -- resampling, smoothing, the grade itself -- and the browser only maps the finished ``grade_pct`` values to colours. The method follows the spike (docs/design_reference/profile_storms_spike_v0_1.md section b): the route's own elevations are resampled every 10 m, smoothed by a moving average over a 600 m window, and the grade is the central difference of the smoothed curve across the same window, in whole percent (BRouter elevations agree with the Copernicus DEM to about 2 pp at this window; a finer figure would claim precision the data does not have). Nothing is fabricated: a point without an elevation stays without one, a grade whose window reaches a missing elevation is ``None``, and trail marking (colours, names, POI) is explicitly UNAVAILABLE -- Phase 1 has no offline OSM relation or node data for any region (spike section a). Grade is a measurement of terrain slope, never a difficulty or safety rating. """ from __future__ import annotations import math from collections.abc import Sequence from typing import Any from mountain_twin.route_analysis import prepare_route ROUTE_PROFILE_CONTRACT = "route_profile_v0_1" GRADE_WINDOW_M = 600.0 RESAMPLE_STEP_M = 10.0 OUTPUT_STEP_M = 50.0 MAX_OUTPUT_SAMPLES = 1500 # Near the start and the end the window is cut to the route; a grade needs # its two smoothed ends at least a quarter window apart (less is not a grade). MIN_GRADE_SPAN_M = GRADE_WINDOW_M / 4 TRAIL_MARKING_UNAVAILABLE = { "state": "UNAVAILABLE", "reason": "NO_TRAIL_MARKING_DATA", "explanation": "brak danych o oznakowaniu", } GRADE_METHOD = { "window_m": GRADE_WINDOW_M, "resample_step_m": RESAMPLE_STEP_M, "smoothing": "MOVING_AVERAGE", "grade": "CENTRAL_DIFFERENCE_OF_SMOOTHED_ELEVATION", "precision": "WHOLE_PERCENT", "elevation_source": "ROUTE_GEOMETRY", } GRADE_LIMITATIONS = ( "Moduł nachylenia, wygładzony (okno ok. 600 m); dokładność ogranicza model wysokości.", "To pomiar spadku terenu, nie ocena trudności ani bezpieczeństwa.", ) def _interpolate(distances, elevations, at): """Elevation at ``at`` metres by linear interpolation between the two route points around it -- ``None`` if either has no elevation.""" low, high = 0, len(distances) - 1 while high - low > 1: middle = (low + high) // 2 if distances[middle] <= at: low = middle else: high = middle a, b = elevations[low], elevations[high] if a is None or b is None: return None span = distances[high] - distances[low] if span <= 0: return a return a + (b - a) * (at - distances[low]) / span def _moving_average(values, half): """Mean over [i - half, i + half] (cut to the series), ``None`` if any value in the window is missing.""" result = [] prefix, missing = [0.0], [0] for value in values: prefix.append(prefix[-1] + (value or 0.0)) missing.append(missing[-1] + (value is None)) for index in range(len(values)): low, high = max(0, index - half), min(len(values) - 1, index + half) if missing[high + 1] - missing[low]: result.append(None) else: result.append((prefix[high + 1] - prefix[low]) / (high - low + 1)) return result def _window_centres(grid, half): """The mean distance of each (cut) smoothing window: where its average elevation actually stands -- off the grid point near the route's ends.""" return _moving_average(grid, half) def grade_profile( distances_m: Sequence[float], elevations_m: Sequence[float | None] ) -> list[dict[str, Any]]: """Samples ``{distance_m, elevation_m, grade_pct}`` every OUTPUT_STEP_M (coarser for a long route, at most MAX_OUTPUT_SAMPLES) along one route whose cumulative distances and elevations are given point by point.""" if len(distances_m) != len(elevations_m) or len(distances_m) < 2: raise ValueError("a profile needs at least two route points") length = distances_m[-1] if length <= 0: raise ValueError("a profile needs a route with a length") grid_count = int(length // RESAMPLE_STEP_M) + 1 grid = [min(length, index * RESAMPLE_STEP_M) for index in range(grid_count)] if grid[-1] < length: grid.append(length) raw = [_interpolate(distances_m, elevations_m, at) for at in grid] half = int(round(GRADE_WINDOW_M / 2 / RESAMPLE_STEP_M)) smooth = _moving_average(raw, half) centres = _window_centres(grid, half) def grade_at(index): # Central difference across the window, divided by the distance # between the two smoothed values' own centres (exact for an even # slope, also where the window is cut at the route's ends). low, high = max(0, index - half), min(len(grid) - 1, index + half) span = centres[high] - centres[low] if span < MIN_GRADE_SPAN_M - 1 or smooth[low] is None or smooth[high] is None: return None return int(round((smooth[high] - smooth[low]) / span * 100)) step = max(OUTPUT_STEP_M, length / (MAX_OUTPUT_SAMPLES - 1)) stride = max(1, math.ceil(step / RESAMPLE_STEP_M - 1e-9)) indexes = list(range(0, len(grid), stride)) if indexes[-1] != len(grid) - 1: indexes.append(len(grid) - 1) return [ { "distance_m": round(grid[index], 1), "elevation_m": None if raw[index] is None else round(raw[index], 1), "grade_pct": grade_at(index), "trail": None, } for index in indexes ] def _movement_grade_profile(prepared) -> list[dict[str, Any]]: """grade_profile() per movement stretch (AV-064): a transfer (dojazd) is its own segment in the movement view, so each stretch is profiled on its own and a sample without height marks the break -- no grade is ever computed across a train ride.""" blocks: list[list[Any]] = [] for point in prepared: if blocks and blocks[-1][-1].segment_index == point.segment_index: blocks[-1].append(point) else: blocks.append([point]) blocks = [ block for block in blocks if len(block) >= 2 and block[-1].cumulative_distance_m > block[0].cumulative_distance_m ] if len(blocks) <= 1: points = blocks[0] if blocks else prepared return grade_profile( [point.cumulative_distance_m for point in points], [point.gpx_elevation_m for point in points], ) samples: list[dict[str, Any]] = [] for block in blocks: start = block[0].cumulative_distance_m part = grade_profile( [point.cumulative_distance_m - start for point in block], [point.gpx_elevation_m for point in block], ) if samples: samples.append( { "distance_m": round(start, 1), "elevation_m": None, "grade_pct": None, "trail": None, "transfer": True, } ) samples.extend( {**sample, "distance_m": round(sample["distance_m"] + start, 1)} for sample in part ) return samples def route_profile_document( route_points: Sequence[Any], *, days: Sequence[Any] = (), camps: Sequence[Any] = (), elevation_source: str = "ROUTE_GEOMETRY", ) -> dict[str, Any]: """The whole route's profile with its labelled points. ``elevation_coverage`` says where the elevations come from (``elevation_source``: ``BROUTER`` for a route drawn in Planning Workspace, ``ROUTE_GEOMETRY`` for a preserved reference geometry) and how many route points have none -- those stay missing (UNAVAILABLE), never filled in (docs/design_reference/elevation_gaps_spike_v0_1.md). ``days`` are DaySegment-like (``day_number``, ``start_point_index``, ``end_point_index``); ``camps`` are CampMarker-like (``route_point_index``, ``label``). Labels are only what the route and the plan know: start, finish, each camp under the user's own label, and the highest point of the route and of each day. No place names (no OSM data in Phase 1).""" prepared = prepare_route(route_points).points distances = [point.cumulative_distance_m for point in prepared] elevations = [point.gpx_elevation_m for point in prepared] distance_by_index = {point.point_index: point.cumulative_distance_m for point in prepared} elevation_by_index = {point.point_index: point.gpx_elevation_m for point in prepared} def highest(first_index, last_index): candidates = [ point for point in prepared if first_index <= point.point_index <= last_index and point.gpx_elevation_m is not None ] if not candidates: return None top = max(candidates, key=lambda point: point.gpx_elevation_m) return { "distance_m": round(top.cumulative_distance_m, 1), "elevation_m": round(top.gpx_elevation_m, 1), } def label(kind, point_index, name=None): elevation = elevation_by_index.get(point_index) return { "kind": kind, "name": name, "distance_m": round(distance_by_index[point_index], 1), "elevation_m": None if elevation is None else round(elevation, 1), } first_index, last_index = prepared[0].point_index, prepared[-1].point_index labels = [label("START", first_index), label("FINISH", last_index)] labels += [ label("CAMP", camp.route_point_index, camp.label) for camp in camps if camp.route_point_index in distance_by_index ] return { "contract": ROUTE_PROFILE_CONTRACT, "length_m": round(distances[-1], 1), "samples": _movement_grade_profile(prepared), "labels": labels, "highest_point": highest(first_index, last_index), "days": [ { "day_number": day.day_number, "start_distance_m": round(distance_by_index[day.start_point_index], 1), "end_distance_m": round(distance_by_index[day.end_point_index], 1), "highest_point": highest(day.start_point_index, day.end_point_index), } for day in days ], "elevation_coverage": { "source": elevation_source, "points": len(prepared), "unavailable_points": sum(1 for elevation in elevations if elevation is None), }, "grade_method": dict(GRADE_METHOD), "grade_limitations": list(GRADE_LIMITATIONS), "trail_marking": dict(TRAIL_MARKING_UNAVAILABLE), "points_of_interest": { "state": "UNAVAILABLE", "reason": "NO_OFFLINE_POI_DATA", "explanation": "brak danych o punktach (woda, jedzenie, pomoc)", }, } __all__ = [ "GRADE_WINDOW_M", "MAX_OUTPUT_SAMPLES", "ROUTE_PROFILE_CONTRACT", "grade_profile", "route_profile_document", ]