RadioPropagationApi/api/app/services/terrain.py

from __future__ import annotations

import numpy as np
from sqlalchemy.exc import SQLAlchemyError
from sqlalchemy.orm import Session

from app.config import get_settings
from app.core import dem
from app.core.diffraction import deygout
from app.core.fresnel import LosSample, los_analysis
from app.core.geo import GeoPoint, PathPoint, linestring_geojson, sample_path
from app.core.surface import SurfaceProfile, build_surface_profile
from app.models.terrain import (
    ElevationResponse,
    LosRequest,
    LosResponse,
    Obstruction,
    ProfileSample,
    TerrainProfileRequest,
    TerrainProfileResponse,
)
from app.services import buildings as buildings_service
from app.services.vegetation import vegetation_loss_along


def elevation_at(lat: float, lon: float, surface: str) -> ElevationResponse:
    settings = get_settings()
    elevation_m = dem.elevation_at(lat, lon, surface=surface, dem_path=settings.dem_path)
    return ElevationResponse(lat=lat, lon=lon, elevation_m=elevation_m, surface=surface)


def terrain_profile(
    request: TerrainProfileRequest,
    db: Session | None = None,
) -> TerrainProfileResponse:
    points = sample_path(
        GeoPoint(lat=request.start.lat, lon=request.start.lon),
        GeoPoint(lat=request.end.lat, lon=request.end.lon),
        request.samples,
    )
    profile = surface_profile_from_points(
        points,
        include_buildings=request.include_buildings,
        include_canopy=request.include_canopy,
        db=db,
    )
    return TerrainProfileResponse(
        distance_m=profile.distance_m,
        samples=[ProfileSample(**sample.__dict__) for sample in profile.samples],
        path_geojson=linestring_geojson(points),
    )


def los(
    request: LosRequest,
    db: Session | None = None,
) -> LosResponse:
    profile_request = TerrainProfileRequest(
        start=request.tx,
        end=request.rx,
        samples=request.samples,
        include_buildings=request.include_buildings,
        include_canopy=request.include_canopy,
        k_factor=request.k_factor,
    )
    points = sample_path(
        GeoPoint(lat=request.tx.lat, lon=request.tx.lon),
        GeoPoint(lat=request.rx.lat, lon=request.rx.lon),
        request.samples,
    )
    surface_profile = surface_profile_from_points(
        points,
        include_buildings=request.include_buildings,
        include_canopy=request.include_canopy,
        db=db,
    )
    result = los_analysis(
        surface_profile,
        request.tx.height_agl,
        request.rx.height_agl,
        request.frequency_mhz * 1_000_000,
        request.fresnel_clearance,
        request.k_factor,
    )
    diffraction_loss = deygout(
        surface_profile,
        request.tx.height_agl,
        request.rx.height_agl,
        request.frequency_mhz * 1_000_000,
    )
    vegetation_loss = vegetation_loss_along(
        points,
        freq_hz=request.frequency_mhz * 1_000_000,
        include_vegetation=request.include_vegetation,
    )

    def convert(sample: LosSample) -> Obstruction:
        return Obstruction(
            distance_m=sample.distance_m,
            clearance_m=sample.clearance_m,
            fresnel_radius_m=sample.fresnel_radius_m,
            required_clearance_m=sample.required_clearance_m,
            type=sample.obstruction_type,
        )

    geometric_obstructions = [
        sample for sample in result.obstructions if sample.clearance_m < 0
    ]
    building_obstructions = [
        sample
        for sample in geometric_obstructions
        if sample.obstruction_type == "building"
    ]

    return LosResponse(
        los_clear=result.los_clear,
        geometric_los=result.geometric_los,
        first_fresnel_clearance_pct=result.first_fresnel_clearance_pct,
        worst_obstruction=convert(result.worst_sample) if result.worst_sample else None,
        obstructions=[convert(sample) for sample in result.obstructions],
        geometric_obstructions=[convert(sample) for sample in geometric_obstructions],
        fresnel_violations_count=len(result.obstructions),
        geometric_obstructions_count=len(geometric_obstructions),
        building_obstructions_count=len(building_obstructions),
        diffraction_loss_db=diffraction_loss,
        vegetation_loss_db=vegetation_loss,
        profile_ref=f"inline:{profile_request.samples}",
    )


def _building_heights(
    points: list[PathPoint],
    include_buildings: bool,
    db: Session | None,
) -> list[float] | None:
    if not include_buildings or db is None:
        return None
    try:
        return buildings_service.building_heights_along(points, db)
    except SQLAlchemyError:
        return None


def surface_profile_from_points(
    points: list[PathPoint],
    include_buildings: bool,
    include_canopy: bool,
    db: Session | None = None,
) -> SurfaceProfile:
    try:
        ground_elevations = dem.elevations_along(points, dem_path=get_settings().dem_path).tolist()
    except dem.DemNotConfiguredError:
        ground_elevations = np.zeros(len(points), dtype=float).tolist()

    building_heights = _building_heights(points, include_buildings, db)
    return build_surface_profile(
        points,
        ground_elevations=ground_elevations,
        building_heights=building_heights,
        include_buildings=include_buildings,
        include_canopy=include_canopy,
    )