first commit

api/app/core/surface.py

@@ -0,0 +1,75 @@
+from __future__ import annotations
+
+from dataclasses import dataclass
+from typing import Literal
+
+from app.core.geo import PathPoint
+
+ObstructionKind = Literal["terrain", "building", "canopy"]
+
+
+@dataclass(frozen=True)
+class SurfaceSample:
+    i: int
+    lat: float
+    lon: float
+    distance_m: float
+    ground_m: float
+    building_m: float
+    canopy_m: float
+    surface_m: float
+    curvature_corr_m: float = 0.0
+
+    @property
+    def dominant_obstruction(self) -> ObstructionKind:
+        if self.building_m > 0 and self.ground_m + self.building_m >= self.surface_m:
+            return "building"
+        if self.canopy_m > 0 and self.ground_m + self.canopy_m >= self.surface_m:
+            return "canopy"
+        return "terrain"
+
+
+@dataclass(frozen=True)
+class SurfaceProfile:
+    distance_m: float
+    samples: list[SurfaceSample]
+
+
+def build_surface_profile(
+    points: list[PathPoint],
+    ground_elevations: list[float] | None = None,
+    building_heights: list[float] | None = None,
+    canopy_heights: list[float] | None = None,
+    include_buildings: bool = True,
+    include_canopy: bool = True,
+) -> SurfaceProfile:
+    if not points:
+        raise ValueError("points must not be empty")
+
+    count = len(points)
+    ground_elevations = ground_elevations or [0.0] * count
+    building_heights = building_heights or [0.0] * count
+    canopy_heights = canopy_heights or [0.0] * count
+    if not (len(ground_elevations) == len(building_heights) == len(canopy_heights) == count):
+        raise ValueError("profile arrays must match points length")
+
+    samples: list[SurfaceSample] = []
+    for i, point in enumerate(points):
+        ground_m = float(ground_elevations[i])
+        building_m = float(building_heights[i]) if include_buildings else 0.0
+        canopy_m = float(canopy_heights[i]) if include_canopy else 0.0
+        surface_m = ground_m + max(building_m, canopy_m)
+        samples.append(
+            SurfaceSample(
+                i=i,
+                lat=point.lat,
+                lon=point.lon,
+                distance_m=point.distance_m,
+                ground_m=ground_m,
+                building_m=building_m,
+                canopy_m=canopy_m,
+                surface_m=surface_m,
+            )
+        )
+
+    return SurfaceProfile(distance_m=points[-1].distance_m, samples=samples)