added formulas

API.md

@@ -7,8 +7,8 @@
 
 ## Состояние Реализации
 
-- Работает: healthcheck, OpenAPI, DEM elevation/profile, OSM buildings query, landcover/canopy path sampling, terrain LOS/Fresnel/diffraction с учётом DEM и buildings, link budget (`manual`, `itm`, `p452`) с P.676 atmospheric loss и P.833 vegetation, antenna omni/sector/file, async coverage/viewshed jobs через Redis + Celery.
+- Работает: healthcheck, OpenAPI, DEM elevation/profile, OSM buildings query, landcover/canopy path sampling, terrain LOS/Fresnel/diffraction с учётом DEM, buildings и canopy surface layer, link budget (`manual`, `itm`, `p452`) с P.676 atmospheric loss, P.833 vegetation и canopy препятствиями, antenna omni/sector/file, async coverage/viewshed jobs через Redis + Celery.
-- Частично работает: coverage raster export (`format=geotiff|png`), canopy bootstrap script (landcover загружается, canopy — вручную).
+- Частично работает: coverage raster export (`format=geotiff|png`), canopy bootstrap script (данные можно качать автоматически, но объём большой).
 - Пока не реализовано: полный P.1812 с ITU digital maps (используется terrain-based ITM + clutter), WhiteboxTools viewshed (используется `gdal_viewshed`).
 
 ## Общие Правила
@@ -336,9 +336,10 @@ Feature properties:
 
 ### `POST /api/v1/link/budget`
 
-Считает point-to-point budget. Сейчас рабочий режим: `model="manual"`.
+Считает point-to-point budget. Рабочие модели: `manual`, `itm`, `p452`.
-Он учитывает FSPL, terrain/buildings diffraction по DEM/PostGIS profile,
+`manual` использует FSPL + Bullington/P.526 + P.676 + vegetation. `itm`
-и vegetation attenuation по WorldCover при `include_vegetation=true`.
+использует Longley-Rice ITM + P.676 + vegetation. `p452` использует pycraf
+P.452 + vegetation.
 
 ```bash
 curl -s -X POST http://localhost:5603/api/v1/link/budget \
@@ -359,6 +360,9 @@ curl -s -X POST http://localhost:5603/api/v1/link/budget \
 - `distance_km`
 - `fspl_db`
 - `diffraction_db`
+- `propagation_model_loss_db`
+- `propagation_excess_db`
+- `knife_edge_diffraction_db`
 - `vegetation_db`
 - `atmospheric_db`
 - `total_loss_db`
@@ -379,9 +383,12 @@ curl -s -X POST http://localhost:5603/api/v1/link/budget \
 P.833 vegetation attenuation в `vegetation_db`. Если landcover raster отсутствует,
 значение остаётся `0`, чтобы link budget продолжал работать.
 
-При `include_buildings=true` API также строит terrain surface profile с DEM и
+API строит terrain surface profile с DEM, OSM buildings и canopy, считает
-OSM buildings, считает Fresnel/LOS и добавляет Bullington-style diffraction loss
+Fresnel/LOS и Bullington-style diffraction diagnostic loss. Для `model=manual`
-в `diffraction_db`. Поле `fresnel_clear` берётся из этого же анализа.
+`diffraction_db` равен `knife_edge_diffraction_db`. Для `model=itm|p452`
+`diffraction_db` оставлен как active propagation excess над FSPL; для новых
+интеграций используй явные поля `propagation_excess_db` и
+`knife_edge_diffraction_db`.
 
 `link_viable` считается по фактическому `fade_margin_db > 0` после применённых
 потерь. Состояние зоны Френеля не выключает линк напрямую, а отдаётся отдельным
@@ -394,8 +401,11 @@ OSM buildings, считает Fresnel/LOS и добавляет Bullington-style
 Диагностические поля `geometric_los`, `worst_obstruction` и счётчики obstruction
 помогают понять, почему качество Френеля ухудшилось.
 
-Модели `p452` и `itm` используют ITU-R P.452 (pycraf) и Longley-Rice ITM (itmlogic) поверх DEM profile.
+Модели `p452` и `itm` используют ITU-R P.452 (pycraf) и Longley-Rice ITM
-Поле `atmospheric_db` считается по ITU-R P.676 (itur). Флаг `include_canopy` учитывает canopy raster в surface profile.
+(itmlogic) поверх DEM profile. Поле `atmospheric_db` считается по ITU-R P.676
+(itur) для `manual` и `itm`; для `p452` оно равно `0`, чтобы не дублировать
+атмосферные потери, уже входящие в P.452. Флаг `include_canopy` учитывает canopy
+raster в surface profile.
 
 ## Antenna
 

api/app/core/diffraction.py

@@ -2,7 +2,7 @@ from __future__ import annotations
 
 from math import log10, sqrt
 
-from app.core.fresnel import wavelength
+from app.core.fresnel import earth_bulge, wavelength
 from app.core.surface import SurfaceProfile
 
 
@@ -24,6 +24,7 @@ def bullington_loss(
     tx_height_agl: float,
     rx_height_agl: float,
     freq_hz: float,
+    k: float = 1.333,
 ) -> float:
     """Bullington-style equivalent edge loss for a terrain profile.
 
@@ -45,17 +46,34 @@ def bullington_loss(
         d1 = sample.distance_m
         d2 = total_distance - d1
         path_height = tx_elevation + (rx_elevation - tx_elevation) * (d1 / total_distance)
-        h = sample.surface_m - path_height
+        h = sample.surface_m + earth_bulge(d1, d2, k) - path_height
         max_v = max(max_v, knife_edge_v(h, d1, d2, freq_hz))
 
     return knife_edge_loss(max_v)
 
 
+def bullington_equivalent_loss(
+    profile: SurfaceProfile,
+    tx_height_agl: float,
+    rx_height_agl: float,
+    freq_hz: float,
+    k: float = 1.333,
+) -> float:
+    return bullington_loss(profile, tx_height_agl, rx_height_agl, freq_hz, k=k)
+
+
 def deygout(
     profile: SurfaceProfile,
     tx_height_agl: float,
     rx_height_agl: float,
     freq_hz: float,
+    k: float = 1.333,
 ) -> float:
-    """Compatibility wrapper; use Bullington equivalent loss for multi-edge profiles."""
+    """Compatibility wrapper; this is Bullington-style, not multi-edge Deygout."""
-    return bullington_loss(profile, tx_height_agl, rx_height_agl, freq_hz)
+    return bullington_equivalent_loss(
+        profile,
+        tx_height_agl,
+        rx_height_agl,
+        freq_hz,
+        k=k,
+    )

api/app/models/link.py

@@ -37,4 +37,7 @@ class LinkBudgetResponse(BaseModel):
     geometric_obstructions_count: int
     building_obstructions_count: int
     worst_obstruction: Obstruction | None
+    propagation_model_loss_db: float
+    propagation_excess_db: float
+    knife_edge_diffraction_db: float
     data_sources: DataSources = Field(default_factory=DataSources)

api/app/services/link.py

@@ -2,7 +2,7 @@ from sqlalchemy.orm import Session
 
 from app.config import get_settings
 from app.core.atmosphere import atmospheric_loss_between
-from app.core.diffraction import deygout
+from app.core.diffraction import bullington_equivalent_loss
 from app.core.fresnel import LosSample, los_analysis
 from app.core.geo import GeoPoint, sample_path
 from app.core.propagation import itm_loss, manual_link_budget, p452_loss
@@ -34,11 +34,12 @@ def link_budget(request: LinkBudgetRequest, db: Session | None = None) -> LinkBu
         freq_hz=freq_hz,
         k=request.k_factor,
     )
-    diffraction_db = deygout(
+    knife_edge_diffraction_db = bullington_equivalent_loss(
         surface_profile,
         tx_height_agl=request.tx.height_agl,
         rx_height_agl=request.rx.height_agl,
         freq_hz=freq_hz,
+        k=request.k_factor,
     )
     vegetation_db = vegetation_loss_along(
         points,
@@ -47,55 +48,46 @@ def link_budget(request: LinkBudgetRequest, db: Session | None = None) -> LinkBu
     )
     atmospheric_db = atmospheric_loss_between(tx, rx, freq_hz)
 
-    active_diffraction_db = diffraction_db
+    fspl_db = manual_link_budget(
+        tx,
+        rx,
+        request.tx.power_dbm,
+        request.tx.gain_dbi,
+        request.rx.gain_dbi,
+        request.rx.sensitivity_dbm,
+        request.frequency_mhz,
+    ).fspl_db
+    active_diffraction_db = knife_edge_diffraction_db
     active_atmospheric_db = atmospheric_db
+    propagation_model_loss_db = fspl_db + knife_edge_diffraction_db
+    propagation_excess_db = knife_edge_diffraction_db
 
     if request.model == "manual":
-        propagation_db = 0.0
+        pass
     elif request.model == "itm":
-        propagation_db = (
+        propagation_model_loss_db = itm_loss(
-            itm_loss(
+            tx,
-                tx,
+            rx,
-                rx,
+            tx_height_agl=request.tx.height_agl,
-                tx_height_agl=request.tx.height_agl,
+            rx_height_agl=request.rx.height_agl,
-                rx_height_agl=request.rx.height_agl,
+            freq_mhz=request.frequency_mhz,
-                freq_mhz=request.frequency_mhz,
+            elevation_profile_m=elevation_profile,
-                elevation_profile_m=elevation_profile,
+            dem_path=str(settings.dem_path),
-                dem_path=str(settings.dem_path),
-            )
-            - manual_link_budget(
-                tx,
-                rx,
-                request.tx.power_dbm,
-                request.tx.gain_dbi,
-                request.rx.gain_dbi,
-                request.rx.sensitivity_dbm,
-                request.frequency_mhz,
-            ).fspl_db
         )
-        active_diffraction_db = propagation_db
+        propagation_excess_db = propagation_model_loss_db - fspl_db
+        active_diffraction_db = propagation_excess_db
     elif request.model == "p452":
-        propagation_db = (
+        propagation_model_loss_db = p452_loss(
-            p452_loss(
+            tx,
-                tx,
+            rx,
-                rx,
+            tx_height_agl=request.tx.height_agl,
-                tx_height_agl=request.tx.height_agl,
+            rx_height_agl=request.rx.height_agl,
-                rx_height_agl=request.rx.height_agl,
+            freq_mhz=request.frequency_mhz,
-                freq_mhz=request.frequency_mhz,
+            elevation_profile_m=elevation_profile,
-                elevation_profile_m=elevation_profile,
+            dem_path=str(settings.dem_path),
-                dem_path=str(settings.dem_path),
-            )
-            - manual_link_budget(
-                tx,
-                rx,
-                request.tx.power_dbm,
-                request.tx.gain_dbi,
-                request.rx.gain_dbi,
-                request.rx.sensitivity_dbm,
-                request.frequency_mhz,
-            ).fspl_db
         )
-        active_diffraction_db = propagation_db
+        propagation_excess_db = propagation_model_loss_db - fspl_db
+        active_diffraction_db = propagation_excess_db
         active_atmospheric_db = 0.0
     else:
         raise NotImplementedError(f"{request.model} link model is not supported")
@@ -131,6 +123,9 @@ def link_budget(request: LinkBudgetRequest, db: Session | None = None) -> LinkBu
         geometric_obstructions_count=len(geometric_obstructions),
         building_obstructions_count=len(building_obstructions),
         worst_obstruction=_convert_obstruction(los_result.worst_sample),
+        propagation_model_loss_db=propagation_model_loss_db,
+        propagation_excess_db=propagation_excess_db,
+        knife_edge_diffraction_db=knife_edge_diffraction_db,
         data_sources=DataSources(
             dem=bool(raster_files(settings.dem_path)),
             buildings=request.include_buildings and db is not None,

api/app/services/terrain.py

@@ -10,10 +10,10 @@ 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.diffraction import bullington_equivalent_loss
 from app.core.fresnel import LosSample, los_analysis
 from app.core.geo import GeoPoint, PathPoint, linestring_geojson, sample_path
-from app.core.raster_sampling import raster_files
+from app.core.raster_sampling import RasterNotConfiguredError, raster_files, sample_along
 from app.core.surface import SurfaceProfile, build_surface_profile
 from app.models.common import DataSources
 from app.models.terrain import (
@@ -67,6 +67,8 @@ def terrain_profile(
             dem_used=bool(raster_files(get_settings().dem_path)),
             include_buildings=request.include_buildings,
             db=db,
+            canopy_used=request.include_canopy
+            and bool(raster_files(get_settings().canopy_path)),
         ),
     )
 
@@ -103,11 +105,12 @@ def los(
         request.fresnel_clearance,
         request.k_factor,
     )
-    diffraction_loss = deygout(
+    diffraction_loss = bullington_equivalent_loss(
         surface_profile,
         request.tx.height_agl,
         request.rx.height_agl,
         request.frequency_mhz * 1_000_000,
+        k=request.k_factor,
     )
     vegetation_loss = vegetation_loss_along(
         points,
@@ -247,6 +250,24 @@ def _building_heights(
         return None
 
 
+def _canopy_heights(points: list[PathPoint], include_canopy: bool) -> list[float] | None:
+    if not include_canopy:
+        return None
+    try:
+        values = sample_along(
+            points,
+            get_settings().canopy_path,
+            "canopy",
+            require_all=False,
+        )
+    except RasterNotConfiguredError:
+        return None
+
+    values = np.nan_to_num(values, nan=0.0, posinf=0.0, neginf=0.0)
+    values = np.maximum(values, 0.0)
+    return values.tolist()
+
+
 def surface_profile_from_points(
     points: list[PathPoint],
     include_buildings: bool,
@@ -259,10 +280,12 @@ def surface_profile_from_points(
         ground_elevations = np.zeros(len(points), dtype=float).tolist()
 
     building_heights = _building_heights(points, include_buildings, db)
+    canopy_heights = _canopy_heights(points, include_canopy)
     return build_surface_profile(
         points,
         ground_elevations=ground_elevations,
         building_heights=building_heights,
+        canopy_heights=canopy_heights,
         include_buildings=include_buildings,
         include_canopy=include_canopy,
     )

api/tests/test_api.py

@@ -147,4 +147,7 @@ def test_manual_link_budget_endpoint() -> None:
     assert data["fspl_db"] > 0
     assert "link_viable" in data
     assert data["fresnel_quality"] in {"clear", "partial", "blocked"}
+    assert "propagation_model_loss_db" in data
+    assert "propagation_excess_db" in data
+    assert "knife_edge_diffraction_db" in data
     assert "data_sources" in data

api/tests/test_core.py

@@ -71,7 +71,7 @@ def test_bullington_loss_uses_dominant_equivalent_edge() -> None:
             ),
         ],
     )
-    v = knife_edge_v(20, 5_000, 5_000, 433_000_000)
+    v = knife_edge_v(20 + earth_bulge(5_000, 5_000, 1.333), 5_000, 5_000, 433_000_000)
 
     assert isclose(bullington_loss(profile, 0, 0, 433_000_000), knife_edge_loss(v))
 

api/tests/test_vegetation_integration.py

@@ -94,6 +94,9 @@ def test_link_budget_includes_diffraction_loss(monkeypatch) -> None:
     result = link_service.link_budget(request)
 
     assert result.diffraction_db > 0
+    assert result.knife_edge_diffraction_db == result.diffraction_db
+    assert result.propagation_excess_db == result.diffraction_db
+    assert isclose(result.propagation_model_loss_db, result.fspl_db + result.diffraction_db)
     assert isclose(
         result.total_loss_db,
         result.fspl_db + result.diffraction_db + result.atmospheric_db,
@@ -156,7 +159,7 @@ def test_link_viable_depends_on_margin_not_fresnel_quality(monkeypatch) -> None:
 def test_itm_link_budget_does_not_double_count_diffraction(monkeypatch) -> None:
     monkeypatch.setattr(link_service, "vegetation_loss_along", lambda *args, **kwargs: 3.0)
     monkeypatch.setattr(link_service, "atmospheric_loss_between", lambda *args, **kwargs: 2.0)
-    monkeypatch.setattr(link_service, "deygout", lambda *args, **kwargs: 99.0)
+    monkeypatch.setattr(link_service, "bullington_equivalent_loss", lambda *args, **kwargs: 99.0)
     monkeypatch.setattr(
         link_service,
         "surface_profile_from_points",
@@ -199,6 +202,9 @@ def test_itm_link_budget_does_not_double_count_diffraction(monkeypatch) -> None:
     result = link_service.link_budget(request)
 
     assert result.diffraction_db == 12.0
+    assert result.knife_edge_diffraction_db == 99.0
+    assert result.propagation_model_loss_db == result.fspl_db + 12.0
+    assert result.propagation_excess_db == 12.0
     assert result.atmospheric_db == 2.0
     assert isclose(result.total_loss_db, result.fspl_db + 12.0 + 3.0 + 2.0)
 
@@ -206,7 +212,7 @@ def test_itm_link_budget_does_not_double_count_diffraction(monkeypatch) -> None:
 def test_p452_link_budget_does_not_double_count_diffraction_or_atmosphere(monkeypatch) -> None:
     monkeypatch.setattr(link_service, "vegetation_loss_along", lambda *args, **kwargs: 3.0)
     monkeypatch.setattr(link_service, "atmospheric_loss_between", lambda *args, **kwargs: 99.0)
-    monkeypatch.setattr(link_service, "deygout", lambda *args, **kwargs: 88.0)
+    monkeypatch.setattr(link_service, "bullington_equivalent_loss", lambda *args, **kwargs: 88.0)
     monkeypatch.setattr(
         link_service,
         "surface_profile_from_points",
@@ -249,5 +255,8 @@ def test_p452_link_budget_does_not_double_count_diffraction_or_atmosphere(monkey
     result = link_service.link_budget(request)
 
     assert result.diffraction_db == 14.0
+    assert result.knife_edge_diffraction_db == 88.0
+    assert result.propagation_model_loss_db == result.fspl_db + 14.0
+    assert result.propagation_excess_db == 14.0
     assert result.atmospheric_db == 0.0
     assert isclose(result.total_loss_db, result.fspl_db + 14.0 + 3.0)

api/tests/test_viewshed_and_canopy.py

@@ -2,13 +2,15 @@ from __future__ import annotations
 
 import sys
 from pathlib import Path
+from types import SimpleNamespace
 
 import numpy as np
 import rasterio
 from rasterio.transform import from_origin
 
-from app.core.geo import GeoPoint
+from app.core.geo import GeoPoint, PathPoint
 from app.core.viewshed import _find_dem_tiles, _radius_bbox
+from app.services import terrain as terrain_service
 
 sys.path.append(str(Path(__file__).resolve().parents[2] / "scripts"))
 import bootstrap_canopy  # noqa: E402
@@ -16,6 +18,10 @@ import bootstrap_canopy  # noqa: E402
 
 def _write_dem(path: Path, west: float) -> None:
     data = np.ones((10, 10), dtype="float32")
+    _write_raster(path, west, data)
+
+
+def _write_raster(path: Path, west: float, data: np.ndarray) -> None:
     transform = from_origin(west, 60.0, 0.01, 0.01)
     with rasterio.open(
         path,
@@ -32,6 +38,10 @@ def _write_dem(path: Path, west: float) -> None:
         dataset.write(data, 1)
 
 
+def _settings(dem_path: Path, canopy_path: Path) -> SimpleNamespace:
+    return SimpleNamespace(dem_path=dem_path, canopy_path=canopy_path)
+
+
 def test_viewshed_dem_search_uses_radius_bbox(tmp_path: Path) -> None:
     _write_dem(tmp_path / "west.tif", 30.0)
     _write_dem(tmp_path / "east.tif", 30.1)
@@ -42,6 +52,79 @@ def test_viewshed_dem_search_uses_radius_bbox(tmp_path: Path) -> None:
     assert {tile.name for tile in tiles} == {"west.tif", "east.tif"}
 
 
+def test_surface_profile_uses_canopy_heights(monkeypatch, tmp_path: Path) -> None:
+    dem_dir = tmp_path / "dem"
+    canopy_dir = tmp_path / "canopy"
+    dem_dir.mkdir()
+    canopy_dir.mkdir()
+    _write_raster(dem_dir / "dem.tif", 30.0, np.full((10, 10), 2.0, dtype="float32"))
+    _write_raster(
+        canopy_dir / "canopy.tif",
+        30.0,
+        np.full((10, 10), 12.0, dtype="float32"),
+    )
+    monkeypatch.setattr(terrain_service, "get_settings", lambda: _settings(dem_dir, canopy_dir))
+    points = [
+        PathPoint(lat=59.95, lon=30.05, distance_m=0),
+        PathPoint(lat=59.94, lon=30.06, distance_m=1_000),
+    ]
+
+    profile = terrain_service.surface_profile_from_points(
+        points,
+        include_buildings=False,
+        include_canopy=True,
+    )
+
+    assert [sample.ground_m for sample in profile.samples] == [2.0, 2.0]
+    assert [sample.canopy_m for sample in profile.samples] == [12.0, 12.0]
+    assert [sample.surface_m for sample in profile.samples] == [14.0, 14.0]
+
+
+def test_surface_profile_ignores_canopy_when_disabled(monkeypatch, tmp_path: Path) -> None:
+    dem_dir = tmp_path / "dem"
+    canopy_dir = tmp_path / "canopy"
+    dem_dir.mkdir()
+    canopy_dir.mkdir()
+    _write_raster(dem_dir / "dem.tif", 30.0, np.full((10, 10), 2.0, dtype="float32"))
+    _write_raster(
+        canopy_dir / "canopy.tif",
+        30.0,
+        np.full((10, 10), 12.0, dtype="float32"),
+    )
+    monkeypatch.setattr(terrain_service, "get_settings", lambda: _settings(dem_dir, canopy_dir))
+    points = [PathPoint(lat=59.95, lon=30.05, distance_m=0)]
+
+    profile = terrain_service.surface_profile_from_points(
+        points,
+        include_buildings=False,
+        include_canopy=False,
+    )
+
+    assert profile.samples[0].canopy_m == 0.0
+    assert profile.samples[0].surface_m == 2.0
+
+
+def test_surface_profile_treats_missing_canopy_as_optional(
+    monkeypatch,
+    tmp_path: Path,
+) -> None:
+    dem_dir = tmp_path / "dem"
+    canopy_dir = tmp_path / "canopy"
+    dem_dir.mkdir()
+    _write_raster(dem_dir / "dem.tif", 30.0, np.full((10, 10), 2.0, dtype="float32"))
+    monkeypatch.setattr(terrain_service, "get_settings", lambda: _settings(dem_dir, canopy_dir))
+    points = [PathPoint(lat=59.95, lon=30.05, distance_m=0)]
+
+    profile = terrain_service.surface_profile_from_points(
+        points,
+        include_buildings=False,
+        include_canopy=True,
+    )
+
+    assert profile.samples[0].canopy_m == 0.0
+    assert profile.samples[0].surface_m == 2.0
+
+
 def test_canopy_selects_intersecting_tiles_from_geojson_index() -> None:
     index = {
         "type": "FeatureCollection",