"""Provider-labelled, explicit-time route surface-illumination analysis.""" from __future__ import annotations import time from dataclasses import asdict, dataclass from datetime import datetime from enum import Enum from typing import Callable, Iterable, Sequence from mountain_twin.exposure import RoutePoint from mountain_twin.route_analysis import prepare_route from mountain_twin.solar.states import ( AstronomicalState, TerrainSolarVisibility, classify_astronomical_state, ) from mountain_twin.solar.sun_engine import solar_position from mountain_twin.terrain.surface_reference import continuous_horizon, geodesic_path class VisibilityStatus(str, Enum): """A route-surface solar result; UNKNOWN is a deliberate product state.""" DIRECT = "DIRECT" SHADOW = "SHADOW" NOT_APPLICABLE = "NOT_APPLICABLE" UNKNOWN = "UNKNOWN" class SemanticType(str, Enum): """Names deliberately separate target and observer interpretations.""" SURFACE_ILLUMINATION = "surface_illumination" HUMAN_SUN_EXPOSURE = "human_sun_exposure" CAMERA_SUN_VISIBILITY = "camera_sun_visibility" TARGET_ILLUMINATION = "target_illumination" class ReasonCode(str, Enum): SUN_BELOW_ASTRONOMICAL_HORIZON = "SUN_BELOW_ASTRONOMICAL_HORIZON" TERRAIN_VISIBILITY_NOT_APPLICABLE = "TERRAIN_VISIBILITY_NOT_APPLICABLE" HORIZON_SUPPORT_INCOMPLETE = "HORIZON_SUPPORT_INCOMPLETE" HORIZON_NUMERICAL_CONVERGENCE_UNMET = "HORIZON_NUMERICAL_CONVERGENCE_UNMET" HORIZON_RANGE_NOT_CONVERGED = "HORIZON_RANGE_NOT_CONVERGED" GLOBAL_LIMITATIONS = ( "GLO30_DSM_NOT_BARE_EARTH_DTM", "GLO30_FINE_NEAR_FIELD_HUMAN_OCCLUDERS_NOT_RESOLVED", "HORIZON_RANGE_5KM_NOT_GLOBALLY_CONVERGED", "NO_CLOUDS_WEATHER_VEGETATION_TRANSMISSIVITY_OR_BUILDING_MODEL", "GEOMETRIC_SURFACE_ILLUMINATION_NOT_HUMAN_EXPOSURE_OR_SAFETY_ASSESSMENT", ) @dataclass(frozen=True) class TerrainContext: """Provider metadata shared by all compact point results in one analysis.""" context_id: str provider: str product: str | None surface_semantics: str native_resolution: tuple[float, float] | None horizontal_crs: str horizontal_units: tuple[str, str] vertical_reference: str | None vertical_reference_verified: bool @dataclass(frozen=True) class RoutePointSolarExposureResult: """One point at one explicitly supplied scenario instant. ``terrain_context_id`` resolves to the dataset-level TerrainContext rather than repeating provenance and limitations 7,575 times. """ route_id: str point_index: int latitude: float longitude: float route_distance_m: float route_elevation_m: float | None scenario_datetime: str timezone: str solar_azimuth_deg: float solar_elevation_deg: float status: VisibilityStatus horizon_angle_deg: float | None controlling_distance_m: float | None observer_height_m: float observer_model: str semantic_type: SemanticType terrain_context_id: str analysis_method: str horizon_range_m: float reason_codes: tuple[str, ...] ranges_tested_m: tuple[float, ...] = () converged_at_m: float | None = None range_convergence_state: str = "not_evaluated" astronomical_state: AstronomicalState | None = None terrain_solar_visibility: TerrainSolarVisibility | None = None @dataclass(frozen=True) class RuntimeMetrics: total_seconds: float horizon_seconds: float horizon_calculations: int @property def seconds_per_horizon(self) -> float: return self.horizon_seconds / self.horizon_calculations def scenario_instants( date: str = "2026-07-02", timezone_name: str = "Europe/Paris" ) -> tuple[datetime, ...]: """Return only supplied scenario instants; route/GPX timestamps never enter.""" from mountain_twin.exposure import resolve_local_time return tuple( resolve_local_time(f"{date}T{hour:02}:00", timezone_name) for hour in (6, 9, 12, 15, 18) ) def expected_result_count(points: Sequence[RoutePoint], instants: Sequence[datetime]) -> int: return len(points) * len(instants) def terrain_context_from_surface(surface) -> TerrainContext: info = surface.dem.info metadata = surface.dem.metadata return TerrainContext( context_id="copernicus_glo30_continuous_bilinear_surface_v1", provider=info.source, product=info.product, surface_semantics=metadata.get("source_type", "DSM; semantics require verification"), native_resolution=info.horizontal_resolution, horizontal_crs=info.crs, horizontal_units=info.horizontal_units, vertical_reference=info.vertical_reference, vertical_reference_verified=info.vertical_reference_verified, ) def classify_visibility( solar_elevation_deg: float, horizon: dict, horizon_range_m: float, range_boundary_margin_m: float, ) -> tuple[VisibilityStatus, tuple[str, ...]]: """Classify finite-range terrain geometry without treating UNKNOWN as null. A known blocking horizon remains SHADOW even near the finite range boundary. A direct result at that boundary is UNKNOWN because unmodelled farther terrain could still block the sun. This is a finite-model policy, not a claim that every non-edge direct horizon has converged beyond 5 km. """ if solar_elevation_deg < 0: return VisibilityStatus.NOT_APPLICABLE, ( ReasonCode.TERRAIN_VISIBILITY_NOT_APPLICABLE.value, ) if horizon["horizon_status"] == "numerical_convergence_unmet": return VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_NUMERICAL_CONVERGENCE_UNMET.value,) if horizon["horizon_status"] != "complete" or horizon["horizon_angle_deg"] is None: return VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_SUPPORT_INCOMPLETE.value,) distance = horizon["maximum"]["distance_m"] range_edge = distance >= horizon_range_m - range_boundary_margin_m blocked = horizon["horizon_angle_deg"] >= solar_elevation_deg if blocked: reasons = (ReasonCode.HORIZON_RANGE_NOT_CONVERGED.value,) if range_edge else () return VisibilityStatus.SHADOW, reasons if range_edge: return VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_RANGE_NOT_CONVERGED.value,) return VisibilityStatus.DIRECT, () def analyze_surface_illumination( points: Sequence[RoutePoint], surface, instants: Sequence[datetime], *, horizon_range_m: float = 5000, range_boundary_margin_m: float = 100, horizon_angle_tolerance_deg: float = 1e-5, solar_fn: Callable[[datetime, float, float], tuple[float, float]] = solar_position, ) -> tuple[tuple[RoutePointSolarExposureResult, ...], RuntimeMetrics]: """Evaluate route-surface illumination with the continuous bilinear model.""" if horizon_range_m <= 0 or range_boundary_margin_m < 0: raise ValueError("invalid horizon range policy") if any(instant.utcoffset() is None for instant in instants): raise ValueError("scenario instants must be timezone-aware") prepared = prepare_route(points).points if len(prepared) != len(points): raise ValueError("prepared route point count mismatch") context = terrain_context_from_surface(surface) started, horizon_seconds, calculations = time.perf_counter(), 0.0, 0 results = [] for instant in instants: timezone_name = getattr(instant.tzinfo, "key", str(instant.tzinfo)) for point, route_point in zip(points, prepared): solar_elevation, solar_azimuth = solar_fn(instant, point.latitude, point.longitude) horizon = { "horizon_status": "not_evaluated", "horizon_angle_deg": None, "maximum": None, } if solar_elevation > 0: ray = geodesic_path(surface, point.latitude, point.longitude, solar_azimuth) began = time.perf_counter() try: horizon = continuous_horizon( surface, ray, horizon_range_m, horizon_angle_tolerance_deg, observer_height_m=0, ) except ValueError as error: if str(error) != "continuous reference did not meet angle tolerance": raise horizon = { "horizon_status": "numerical_convergence_unmet", "horizon_angle_deg": None, "maximum": None, } horizon_seconds += time.perf_counter() - began calculations += 1 status, reasons = classify_visibility( solar_elevation, horizon, horizon_range_m, range_boundary_margin_m ) maximum = horizon["maximum"] results.append( RoutePointSolarExposureResult( route_id=point.route_id, point_index=point.point_index, latitude=point.latitude, longitude=point.longitude, route_distance_m=route_point.cumulative_distance_m, route_elevation_m=point.elevation_m, scenario_datetime=instant.isoformat(), timezone=timezone_name, solar_azimuth_deg=solar_azimuth, solar_elevation_deg=solar_elevation, status=status, horizon_angle_deg=horizon["horizon_angle_deg"], controlling_distance_m=maximum["distance_m"] if maximum else None, observer_height_m=0.0, observer_model="provider_surface_bilinear_zero_height", semantic_type=SemanticType.SURFACE_ILLUMINATION, terrain_context_id=context.context_id, analysis_method="continuous_bilinear_grid_post_horizon_v1", horizon_range_m=horizon_range_m, reason_codes=reasons, astronomical_state=classify_astronomical_state(solar_elevation), terrain_solar_visibility=TerrainSolarVisibility(status.value), ) ) return tuple(results), RuntimeMetrics( total_seconds=time.perf_counter() - started, horizon_seconds=horizon_seconds, horizon_calculations=calculations, ) def _geometric_status(horizon: dict, solar_elevation_deg: float) -> VisibilityStatus | None: """Return direct/shadow geometry even when finite-range policy made it UNKNOWN.""" if horizon.get("horizon_status") != "complete" or horizon.get("horizon_angle_deg") is None: return None return ( VisibilityStatus.SHADOW if horizon["horizon_angle_deg"] >= solar_elevation_deg else VisibilityStatus.DIRECT ) def analyze_surface_illumination_adaptive( points: Sequence[RoutePoint], surfaces: dict[int, object], instants: Sequence[datetime], *, ranges_m: Sequence[int] = (5000, 7500, 10000), range_boundary_margin_m: float = 100, horizon_angle_tolerance_deg: float = 1e-5, angle_convergence_tolerance_deg: float = 0.25, solar_fn: Callable[[datetime, float, float], tuple[float, float]] = solar_position, ) -> tuple[tuple[RoutePointSolarExposureResult, ...], RuntimeMetrics]: """Evaluate finite terrain range adaptively, stopping at measured convergence. A range is converged when the controller is at least ``range_boundary_margin_m`` inside the tested extent and, after the first range, the geometric class and horizon angle are stable within ``angle_convergence_tolerance_deg``. A final unresolved range is UNKNOWN, including a known blocking feature at the boundary. """ ranges = tuple(int(value) for value in ranges_m) if ranges != tuple(sorted(ranges)) or not ranges or any(value <= 0 for value in ranges): raise ValueError("ranges must be positive and ascending") if range_boundary_margin_m < 0 or angle_convergence_tolerance_deg <= 0: raise ValueError("invalid adaptive range policy") if set(ranges) - set(surfaces): raise ValueError("a surface is required for every adaptive range") if any(instant.utcoffset() is None for instant in instants): raise ValueError("scenario instants must be timezone-aware") prepared = prepare_route(points).points if len(prepared) != len(points): raise ValueError("prepared route point count mismatch") baseline_context = terrain_context_from_surface(surfaces[ranges[0]]) started, horizon_seconds, calculations = time.perf_counter(), 0.0, 0 results = [] for instant in instants: timezone_name = getattr(instant.tzinfo, "key", str(instant.tzinfo)) for point, route_point in zip(points, prepared): solar_elevation, solar_azimuth = solar_fn(instant, point.latitude, point.longitude) tested, previous = [], None final = None converged_at = None state = "not_applicable" if solar_elevation < 0 else "not_converged" if solar_elevation > 0: for range_m in ranges: surface = surfaces[range_m] ray = geodesic_path(surface, point.latitude, point.longitude, solar_azimuth) began = time.perf_counter() try: horizon = continuous_horizon( surface, ray, range_m, horizon_angle_tolerance_deg, observer_height_m=0, ) except ValueError as error: if str(error) != "continuous reference did not meet angle tolerance": raise horizon = { "horizon_status": "numerical_convergence_unmet", "horizon_angle_deg": None, "maximum": None, } horizon_seconds += time.perf_counter() - began calculations += 1 tested.append(range_m) final = (range_m, horizon) if horizon["horizon_status"] == "numerical_convergence_unmet": state = "numerical_unknown" break if horizon["horizon_status"] != "complete": state = "support_unknown" previous = horizon continue geometric = _geometric_status(horizon, solar_elevation) distance = horizon["maximum"]["distance_m"] inside = distance < range_m - range_boundary_margin_m angle_stable = previous is None or ( _geometric_status(previous, solar_elevation) == geometric and abs(horizon["horizon_angle_deg"] - previous["horizon_angle_deg"]) <= angle_convergence_tolerance_deg ) if inside and angle_stable: converged_at, state = range_m, "converged" break previous = horizon range_m, horizon = final or ( ranges[0], {"horizon_status": "not_evaluated", "horizon_angle_deg": None, "maximum": None}, ) if solar_elevation < 0: status, reasons = ( VisibilityStatus.NOT_APPLICABLE, (ReasonCode.TERRAIN_VISIBILITY_NOT_APPLICABLE.value,), ) elif state == "numerical_unknown": status, reasons = ( VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_NUMERICAL_CONVERGENCE_UNMET.value,), ) elif state == "support_unknown": status, reasons = ( VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_SUPPORT_INCOMPLETE.value,), ) elif converged_at is None: status, reasons = ( VisibilityStatus.UNKNOWN, (ReasonCode.HORIZON_RANGE_NOT_CONVERGED.value,), ) else: status, reasons = classify_visibility( solar_elevation, horizon, range_m, range_boundary_margin_m ) maximum = horizon.get("maximum") or {} results.append( RoutePointSolarExposureResult( route_id=point.route_id, point_index=point.point_index, latitude=point.latitude, longitude=point.longitude, route_distance_m=route_point.cumulative_distance_m, route_elevation_m=point.elevation_m, scenario_datetime=instant.isoformat(), timezone=timezone_name, solar_azimuth_deg=solar_azimuth, solar_elevation_deg=solar_elevation, status=status, horizon_angle_deg=horizon.get("horizon_angle_deg"), controlling_distance_m=maximum.get("distance_m"), observer_height_m=0.0, observer_model="provider_surface_bilinear_zero_height", semantic_type=SemanticType.SURFACE_ILLUMINATION, terrain_context_id=baseline_context.context_id, analysis_method="adaptive_continuous_bilinear_grid_post_horizon_v1", horizon_range_m=range_m, reason_codes=reasons, ranges_tested_m=tuple(tested), converged_at_m=converged_at, range_convergence_state=state, astronomical_state=classify_astronomical_state(solar_elevation), terrain_solar_visibility=TerrainSolarVisibility(status.value), ) ) return tuple(results), RuntimeMetrics( total_seconds=time.perf_counter() - started, horizon_seconds=horizon_seconds, horizon_calculations=calculations, ) def _sections(records: Iterable[RoutePointSolarExposureResult]) -> list[dict]: rows = list(records) if not rows: return [] sections, first = [], rows[0] previous = first for row in rows[1:]: if row.status != first.status: sections.append(_section(first, previous)) first = row previous = row sections.append(_section(first, previous)) return sections def _section(first: RoutePointSolarExposureResult, last: RoutePointSolarExposureResult) -> dict: return { "status": first.status.value, "start_point_index": first.point_index, "end_point_index": last.point_index, "start_distance_m": first.route_distance_m, "end_distance_m": last.route_distance_m, "length_m": last.route_distance_m - first.route_distance_m, } def scenario_summaries( results: Sequence[RoutePointSolarExposureResult], *, major_section_min_m: float = 500 ) -> list[dict]: """Spatial-snapshot summaries; they intentionally do not estimate time in sun.""" by_scenario: dict[str, list[RoutePointSolarExposureResult]] = {} for result in results: by_scenario.setdefault(result.scenario_datetime, []).append(result) summaries = [] for scenario, rows in by_scenario.items(): rows.sort(key=lambda row: row.point_index) total = len(rows) counts = { status.value: sum(row.status is status for row in rows) for status in VisibilityStatus } sections = _sections(rows) direct_shadow = [ section for section in sections if section["status"] != VisibilityStatus.UNKNOWN.value ] major = [section for section in direct_shadow if section["length_m"] >= major_section_min_m] summaries.append( { "scenario_datetime": scenario, "timezone": rows[0].timezone, "counts": counts, "percentages": {name: count * 100 / total for name, count in counts.items()}, "longest_contiguous_direct_section_m": max( (section["length_m"] for section in sections if section["status"] == "DIRECT"), default=0.0, ), "longest_contiguous_shadow_section_m": max( (section["length_m"] for section in sections if section["status"] == "SHADOW"), default=0.0, ), "longest_contiguous_unknown_section_m": max( (section["length_m"] for section in sections if section["status"] == "UNKNOWN"), default=0.0, ), "visibility_transitions": len(sections) - 1, "major_direct_shadow_sections": major, } ) return summaries def analysis_document( points: Sequence[RoutePoint], results: Sequence[RoutePointSolarExposureResult], surface, instants: Sequence[datetime], *, horizon_range_m: float, range_boundary_margin_m: float, ) -> dict: """Build deterministic structured output; measured timing is deliberately separate.""" prepared = prepare_route(points).points context = terrain_context_from_surface(surface) route_elevations = [point.elevation_m for point in points if point.elevation_m is not None] adaptive = any(result.ranges_tested_m for result in results) range_counts = {str(value): 0 for value in (5000, 7500, 10000)} evaluation_counts = {str(value): 0 for value in (5000, 7500, 10000)} for result in results: if result.ranges_tested_m: range_counts[str(result.ranges_tested_m[-1])] += 1 for value in result.ranges_tested_m: evaluation_counts[str(value)] += 1 return { "contract_version": "route_solar_exposure_v0_1", "analysis_id": ( "tmb_day_01_glo30_surface_illumination_adaptive_20260702_v0_1" if adaptive else "tmb_day_01_glo30_surface_illumination_20260702_v0_1" ), "time_semantics": "explicit_simultaneous_scenarios_not_gpx_or_hiking_times", "route": { "route_id": "tmb_day_01", "point_count": len(points), "distance_m": prepared[-1].cumulative_distance_m, "route_elevation_source": "GPX-derived; vertical reference unknown/unverified", "minimum_route_elevation_m": min(route_elevations), "maximum_route_elevation_m": max(route_elevations), }, "scenario_datetimes": [instant.isoformat() for instant in instants], "semantic_type": SemanticType.SURFACE_ILLUMINATION.value, "observer_height_m": 0.0, "observer_model": "provider_surface_bilinear_zero_height", "terrain": asdict(context), "analysis_method": ( "adaptive_continuous_bilinear_grid_post_horizon_v1" if adaptive else "continuous_bilinear_grid_post_horizon_v1" ), "horizon_range_m": max( (result.horizon_range_m for result in results), default=horizon_range_m ), "range_boundary_margin_m": range_boundary_margin_m, "range_policy": { "ranges_m": [5000, 7500, 10000] if adaptive else [horizon_range_m], "stop_when": "complete support, stable geometric classification and <=0.25 degree angle change; first range requires controller >100 m inside boundary", "unresolved_at_final_range": "UNKNOWN with HORIZON_RANGE_NOT_CONVERGED", "direct_at_range_edge": "UNKNOWN with HORIZON_RANGE_NOT_CONVERGED", "shadow_at_range_edge": "SHADOW with HORIZON_RANGE_NOT_CONVERGED; known terrain already blocks sun", }, "range_statistics": range_counts, "final_range_counts": range_counts, "solver_evaluation_counts": evaluation_counts, "limitations": list(GLOBAL_LIMITATIONS), "expected_result_count": expected_result_count(points, instants), "results": [asdict(result) for result in results], "scenario_summaries": scenario_summaries(results), }