first commit

app/services/geodesy.py

@@ -0,0 +1,73 @@
+from __future__ import annotations
+
+import math
+
+from pyproj import Geod
+
+_WGS84 = Geod(ellps="WGS84")
+
+
+def sample_line(
+  from_lat: float,
+  from_lng: float,
+  to_lat: float,
+  to_lng: float,
+  step_m: float,
+) -> list[tuple[float, float, float]]:
+  """Return (lat, lng, distance_from_start_m) points along a geodesic line."""
+  _, _, total_distance = _WGS84.inv(from_lng, from_lat, to_lng, to_lat)
+  total_distance = abs(total_distance)
+
+  if total_distance == 0:
+    return [(from_lat, from_lng, 0.0)]
+
+  num_steps = max(1, int(math.ceil(total_distance / step_m)))
+  actual_step = total_distance / num_steps
+
+  azimuth, _, _ = _WGS84.inv(from_lng, from_lat, to_lng, to_lat)
+  points: list[tuple[float, float, float]] = []
+  for step_index in range(num_steps + 1):
+    distance = min(step_index * actual_step, total_distance)
+    if step_index == num_steps:
+      lat, lng = to_lat, to_lng
+    else:
+      lng, lat, _ = _WGS84.fwd(from_lng, from_lat, azimuth, distance)
+    points.append((lat, lng, distance))
+
+  return points
+
+
+def line_distance_m(from_lat: float, from_lng: float, to_lat: float, to_lng: float) -> float:
+  _, _, distance = _WGS84.inv(from_lng, from_lat, to_lng, to_lat)
+  return abs(distance)
+
+
+def los_height_at(distance_m: float, from_antenna_asl: float, to_antenna_asl: float, total_distance_m: float) -> float:
+  if total_distance_m <= 0:
+    return from_antenna_asl
+  fraction = distance_m / total_distance_m
+  return from_antenna_asl + (to_antenna_asl - from_antenna_asl) * fraction
+
+
+def grid_in_circle(
+  origin_lat: float,
+  origin_lng: float,
+  radius_m: float,
+  resolution_m: float,
+) -> list[tuple[float, float]]:
+  """Generate grid points inside a circle around origin."""
+  points: list[tuple[float, float]] = []
+  steps = int(math.ceil(radius_m / resolution_m))
+
+  for row in range(-steps, steps + 1):
+    for col in range(-steps, steps + 1):
+      easting = col * resolution_m
+      northing = row * resolution_m
+      if math.hypot(easting, northing) > radius_m:
+        continue
+
+      lng, lat, _ = _WGS84.fwd(origin_lng, origin_lat, 90.0, easting)
+      lng, lat, _ = _WGS84.fwd(lng, lat, 0.0, northing)
+      points.append((lat, lng))
+
+  return points