added formulas

api/app/core/coverage.py

@@ -1,17 +1,28 @@
 from __future__ import annotations
 
+from pathlib import Path
 from typing import Any
+from uuid import uuid4
 
+from affine import Affine
+import numpy as np
+import rasterio
 from pyproj import Geod
+from pyproj import Transformer
+from rasterio.transform import from_origin
+from sqlalchemy.orm import Session
 
 from app.core.antenna import AntennaPattern, gain
+from app.core.diffraction import bullington_equivalent_loss
 from app.core.geo import GeoPoint, sample_path
 from app.core.propagation import fspl, itm_loss, p1812_field
 from app.models.antenna import AntennaPatternSpec
 from app.models.coverage import CoverageRequest
+from app.services.terrain import surface_profile_from_points
 from app.services.vegetation import vegetation_loss_along
 
 _GEOD = Geod(ellps="WGS84")
+_NODATA_DBM = -9999.0
 
 
 def _antenna_pattern(spec: AntennaPatternSpec) -> AntennaPattern:
@@ -66,10 +77,13 @@ def _rx_power_dbm(
     distance_km: float,
     azimuth_deg: float,
     elevation_profile: list[float] | None = None,
+    db: Session | None = None,
+    rx: GeoPoint | None = None,
 ) -> float:
     tx = GeoPoint(lat=request.tx.lat, lon=request.tx.lon)
-    rx_point = _GEOD.fwd(request.tx.lon, request.tx.lat, azimuth_deg, distance_km * 1000)
-    rx = GeoPoint(lat=rx_point[1], lon=rx_point[0])
+    if rx is None:
+        rx_point = _GEOD.fwd(request.tx.lon, request.tx.lat, azimuth_deg, distance_km * 1000)
+        rx = GeoPoint(lat=rx_point[1], lon=rx_point[0])
     pattern = _antenna_pattern(request.antenna)
     antenna_gain = gain(pattern, azimuth_deg, 0.0)
     path_loss = _path_loss_db(request, tx, rx, distance_km, azimuth_deg, elevation_profile)
@@ -81,14 +95,251 @@ def _rx_power_dbm(
             freq_hz=request.tx.frequency_mhz * 1_000_000,
             include_vegetation=True,
         )
+    surface_obstruction_db = _surface_obstruction_loss_db(request, tx, rx, db)
     eirp = request.tx.power_dbm + antenna_gain
-    return eirp - path_loss - vegetation_db + request.rx.gain_dbi
+    return eirp - path_loss - vegetation_db - surface_obstruction_db + request.rx.gain_dbi
+
+
+def _surface_obstruction_loss_db(
+    request: CoverageRequest,
+    tx: GeoPoint,
+    rx: GeoPoint,
+    db: Session | None,
+) -> float:
+    if not request.include_buildings and not request.include_canopy:
+        return 0.0
+
+    points = sample_path(tx, rx, 64)
+    surface_profile = surface_profile_from_points(
+        points,
+        include_buildings=request.include_buildings,
+        include_canopy=request.include_canopy,
+        db=db,
+    )
+    surface_loss = bullington_equivalent_loss(
+        surface_profile,
+        tx_height_agl=request.tx.height_agl,
+        rx_height_agl=request.rx.height_agl,
+        freq_hz=request.tx.frequency_mhz * 1_000_000,
+    )
+    if request.model == "fspl":
+        return surface_loss
+
+    terrain_profile = surface_profile_from_points(
+        points,
+        include_buildings=False,
+        include_canopy=False,
+        db=None,
+    )
+    terrain_loss = bullington_equivalent_loss(
+        terrain_profile,
+        tx_height_agl=request.tx.height_agl,
+        rx_height_agl=request.rx.height_agl,
+        freq_hz=request.tx.frequency_mhz * 1_000_000,
+    )
+    return max(0.0, surface_loss - terrain_loss)
+
+
+def _utm_epsg(lon: float, lat: float) -> int:
+    zone = int((lon + 180) // 6) + 1
+    return 32600 + zone if lat >= 0 else 32700 + zone
+
+
+def _rx_power_at_point(
+    request: CoverageRequest,
+    lon: float,
+    lat: float,
+    distance_m: float,
+    azimuth_deg: float,
+    dem_path: str | None,
+    db: Session | None,
+) -> float:
+    distance_km = max(distance_m, 1.0) / 1000.0
+    elevation_profile = None
+    if request.model in {"itm", "p1812"} and dem_path is not None:
+        tx = GeoPoint(lat=request.tx.lat, lon=request.tx.lon)
+        rx = GeoPoint(lat=lat, lon=lon)
+        try:
+            from app.core.dem import elevations_along
+
+            points = sample_path(tx, rx, 64)
+            elevation_profile = elevations_along(points, dem_path=dem_path).tolist()
+        except Exception:
+            elevation_profile = None
+    rx = GeoPoint(lat=lat, lon=lon)
+    return _rx_power_dbm(request, distance_km, azimuth_deg, elevation_profile, db=db, rx=rx)
+
+
+def _coverage_grid(
+    request: CoverageRequest,
+    *,
+    dem_path: str | None = None,
+    db: Session | None = None,
+) -> tuple[np.ndarray, dict[str, Any]]:
+    resolution_m = request.range_step_m
+    radius_m = request.radius_m
+    width = max(1, int(np.ceil((radius_m * 2) / resolution_m)))
+    height = width
+    epsg = _utm_epsg(request.tx.lon, request.tx.lat)
+    to_utm = Transformer.from_crs("EPSG:4326", f"EPSG:{epsg}", always_xy=True)
+    to_wgs84 = Transformer.from_crs(f"EPSG:{epsg}", "EPSG:4326", always_xy=True)
+    tx_x, tx_y = to_utm.transform(request.tx.lon, request.tx.lat)
+    west = tx_x - (width * resolution_m) / 2
+    north = tx_y + (height * resolution_m) / 2
+    data = np.full((height, width), _NODATA_DBM, dtype="float32")
+
+    for row in range(height):
+        y = north - (row + 0.5) * resolution_m
+        for col in range(width):
+            x = west + (col + 0.5) * resolution_m
+            dx = x - tx_x
+            dy = y - tx_y
+            distance_m = float(np.hypot(dx, dy))
+            if distance_m > radius_m:
+                continue
+            lon, lat = to_wgs84.transform(x, y)
+            azimuth, _, geodesic_distance_m = _GEOD.inv(
+                request.tx.lon,
+                request.tx.lat,
+                lon,
+                lat,
+            )
+            data[row, col] = _rx_power_at_point(
+                request,
+                lon,
+                lat,
+                geodesic_distance_m,
+                azimuth % 360,
+                dem_path,
+                db,
+            )
+
+    transform = from_origin(west, north, resolution_m, resolution_m)
+    finite = data[data != _NODATA_DBM]
+    metadata: dict[str, Any] = {
+        "crs": f"EPSG:{epsg}",
+        "width": width,
+        "height": height,
+        "resolution_m": resolution_m,
+        "radius_m": radius_m,
+        "nodata": _NODATA_DBM,
+        "bounds": {
+            "west": west,
+            "south": north - height * resolution_m,
+            "east": west + width * resolution_m,
+            "north": north,
+        },
+        "transform": [transform.a, transform.b, transform.c, transform.d, transform.e, transform.f],
+        "value_units": "dBm",
+        "valid_pixels": int(finite.size),
+        "include_buildings": request.include_buildings,
+        "include_canopy": request.include_canopy,
+        "surface_obstruction": request.include_buildings or request.include_canopy,
+    }
+    if finite.size:
+        metadata["min_dbm"] = float(np.min(finite))
+        metadata["max_dbm"] = float(np.max(finite))
+    return data, metadata
+
+
+def _write_geotiff(
+    data: np.ndarray,
+    metadata: dict[str, Any],
+    path: Path,
+) -> None:
+    transform = Affine(*metadata["transform"])
+    with rasterio.open(
+        path,
+        "w",
+        driver="GTiff",
+        height=data.shape[0],
+        width=data.shape[1],
+        count=1,
+        dtype="float32",
+        crs=metadata["crs"],
+        transform=transform,
+        nodata=_NODATA_DBM,
+        compress="deflate",
+    ) as dataset:
+        dataset.write(data, 1)
+        dataset.update_tags(
+            model=metadata["model"],
+            frequency_mhz=str(metadata["frequency_mhz"]),
+            value_units="dBm",
+        )
+
+
+def _write_png_preview(data: np.ndarray, metadata: dict[str, Any], path: Path) -> dict[str, Any]:
+    finite_mask = data != _NODATA_DBM
+    preview = np.zeros(data.shape, dtype="uint8")
+    if finite_mask.any():
+        finite = data[finite_mask]
+        min_value = float(np.min(finite))
+        max_value = float(np.max(finite))
+        if max_value > min_value:
+            scaled = 1 + ((data[finite_mask] - min_value) / (max_value - min_value) * 254)
+            preview[finite_mask] = scaled.astype("uint8")
+        else:
+            preview[finite_mask] = 255
+    else:
+        min_value = None
+        max_value = None
+
+    with rasterio.open(
+        path,
+        "w",
+        driver="PNG",
+        height=preview.shape[0],
+        width=preview.shape[1],
+        count=1,
+        dtype="uint8",
+        crs=metadata["crs"],
+        transform=Affine(*metadata["transform"]),
+    ) as dataset:
+        dataset.write(preview, 1)
+    return {"png_min_dbm": min_value, "png_max_dbm": max_value}
+
+
+def _raster_export(
+    request: CoverageRequest,
+    *,
+    dem_path: str | None,
+    output_dir: Path,
+    db: Session | None,
+) -> dict[str, Any]:
+    output_dir.mkdir(parents=True, exist_ok=True)
+    data, metadata = _coverage_grid(request, dem_path=dem_path, db=db)
+    metadata.update(
+        {
+            "model": request.model,
+            "frequency_mhz": request.tx.frequency_mhz,
+            "tx": {"lat": request.tx.lat, "lon": request.tx.lon},
+            "format": request.format,
+        }
+    )
+    stem = f"coverage-{request.model}-{uuid4().hex[:12]}"
+    if request.format == "geotiff":
+        path = output_dir / f"{stem}.tif"
+        _write_geotiff(data, metadata, path)
+    elif request.format == "png":
+        path = output_dir / f"{stem}.png"
+        metadata.update(_write_png_preview(data, metadata, path))
+    else:
+        raise ValueError(f"Unsupported raster coverage format: {request.format}")
+
+    return {
+        "kind": "coverage_raster",
+        "format": request.format,
+        "uri": str(path),
+        "metadata": metadata,
+    }
 
 
 def _contour_points(
     request: CoverageRequest,
     level_dbm: float,
     dem_path: str | None = None,
+    db: Session | None = None,
 ) -> list[list[float]]:
     coords: list[list[float]] = []
     azimuth = 0.0
@@ -111,7 +362,7 @@ def _contour_points(
                     elevation_profile = elevations_along(points, dem_path=dem_path).tolist()
                 except Exception:
                     elevation_profile = None
-            rx_power = _rx_power_dbm(request, distance_km, azimuth, elevation_profile)
+            rx_power = _rx_power_dbm(request, distance_km, azimuth, elevation_profile, db=db)
             if rx_power >= level_dbm:
                 lon, lat, _ = _GEOD.fwd(request.tx.lon, request.tx.lat, azimuth, distance_m)
                 last_good = [lon, lat]
@@ -130,10 +381,12 @@ def compute_coverage(
     request: CoverageRequest,
     *,
     dem_path: str | None = None,
+    output_dir: str | Path | None = None,
+    db: Session | None = None,
 ) -> dict[str, Any]:
     features: list[dict[str, Any]] = []
     for level in request.levels_dbm:
-        ring = _contour_points(request, level, dem_path=dem_path)
+        ring = _contour_points(request, level, dem_path=dem_path, db=db)
         if len(ring) < 4:
             continue
         features.append(
@@ -143,6 +396,8 @@ def compute_coverage(
                     "level_dbm": level,
                     "model": request.model,
                     "frequency_mhz": request.tx.frequency_mhz,
+                    "include_buildings": request.include_buildings,
+                    "include_canopy": request.include_canopy,
                 },
                 "geometry": {"type": "Polygon", "coordinates": [ring]},
             }
@@ -156,11 +411,12 @@ def compute_coverage(
                 "model": request.model,
                 "radius_m": request.radius_m,
                 "levels_dbm": request.levels_dbm,
+                "include_buildings": request.include_buildings,
+                "include_canopy": request.include_canopy,
+                "surface_obstruction": request.include_buildings or request.include_canopy,
             },
         }
 
-    return {
-        "format": request.format,
-        "feature_count": len(features),
-        "message": "Raster coverage export is not implemented; use format=geojson",
-    }
+    if output_dir is None:
+        raise ValueError("output_dir is required for raster coverage export")
+    return _raster_export(request, dem_path=dem_path, output_dir=Path(output_dir), db=db)