From f4f254116b3c218ae4fc3850f7f6b6bb8b7e9b13 Mon Sep 17 00:00:00 2001 From: grigo Date: Wed, 17 Jun 2026 13:26:35 +0300 Subject: [PATCH] first commit --- .gitignore | 12 +++ Dockerfile | 22 +++++ README.md | 118 ++++++++++++++++++++++++++ app/__init__.py | 0 app/config.py | 19 +++++ app/main.py | 46 ++++++++++ app/models/__init__.py | 0 app/models/schemas.py | 118 ++++++++++++++++++++++++++ app/routers/__init__.py | 0 app/routers/coverage.py | 27 ++++++ app/routers/elevation.py | 36 ++++++++ app/routers/link.py | 30 +++++++ app/services/__init__.py | 0 app/services/buildings.py | 90 ++++++++++++++++++++ app/services/canopy.py | 16 ++++ app/services/coverage.py | 70 ++++++++++++++++ app/services/fresnel.py | 19 +++++ app/services/geodesy.py | 73 ++++++++++++++++ app/services/link_calculator.py | 144 ++++++++++++++++++++++++++++++++ app/services/path_loss.py | 36 ++++++++ app/services/raster_base.py | 121 +++++++++++++++++++++++++++ app/services/terrain.py | 16 ++++ data/canopy/.gitkeep | 0 data/osm/.gitkeep | 0 data/srtm/.gitkeep | 0 docker-compose.yml | 36 ++++++++ requirements.txt | 15 ++++ scripts/download_canopy.py | 68 +++++++++++++++ scripts/download_srtm.py | 80 ++++++++++++++++++ scripts/import_buildings.py | 136 ++++++++++++++++++++++++++++++ sql/init.sql | 10 +++ 31 files changed, 1358 insertions(+) create mode 100644 .gitignore create mode 100644 Dockerfile create mode 100644 README.md create mode 100644 app/__init__.py create mode 100644 app/config.py create mode 100644 app/main.py create mode 100644 app/models/__init__.py create mode 100644 app/models/schemas.py create mode 100644 app/routers/__init__.py create mode 100644 app/routers/coverage.py create mode 100644 app/routers/elevation.py create mode 100644 app/routers/link.py create mode 100644 app/services/__init__.py create mode 100644 app/services/buildings.py create mode 100644 app/services/canopy.py create mode 100644 app/services/coverage.py create mode 100644 app/services/fresnel.py create mode 100644 app/services/geodesy.py create mode 100644 app/services/link_calculator.py create mode 100644 app/services/path_loss.py create mode 100644 app/services/raster_base.py create mode 100644 app/services/terrain.py create mode 100644 data/canopy/.gitkeep create mode 100644 data/osm/.gitkeep create mode 100644 data/srtm/.gitkeep create mode 100644 docker-compose.yml create mode 100644 requirements.txt create mode 100644 scripts/download_canopy.py create mode 100644 scripts/download_srtm.py create mode 100644 scripts/import_buildings.py create mode 100644 sql/init.sql diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..27fcc44 --- /dev/null +++ b/.gitignore @@ -0,0 +1,12 @@ +__pycache__/ +*.py[cod] +.env +.venv/ +venv/ +*.egg-info/ +.pytest_cache/ +data/srtm/*.tif +data/srtm/*.tiff +data/canopy/*.tif +data/canopy/*.tiff +data/osm/*.pbf diff --git a/Dockerfile b/Dockerfile new file mode 100644 index 0000000..a80b2e9 --- /dev/null +++ b/Dockerfile @@ -0,0 +1,22 @@ +FROM python:3.12-slim + +RUN apt-get update && apt-get install -y --no-install-recommends \ + gdal-bin \ + libgdal-dev \ + gcc \ + g++ \ + && rm -rf /var/lib/apt/lists/* + +ENV GDAL_CONFIG=/usr/bin/gdal-config + +WORKDIR /app + +COPY requirements.txt . +RUN pip install --no-cache-dir -r requirements.txt + +COPY app/ ./app/ +COPY sql/ ./sql/ + +EXPOSE 8000 + +CMD ["uvicorn", "app.main:app", "--host", "0.0.0.0", "--port", "8000"] diff --git a/README.md b/README.md new file mode 100644 index 0000000..145a3e9 --- /dev/null +++ b/README.md @@ -0,0 +1,118 @@ +# Radio Link Service + +Микросервис для расчёта радиосвязи между двумя точками с учётом рельефа (SRTM), крон деревьев (ETH Canopy) и зданий OSM (PostGIS). + +## Быстрый старт + +```bash +docker compose up -d --build +``` + +Сервис: http://localhost:8000 +Документация API: http://localhost:8000/docs + +## Подготовка данных (офлайн) + +После загрузки данных сервис работает полностью офлайн. + +### 1. Рельеф SRTM GL1 (30 м) + +```bash +python scripts/download_srtm.py --bbox 59.0 30.0 60.0 31.0 --output data/srtm +``` + +Для OpenTopography может потребоваться API key: + +```bash +python scripts/download_srtm.py --bbox 59.0 30.0 60.0 31.0 --api-key YOUR_KEY +``` + +### 2. Высота крон ETH Global Canopy Height 2020 (10 м) + +```bash +python scripts/download_canopy.py --bbox 59.0 30.0 60.0 31.0 --output data/canopy +``` + +### 3. Здания OSM + +Через bbox (онлайн, osmnx): + +```bash +python scripts/import_buildings.py --bbox 59.0 30.0 60.0 31.0 --truncate +``` + +Через локальный PBF: + +```bash +python scripts/import_buildings.py --pbf data/osm/region.osm.pbf --truncate +``` + +## Примеры запросов + +### Health + +```bash +curl http://localhost:8000/health +``` + +### Расчёт линка + +```bash +curl -X POST http://localhost:8000/v1/link \ + -H "Content-Type: application/json" \ + -d '{ + "from": {"lat": 59.93, "lng": 30.33, "antenna_height_m": 2.0}, + "to": {"lat": 59.95, "lng": 30.40, "antenna_height_m": 2.0}, + "frequency_mhz": 433.5, + "fresnel_clearance": 0.6 + }' +``` + +Опциональные параметры link budget: `tx_power_dbm` (default 20), `rx_sensitivity_dbm` (default -137). + +### Профиль высот (только рельеф) + +```bash +curl "http://localhost:8000/v1/elevation-profile?from_lat=59.93&from_lng=30.33&to_lat=59.95&to_lng=30.40&step_m=30" +``` + +### Зона покрытия + +```bash +curl -X POST http://localhost:8000/v1/coverage \ + -H "Content-Type: application/json" \ + -d '{ + "origin": {"lat": 59.93, "lng": 30.33, "antenna_height_m": 10.0}, + "frequency_mhz": 433.5, + "radius_m": 5000, + "resolution_m": 100 + }' +``` + +Лимит сетки: `MAX_GRID_POINTS=10000` (env). При превышении — HTTP 400. + +## Локальный запуск без Docker + +```bash +pip install -r requirements.txt +export DATABASE_URL=postgresql://radio:radio@localhost:5432/radio +uvicorn app.main:app --reload +``` + +## Структура + +- `app/` — FastAPI приложение +- `data/srtm/` — GeoTIFF тайлы рельефа +- `data/canopy/` — GeoTIFF тайлы крон +- `data/osm/` — PBF файлы OSM +- `scripts/` — загрузка данных и импорт зданий +- `sql/init.sql` — схема PostGIS + +## Алгоритм `/v1/link` + +1. Сэмплирование трассы с шагом 30 м (геодезическая линия) +2. Для каждой точки: SRTM + canopy + max(building height) +3. LOS-линия между антеннами +4. Радиус Френеля: `r = 17.3 * sqrt(d1*d2 / (f_ghz * d))` +5. Clearance: `los_height - total_height - r * fresnel_clearance` +6. FSPL: `20*log10(d_km) + 20*log10(f_mhz) + 27.55` diff --git a/app/__init__.py b/app/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/app/config.py b/app/config.py new file mode 100644 index 0000000..8d067f8 --- /dev/null +++ b/app/config.py @@ -0,0 +1,19 @@ +from pydantic_settings import BaseSettings, SettingsConfigDict + + +class Settings(BaseSettings): + model_config = SettingsConfigDict(env_file=".env", extra="ignore") + + database_url: str = "postgresql://radio:radio@localhost:5432/radio" + srtm_path: str = "data/srtm" + canopy_path: str = "data/canopy" + sample_step_m: float = 30.0 + default_building_height_m: float = 10.0 + default_tx_power_dbm: float = 20.0 + default_rx_sensitivity_dbm: float = -137.0 + max_grid_points: int = 10000 + coverage_concurrency: int = 16 + raster_cache_size: int = 16 + + +settings = Settings() diff --git a/app/main.py b/app/main.py new file mode 100644 index 0000000..f0d53b1 --- /dev/null +++ b/app/main.py @@ -0,0 +1,46 @@ +from __future__ import annotations + +from contextlib import asynccontextmanager + +from fastapi import FastAPI + +from app.config import settings +from app.models.schemas import DataLayersHealth, HealthResponse +from app.routers import coverage, elevation, link +from app.services.buildings import BuildingService +from app.services.canopy import CanopyService +from app.services.coverage import CoverageService +from app.services.link_calculator import LinkCalculator +from app.services.terrain import TerrainService + +terrain_service = TerrainService(settings.srtm_path, cache_size=settings.raster_cache_size) +canopy_service = CanopyService(settings.canopy_path, cache_size=settings.raster_cache_size) +building_service = BuildingService(settings.database_url) +link_calculator = LinkCalculator(terrain_service, canopy_service, building_service) +coverage_service = CoverageService(link_calculator) + + +@asynccontextmanager +async def lifespan(_: FastAPI): + building_service.connect() + yield + building_service.close() + terrain_service.close() + canopy_service.close() + + +app = FastAPI(title="Radio Link Service", version="1.0.0", lifespan=lifespan) +app.include_router(link.router) +app.include_router(coverage.router) +app.include_router(elevation.router) + + +@app.get("/health", response_model=HealthResponse) +async def health() -> HealthResponse: + layers = DataLayersHealth( + srtm=terrain_service.is_available(), + canopy=canopy_service.is_available(), + buildings=building_service.is_available(), + ) + status = "ok" if all(layers.model_dump().values()) else "degraded" + return HealthResponse(status=status, data_layers=layers) diff --git a/app/models/__init__.py b/app/models/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/app/models/schemas.py b/app/models/schemas.py new file mode 100644 index 0000000..d36bfcd --- /dev/null +++ b/app/models/schemas.py @@ -0,0 +1,118 @@ +from __future__ import annotations + +from typing import Any, Literal + +from pydantic import BaseModel, ConfigDict, Field, field_validator + + +class GeoAntenna(BaseModel): + lat: float + lng: float + antenna_height_m: float = Field(default=2.0, ge=0) + + @field_validator("lat") + @classmethod + def validate_lat(cls, value: float) -> float: + if not -90 <= value <= 90: + raise ValueError("lat must be between -90 and 90") + return value + + @field_validator("lng") + @classmethod + def validate_lng(cls, value: float) -> float: + if not -180 <= value <= 180: + raise ValueError("lng must be between -180 and 180") + return value + + +class LinkRequest(BaseModel): + model_config = ConfigDict(populate_by_name=True) + + from_: GeoAntenna = Field(alias="from") + to: GeoAntenna + frequency_mhz: float = Field(gt=0) + fresnel_clearance: float = Field(default=0.6, ge=0, le=1) + tx_power_dbm: float = 20.0 + rx_sensitivity_dbm: float = -137.0 + sample_step_m: float = Field(default=30.0, gt=0) + + +class LinkBudget(BaseModel): + tx_power_dbm: float + rx_sensitivity_dbm: float + margin_db: float + status: Literal["good", "marginal", "poor"] + + +class Obstacle(BaseModel): + lat: float + lng: float + distance_from_start_m: float + terrain_height_m: float + canopy_height_m: float + building_height_m: float + total_height_m: float + fresnel_radius_m: float + clearance_m: float + + +class PathProfilePoint(BaseModel): + distance_m: float + lat: float + lng: float + terrain_m: float + canopy_m: float + building_m: float + total_m: float + los_height_m: float + fresnel_radius_m: float + clearance_m: float + + +class LinkResponse(BaseModel): + distance_m: float + los: bool + fresnel_clear: bool + free_space_loss_db: float + link_budget: LinkBudget + obstacles: list[Obstacle] + path_profile: list[PathProfilePoint] + + +class CoverageRequest(BaseModel): + origin: GeoAntenna + frequency_mhz: float = Field(gt=0) + radius_m: float = Field(gt=0) + resolution_m: float = Field(gt=0) + fresnel_clearance: float = Field(default=0.6, ge=0, le=1) + tx_power_dbm: float = 20.0 + rx_sensitivity_dbm: float = -137.0 + sample_step_m: float = Field(default=30.0, gt=0) + + +class CoverageResponse(BaseModel): + origin: GeoAntenna + coverage_geojson: dict[str, Any] + + +class ElevationProfilePoint(BaseModel): + distance_m: float + lat: float + lng: float + terrain_m: float + + +class ElevationProfileResponse(BaseModel): + distance_m: float + path_profile: list[ElevationProfilePoint] + + +class DataLayersHealth(BaseModel): + srtm: bool + canopy: bool + buildings: bool + + +class HealthResponse(BaseModel): + status: Literal["ok", "degraded"] + data_layers: DataLayersHealth diff --git a/app/routers/__init__.py b/app/routers/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/app/routers/coverage.py b/app/routers/coverage.py new file mode 100644 index 0000000..60a3038 --- /dev/null +++ b/app/routers/coverage.py @@ -0,0 +1,27 @@ +from __future__ import annotations + +from fastapi import APIRouter, Depends, HTTPException, Query + +from app.models.schemas import CoverageRequest, CoverageResponse +from app.services.coverage import CoverageService + +router = APIRouter(prefix="/v1", tags=["coverage"]) + + +def get_coverage_service() -> CoverageService: + from app.main import coverage_service + + return coverage_service + + +@router.post("/coverage", response_model=CoverageResponse) +async def calculate_coverage( + request: CoverageRequest, + service: CoverageService = Depends(get_coverage_service), +) -> CoverageResponse: + try: + geojson = await service.calculate_coverage(request) + except ValueError as exc: + raise HTTPException(status_code=400, detail=str(exc)) from exc + + return CoverageResponse(origin=request.origin, coverage_geojson=geojson) diff --git a/app/routers/elevation.py b/app/routers/elevation.py new file mode 100644 index 0000000..d830952 --- /dev/null +++ b/app/routers/elevation.py @@ -0,0 +1,36 @@ +from __future__ import annotations + +from fastapi import APIRouter, Depends, Query + +from app.models.schemas import ElevationProfilePoint, ElevationProfileResponse +from app.services.link_calculator import LinkCalculator + +router = APIRouter(prefix="/v1", tags=["elevation"]) + + +def get_link_calculator() -> LinkCalculator: + from app.main import link_calculator + + return link_calculator + + +@router.get("/elevation-profile", response_model=ElevationProfileResponse) +async def elevation_profile( + from_lat: float = Query(..., ge=-90, le=90), + from_lng: float = Query(..., ge=-180, le=180), + to_lat: float = Query(..., ge=-90, le=90), + to_lng: float = Query(..., ge=-180, le=180), + step_m: float = Query(default=30.0, gt=0), + calculator: LinkCalculator = Depends(get_link_calculator), +) -> ElevationProfileResponse: + distance_m, profile = await calculator.elevation_profile( + from_lat, + from_lng, + to_lat, + to_lng, + step_m, + ) + return ElevationProfileResponse( + distance_m=distance_m, + path_profile=[ElevationProfilePoint(**point) for point in profile], + ) diff --git a/app/routers/link.py b/app/routers/link.py new file mode 100644 index 0000000..b7bb676 --- /dev/null +++ b/app/routers/link.py @@ -0,0 +1,30 @@ +from __future__ import annotations + +from fastapi import APIRouter, Depends + +from app.models.schemas import LinkRequest, LinkResponse +from app.services.link_calculator import LinkCalculator + +router = APIRouter(prefix="/v1", tags=["link"]) + + +def get_link_calculator() -> LinkCalculator: + from app.main import link_calculator + + return link_calculator + + +@router.post("/link", response_model=LinkResponse) +async def calculate_link( + request: LinkRequest, + calculator: LinkCalculator = Depends(get_link_calculator), +) -> LinkResponse: + return await calculator.calculate_link( + from_pt=request.from_, + to_pt=request.to, + frequency_mhz=request.frequency_mhz, + fresnel_clearance=request.fresnel_clearance, + tx_power_dbm=request.tx_power_dbm, + rx_sensitivity_dbm=request.rx_sensitivity_dbm, + sample_step_m=request.sample_step_m, + ) diff --git a/app/services/__init__.py b/app/services/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/app/services/buildings.py b/app/services/buildings.py new file mode 100644 index 0000000..a2033b3 --- /dev/null +++ b/app/services/buildings.py @@ -0,0 +1,90 @@ +from __future__ import annotations + +import asyncio +from typing import Iterable + +from psycopg.rows import dict_row +from psycopg_pool import ConnectionPool + +from app.config import settings + + +class BuildingService: + def __init__(self, database_url: str | None = None) -> None: + self.database_url = database_url or settings.database_url + self._pool: ConnectionPool | None = None + self._available = False + + def connect(self) -> None: + try: + self._pool = ConnectionPool( + conninfo=self.database_url, + min_size=1, + max_size=10, + kwargs={"row_factory": dict_row}, + open=True, + ) + with self._pool.connection() as conn: + with conn.cursor() as cur: + cur.execute("SELECT to_regclass('public.buildings') IS NOT NULL AS exists") + row = cur.fetchone() + self._available = bool(row and row["exists"]) + except Exception: + self._available = False + if self._pool is not None: + self._pool.close() + self._pool = None + + def close(self) -> None: + if self._pool is not None: + self._pool.close() + self._pool = None + + def is_available(self) -> bool: + return self._available and self._pool is not None + + def _get_max_height_sync(self, lat: float, lng: float) -> float: + if not self.is_available() or self._pool is None: + return 0.0 + + with self._pool.connection() as conn: + with conn.cursor() as cur: + cur.execute( + """ + SELECT COALESCE(MAX(height_m), 0) AS height_m + FROM buildings + WHERE ST_Contains(geom, ST_SetSRID(ST_Point(%s, %s), 4326)) + """, + (lng, lat), + ) + row = cur.fetchone() + return float(row["height_m"]) if row else 0.0 + + def _get_max_heights_sync(self, coordinates: Iterable[tuple[float, float]]) -> list[float]: + coords = list(coordinates) + if not coords: + return [] + if not self.is_available() or self._pool is None: + return [0.0] * len(coords) + + values = [0.0] * len(coords) + with self._pool.connection() as conn: + with conn.cursor() as cur: + for index, (lat, lng) in enumerate(coords): + cur.execute( + """ + SELECT COALESCE(MAX(height_m), 0) AS height_m + FROM buildings + WHERE ST_Contains(geom, ST_SetSRID(ST_Point(%s, %s), 4326)) + """, + (lng, lat), + ) + row = cur.fetchone() + values[index] = float(row["height_m"]) if row else 0.0 + return values + + async def get_max_height(self, lat: float, lng: float) -> float: + return await asyncio.to_thread(self._get_max_height_sync, lat, lng) + + async def get_max_heights(self, coordinates: list[tuple[float, float]]) -> list[float]: + return await asyncio.to_thread(self._get_max_heights_sync, coordinates) diff --git a/app/services/canopy.py b/app/services/canopy.py new file mode 100644 index 0000000..515c840 --- /dev/null +++ b/app/services/canopy.py @@ -0,0 +1,16 @@ +from __future__ import annotations + +from pathlib import Path + +from app.services.raster_base import RasterService + + +class CanopyService(RasterService): + def __init__(self, data_dir: str | Path, cache_size: int = 16) -> None: + super().__init__(data_dir, cache_size=cache_size) + + def get_canopy_height(self, lat: float, lng: float) -> float: + return max(0.0, self.get_value(lat, lng, default=0.0)) + + def get_canopy_heights(self, coordinates: list[tuple[float, float]]) -> list[float]: + return [max(0.0, value) for value in self.get_values(coordinates, default=0.0)] diff --git a/app/services/coverage.py b/app/services/coverage.py new file mode 100644 index 0000000..7301c9f --- /dev/null +++ b/app/services/coverage.py @@ -0,0 +1,70 @@ +from __future__ import annotations + +import asyncio +from typing import Any + +from app.config import settings +from app.models.schemas import CoverageRequest, GeoAntenna +from app.services.geodesy import grid_in_circle, line_distance_m +from app.services.link_calculator import LinkCalculator + + +class CoverageService: + def __init__(self, link_calculator: LinkCalculator) -> None: + self.link_calculator = link_calculator + + async def calculate_coverage(self, request: CoverageRequest) -> dict[str, Any]: + grid_points = grid_in_circle( + request.origin.lat, + request.origin.lng, + request.radius_m, + request.resolution_m, + ) + + if len(grid_points) > settings.max_grid_points: + raise ValueError( + f"Grid has {len(grid_points)} points, exceeds max_grid_points={settings.max_grid_points}" + ) + + semaphore = asyncio.Semaphore(settings.coverage_concurrency) + features: list[dict[str, Any]] = [] + + async def evaluate_point(lat: float, lng: float) -> dict[str, Any]: + async with semaphore: + to_point = GeoAntenna(lat=lat, lng=lng, antenna_height_m=request.origin.antenna_height_m) + result = await self.link_calculator.calculate_link( + from_pt=request.origin, + to_pt=to_point, + frequency_mhz=request.frequency_mhz, + fresnel_clearance=request.fresnel_clearance, + tx_power_dbm=request.tx_power_dbm, + rx_sensitivity_dbm=request.rx_sensitivity_dbm, + sample_step_m=request.sample_step_m, + ) + + covered = result.los and result.fresnel_clear and result.link_budget.status != "poor" + distance = line_distance_m(request.origin.lat, request.origin.lng, lat, lng) + + return { + "type": "Feature", + "geometry": { + "type": "Point", + "coordinates": [round(lng, 6), round(lat, 6)], + }, + "properties": { + "covered": covered, + "los": result.los, + "fresnel_clear": result.fresnel_clear, + "margin_db": result.link_budget.margin_db, + "status": result.link_budget.status, + "distance_m": round(distance, 1), + }, + } + + tasks = [evaluate_point(lat, lng) for lat, lng in grid_points] + features = await asyncio.gather(*tasks) + + return { + "type": "FeatureCollection", + "features": features, + } diff --git a/app/services/fresnel.py b/app/services/fresnel.py new file mode 100644 index 0000000..b1e6bb9 --- /dev/null +++ b/app/services/fresnel.py @@ -0,0 +1,19 @@ +from __future__ import annotations + +import math + + +def fresnel_radius_m(d1_m: float, d2_m: float, d_total_m: float, freq_mhz: float) -> float: + if d_total_m <= 0 or freq_mhz <= 0: + return 0.0 + f_ghz = freq_mhz / 1000.0 + return 17.3 * math.sqrt((d1_m * d2_m) / (f_ghz * d_total_m)) + + +def clearance_m( + los_height_m: float, + total_height_m: float, + fresnel_radius_m: float, + fresnel_factor: float, +) -> float: + return los_height_m - total_height_m - fresnel_radius_m * fresnel_factor diff --git a/app/services/geodesy.py b/app/services/geodesy.py new file mode 100644 index 0000000..da71724 --- /dev/null +++ b/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 diff --git a/app/services/link_calculator.py b/app/services/link_calculator.py new file mode 100644 index 0000000..47e9da4 --- /dev/null +++ b/app/services/link_calculator.py @@ -0,0 +1,144 @@ +from __future__ import annotations + +from dataclasses import dataclass + +from app.models.schemas import ( + GeoAntenna, + LinkBudget, + LinkResponse, + Obstacle, + PathProfilePoint, +) +from app.services.buildings import BuildingService +from app.services.canopy import CanopyService +from app.services.fresnel import clearance_m, fresnel_radius_m +from app.services.geodesy import line_distance_m, los_height_at, sample_line +from app.services.path_loss import compute_link_budget, fspl_db +from app.services.terrain import TerrainService + + +@dataclass +class LinkCalculator: + terrain: TerrainService + canopy: CanopyService + buildings: BuildingService + + async def calculate_link( + self, + from_pt: GeoAntenna, + to_pt: GeoAntenna, + frequency_mhz: float, + fresnel_clearance: float, + tx_power_dbm: float, + rx_sensitivity_dbm: float, + sample_step_m: float, + ) -> LinkResponse: + total_distance = line_distance_m(from_pt.lat, from_pt.lng, to_pt.lat, to_pt.lng) + path_points = sample_line( + from_pt.lat, + from_pt.lng, + to_pt.lat, + to_pt.lng, + sample_step_m, + ) + + coordinates = [(lat, lng) for lat, lng, _ in path_points] + terrains = self.terrain.get_elevations(coordinates) + canopies = self.canopy.get_canopy_heights(coordinates) + building_heights = await self.buildings.get_max_heights(coordinates) + + from_terrain = self.terrain.get_elevation(from_pt.lat, from_pt.lng) + to_terrain = self.terrain.get_elevation(to_pt.lat, to_pt.lng) + from_antenna_asl = from_terrain + from_pt.antenna_height_m + to_antenna_asl = to_terrain + to_pt.antenna_height_m + + profile: list[PathProfilePoint] = [] + obstacles: list[Obstacle] = [] + los_clear = True + fresnel_clear = True + + for index, (lat, lng, distance) in enumerate(path_points): + terrain_m = terrains[index] + canopy_m = canopies[index] + building_m = building_heights[index] + obstacle_height = max(canopy_m, building_m) + total_m = terrain_m + obstacle_height + los_height = los_height_at(distance, from_antenna_asl, to_antenna_asl, total_distance) + d1 = distance + d2 = max(total_distance - distance, 0.0) + fresnel_r = fresnel_radius_m(d1, d2, total_distance, frequency_mhz) + point_clearance = clearance_m(los_height, total_m, fresnel_r, fresnel_clearance) + + if los_height < total_m: + los_clear = False + if point_clearance < 0: + fresnel_clear = False + obstacles.append( + Obstacle( + lat=round(lat, 6), + lng=round(lng, 6), + distance_from_start_m=round(distance, 1), + terrain_height_m=round(terrain_m, 1), + canopy_height_m=round(canopy_m, 1), + building_height_m=round(building_m, 1), + total_height_m=round(total_m, 1), + fresnel_radius_m=round(fresnel_r, 1), + clearance_m=round(point_clearance, 1), + ) + ) + + profile.append( + PathProfilePoint( + distance_m=round(distance, 1), + lat=round(lat, 6), + lng=round(lng, 6), + terrain_m=round(terrain_m, 1), + canopy_m=round(canopy_m, 1), + building_m=round(building_m, 1), + total_m=round(total_m, 1), + los_height_m=round(los_height, 1), + fresnel_radius_m=round(fresnel_r, 1), + clearance_m=round(point_clearance, 1), + ) + ) + + link_budget: LinkBudget = compute_link_budget( + tx_power_dbm, + rx_sensitivity_dbm, + total_distance, + frequency_mhz, + ) + + return LinkResponse( + distance_m=round(total_distance, 1), + los=los_clear, + fresnel_clear=fresnel_clear, + free_space_loss_db=round(fspl_db(total_distance, frequency_mhz), 1), + link_budget=link_budget, + obstacles=obstacles, + path_profile=profile, + ) + + async def elevation_profile( + self, + from_lat: float, + from_lng: float, + to_lat: float, + to_lng: float, + step_m: float, + ) -> tuple[float, list[dict[str, float]]]: + total_distance = line_distance_m(from_lat, from_lng, to_lat, to_lng) + path_points = sample_line(from_lat, from_lng, to_lat, to_lng, step_m) + coordinates = [(lat, lng) for lat, lng, _ in path_points] + terrains = self.terrain.get_elevations(coordinates) + + profile = [ + { + "distance_m": round(distance, 1), + "lat": round(lat, 6), + "lng": round(lng, 6), + "terrain_m": round(terrains[index], 1), + } + for index, (lat, lng, distance) in enumerate(path_points) + ] + return round(total_distance, 1), profile diff --git a/app/services/path_loss.py b/app/services/path_loss.py new file mode 100644 index 0000000..51703de --- /dev/null +++ b/app/services/path_loss.py @@ -0,0 +1,36 @@ +from __future__ import annotations + +import math +from typing import Literal + +from app.models.schemas import LinkBudget + + +def fspl_db(distance_m: float, freq_mhz: float) -> float: + if distance_m <= 0: + return 0.0 + return 20.0 * math.log10(distance_m / 1000.0) + 20.0 * math.log10(freq_mhz) + 27.55 + + +def compute_link_budget( + tx_power_dbm: float, + rx_sensitivity_dbm: float, + distance_m: float, + freq_mhz: float, +) -> LinkBudget: + loss = fspl_db(distance_m, freq_mhz) + margin = tx_power_dbm - rx_sensitivity_dbm - loss + status: Literal["good", "marginal", "poor"] + if margin >= 20: + status = "good" + elif margin >= 10: + status = "marginal" + else: + status = "poor" + + return LinkBudget( + tx_power_dbm=tx_power_dbm, + rx_sensitivity_dbm=rx_sensitivity_dbm, + margin_db=round(margin, 1), + status=status, + ) diff --git a/app/services/raster_base.py b/app/services/raster_base.py new file mode 100644 index 0000000..bec27d1 --- /dev/null +++ b/app/services/raster_base.py @@ -0,0 +1,121 @@ +from __future__ import annotations + +import math +from dataclasses import dataclass +from pathlib import Path +from typing import Iterable + +import rasterio +from cachetools import LRUCache +from rasterio.io import DatasetReader +from rasterio.transform import rowcol + + +@dataclass(frozen=True) +class TileInfo: + path: Path + bounds: tuple[float, float, float, float] # left, bottom, right, top + nodata: float | None + + +class RasterService: + """Base GeoTIFF reader with tile index and LRU cache.""" + + def __init__(self, data_dir: str | Path, cache_size: int = 16) -> None: + self.data_dir = Path(data_dir) + self._tiles: list[TileInfo] = [] + self._cache: LRUCache[str, DatasetReader] = LRUCache(maxsize=cache_size) + self._scan_tiles() + + def _scan_tiles(self) -> None: + if not self.data_dir.exists(): + return + + for path in sorted(self.data_dir.glob("*.tif")) + sorted(self.data_dir.glob("*.tiff")): + try: + with rasterio.open(path) as dataset: + bounds = dataset.bounds + nodata = dataset.nodata + except rasterio.errors.RasterioIOError: + continue + + self._tiles.append( + TileInfo( + path=path, + bounds=(bounds.left, bounds.bottom, bounds.right, bounds.top), + nodata=nodata, + ) + ) + + def is_available(self) -> bool: + return len(self._tiles) > 0 + + def _find_tile(self, lng: float, lat: float) -> TileInfo | None: + for tile in self._tiles: + left, bottom, right, top = tile.bounds + if left <= lng <= right and bottom <= lat <= top: + return tile + return None + + def _open_dataset(self, tile: TileInfo) -> DatasetReader: + key = str(tile.path) + if key not in self._cache: + self._cache[key] = rasterio.open(tile.path) + return self._cache[key] + + def close(self) -> None: + for dataset in self._cache.values(): + dataset.close() + self._cache.clear() + + def _normalize_value(self, value: float, nodata: float | None, default: float) -> float: + if not math.isfinite(value): + return default + if nodata is not None and value == nodata: + return default + return float(value) + + def get_value(self, lat: float, lng: float, default: float = 0.0) -> float: + tile = self._find_tile(lng, lat) + if tile is None: + return default + + dataset = self._open_dataset(tile) + row, col = rowcol(dataset.transform, lng, lat) + if row < 0 or col < 0 or row >= dataset.height or col >= dataset.width: + return default + + value = float(dataset.read(1, window=((row, row + 1), (col, col + 1)))[0, 0]) + return self._normalize_value(value, tile.nodata, default) + + def get_values(self, coordinates: Iterable[tuple[float, float]], default: float = 0.0) -> list[float]: + coords = list(coordinates) + if not coords: + return [] + + results = [default] * len(coords) + by_tile: dict[str, list[tuple[int, float, float]]] = {} + + for index, (lat, lng) in enumerate(coords): + tile = self._find_tile(lng, lat) + if tile is None: + continue + by_tile.setdefault(str(tile.path), []).append((index, lat, lng)) + + for tile_path, items in by_tile.items(): + tile = next(t for t in self._tiles if str(t.path) == tile_path) + dataset = self._open_dataset(tile) + lats = [lat for _, lat, _ in items] + lngs = [lng for _, _, lng in items] + rows, cols = rowcol(dataset.transform, lngs, lats) + + for (index, lat, lng), row, col in zip(items, rows, cols, strict=True): + if row < 0 or col < 0 or row >= dataset.height or col >= dataset.width: + continue + value = float(dataset.read(1, window=((row, row + 1), (col, col + 1)))[0, 0]) + results[index] = self._normalize_value(value, tile.nodata, default) + + return results + + def get_value_at_distance(self, lat: float, lng: float) -> float: + return self.get_value(lat, lng, default=0.0) diff --git a/app/services/terrain.py b/app/services/terrain.py new file mode 100644 index 0000000..958acf1 --- /dev/null +++ b/app/services/terrain.py @@ -0,0 +1,16 @@ +from __future__ import annotations + +from pathlib import Path + +from app.services.raster_base import RasterService + + +class TerrainService(RasterService): + def __init__(self, data_dir: str | Path, cache_size: int = 16) -> None: + super().__init__(data_dir, cache_size=cache_size) + + def get_elevation(self, lat: float, lng: float) -> float: + return self.get_value(lat, lng, default=0.0) + + def get_elevations(self, coordinates: list[tuple[float, float]]) -> list[float]: + return self.get_values(coordinates, default=0.0) diff --git a/data/canopy/.gitkeep b/data/canopy/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/data/osm/.gitkeep b/data/osm/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/data/srtm/.gitkeep b/data/srtm/.gitkeep new file mode 100644 index 0000000..e69de29 diff --git a/docker-compose.yml b/docker-compose.yml new file mode 100644 index 0000000..26b3d70 --- /dev/null +++ b/docker-compose.yml @@ -0,0 +1,36 @@ +services: + postgis: + image: postgis/postgis:16-3.4 + environment: + POSTGRES_USER: radio + POSTGRES_PASSWORD: radio + POSTGRES_DB: radio + ports: + - "5432:5432" + volumes: + - postgis_data:/var/lib/postgresql/data + - ./sql/init.sql:/docker-entrypoint-initdb.d/init.sql:ro + healthcheck: + test: ["CMD-SHELL", "pg_isready -U radio -d radio"] + interval: 5s + timeout: 5s + retries: 10 + + app: + build: . + ports: + - "8000:8000" + environment: + DATABASE_URL: postgresql://radio:radio@postgis:5432/radio + SRTM_PATH: /data/srtm + CANOPY_PATH: /data/canopy + SAMPLE_STEP_M: "30" + MAX_GRID_POINTS: "10000" + volumes: + - ./data:/data:ro + depends_on: + postgis: + condition: service_healthy + +volumes: + postgis_data: diff --git a/requirements.txt b/requirements.txt new file mode 100644 index 0000000..e897703 --- /dev/null +++ b/requirements.txt @@ -0,0 +1,15 @@ +fastapi>=0.110 +uvicorn[standard]>=0.27 +pydantic>=2.0 +pydantic-settings>=2.0 +rasterio>=1.3 +numpy>=1.26 +pyproj>=3.6 +psycopg[binary,pool]>=3.1 +shapely>=2.0 +geojson>=3.1 +cachetools>=5.3 +geopandas>=0.14 +pyrosm>=0.6 +osmnx>=1.9 +httpx>=0.27 diff --git a/scripts/download_canopy.py b/scripts/download_canopy.py new file mode 100644 index 0000000..58c9ec9 --- /dev/null +++ b/scripts/download_canopy.py @@ -0,0 +1,68 @@ +#!/usr/bin/env python3 +"""Download ETH Global Canopy Height 2020 tiles for a bounding box.""" + +from __future__ import annotations + +import argparse +import math +import sys +from pathlib import Path + +import httpx + +# Public tile base used by the ETH Global Canopy Height dataset portal. +CANOPY_BASE_URL = "https://storage.googleapis.com/earthenginepartners-hansen/GCH2020" + + +def tile_names_for_bbox(min_lat: float, min_lng: float, max_lat: float, max_lng: float) -> list[str]: + tiles: list[str] = [] + for lat in range(math.floor(min_lat), math.ceil(max_lat)): + for lng in range(math.floor(min_lng), math.ceil(max_lng)): + tiles.append(f"GCH2020_E{lng:03d}N{lat:02d}.tif") + return tiles + + +def download_tile(tile: str, output_dir: Path) -> None: + output_path = output_dir / tile + if output_path.exists(): + print(f"skip {tile}: already exists") + return + + url = f"{CANOPY_BASE_URL}/{tile}" + print(f"downloading {tile}...") + try: + with httpx.stream("GET", url, timeout=300.0, follow_redirects=True) as response: + response.raise_for_status() + with output_path.open("wb") as file: + for chunk in response.iter_bytes(): + file.write(chunk) + print(f"saved {output_path}") + except httpx.HTTPError as exc: + print(f"failed {tile}: {exc}", file=sys.stderr) + + +def main() -> int: + parser = argparse.ArgumentParser(description="Download ETH Global Canopy Height tiles") + parser.add_argument("--bbox", nargs=4, type=float, metavar=("MIN_LAT", "MIN_LNG", "MAX_LAT", "MAX_LNG")) + parser.add_argument("--tiles", nargs="*", help="Explicit tile filenames") + parser.add_argument("--output", default="data/canopy", help="Output directory") + args = parser.parse_args() + + output_dir = Path(args.output) + output_dir.mkdir(parents=True, exist_ok=True) + + if args.tiles: + tiles = args.tiles + elif args.bbox: + tiles = tile_names_for_bbox(*args.bbox) + else: + parser.error("provide --bbox or --tiles") + + for tile in tiles: + download_tile(tile, output_dir) + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/download_srtm.py b/scripts/download_srtm.py new file mode 100644 index 0000000..7b5d0f9 --- /dev/null +++ b/scripts/download_srtm.py @@ -0,0 +1,80 @@ +#!/usr/bin/env python3 +"""Download SRTM GL1 GeoTIFF tiles for a bounding box via OpenTopography.""" + +from __future__ import annotations + +import argparse +import math +import sys +from pathlib import Path + +import httpx + +OPENTOPOGRAPHY_URL = "https://portal.opentopography.org/API/globaldem" + + +def tile_names_for_bbox(min_lat: float, min_lng: float, max_lat: float, max_lng: float) -> list[str]: + tiles: list[str] = [] + for lat in range(math.floor(min_lat), math.ceil(max_lat)): + ns = "N" if lat >= 0 else "S" + for lng in range(math.floor(min_lng), math.ceil(max_lng)): + ew = "E" if lng >= 0 else "W" + tiles.append(f"{ns}{abs(lat):02d}{ew}{abs(lng):03d}") + return tiles + + +def download_tile(tile: str, output_dir: Path, api_key: str | None) -> None: + params = { + "demtype": "SRTMGL1", + "south": int(tile[1:3]) * (1 if tile[0] == "N" else -1), + "north": int(tile[1:3]) * (1 if tile[0] == "N" else -1) + 1, + "west": int(tile[4:7]) * (1 if tile[3] == "E" else -1), + "east": int(tile[4:7]) * (1 if tile[3] == "E" else -1) + 1, + "outputFormat": "GTiff", + } + if api_key: + params["API_Key"] = api_key + + output_path = output_dir / f"{tile}.tif" + if output_path.exists(): + print(f"skip {tile}: already exists") + return + + print(f"downloading {tile}...") + with httpx.stream("GET", OPENTOPOGRAPHY_URL, params=params, timeout=300.0, follow_redirects=True) as response: + response.raise_for_status() + with output_path.open("wb") as file: + for chunk in response.iter_bytes(): + file.write(chunk) + print(f"saved {output_path}") + + +def main() -> int: + parser = argparse.ArgumentParser(description="Download SRTM GL1 GeoTIFF tiles") + parser.add_argument("--bbox", nargs=4, type=float, metavar=("MIN_LAT", "MIN_LNG", "MAX_LAT", "MAX_LNG")) + parser.add_argument("--tiles", nargs="*", help="Explicit tile names, e.g. N59E030") + parser.add_argument("--output", default="data/srtm", help="Output directory") + parser.add_argument("--api-key", default=None, help="OpenTopography API key") + args = parser.parse_args() + + output_dir = Path(args.output) + output_dir.mkdir(parents=True, exist_ok=True) + + if args.tiles: + tiles = args.tiles + elif args.bbox: + tiles = tile_names_for_bbox(*args.bbox) + else: + parser.error("provide --bbox or --tiles") + + for tile in tiles: + try: + download_tile(tile, output_dir, args.api_key) + except httpx.HTTPError as exc: + print(f"failed {tile}: {exc}", file=sys.stderr) + + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/import_buildings.py b/scripts/import_buildings.py new file mode 100644 index 0000000..e415d35 --- /dev/null +++ b/scripts/import_buildings.py @@ -0,0 +1,136 @@ +#!/usr/bin/env python3 +"""Import OSM building footprints into PostGIS.""" + +from __future__ import annotations + +import argparse +import re +import sys +from pathlib import Path + +import geopandas as gpd +import osmnx as ox +import psycopg +from shapely.geometry import Polygon + + +def parse_height(tags: dict, default_height: float) -> float: + height_raw = tags.get("height") + if height_raw: + match = re.search(r"([\d.]+)", str(height_raw)) + if match: + return float(match.group(1)) + + levels_raw = tags.get("building:levels") + if levels_raw: + match = re.search(r"([\d.]+)", str(levels_raw)) + if match: + return float(match.group(1)) * 3.0 + + return default_height + + +def load_buildings_from_pbf(pbf_path: Path) -> gpd.GeoDataFrame: + try: + import pyrosm + except ImportError as exc: + raise RuntimeError("pyrosm is required for --pbf imports") from exc + + osm = pyrosm.OSM(str(pbf_path)) + buildings = osm.get_buildings() + if buildings is None or buildings.empty: + return gpd.GeoDataFrame(geometry=[], crs="EPSG:4326") + return buildings + + +def load_buildings_from_bbox( + min_lat: float, + min_lng: float, + max_lat: float, + max_lng: float, +) -> gpd.GeoDataFrame: + tags = {"building": True} + return ox.features_from_bbox(bbox=(max_lat, min_lat, min_lng, max_lng), tags=tags) + + +def normalize_polygon(geometry) -> Polygon | None: + if geometry is None or geometry.is_empty: + return None + if geometry.geom_type == "Polygon": + return geometry + if geometry.geom_type == "MultiPolygon": + return max(geometry.geoms, key=lambda geom: geom.area) + return None + + +def import_buildings( + buildings: gpd.GeoDataFrame, + database_url: str, + default_height: float, + truncate: bool, +) -> int: + if buildings.crs is None: + buildings = buildings.set_crs("EPSG:4326") + else: + buildings = buildings.to_crs("EPSG:4326") + + rows: list[tuple[int | None, float, str]] = [] + for _, feature in buildings.iterrows(): + polygon = normalize_polygon(feature.geometry) + if polygon is None: + continue + + tags = feature.to_dict() + height = parse_height(tags, default_height) + osm_id = feature.get("id") + if osm_id is not None: + try: + osm_id = int(osm_id) + except (TypeError, ValueError): + osm_id = None + rows.append((osm_id, height, polygon.wkt)) + + with psycopg.connect(database_url) as conn: + with conn.cursor() as cur: + if truncate: + cur.execute("TRUNCATE TABLE buildings") + for osm_id, height, wkt in rows: + cur.execute( + """ + INSERT INTO buildings (osm_id, height_m, geom) + VALUES (%s, %s, ST_GeomFromText(%s, 4326)) + """, + (osm_id, height, wkt), + ) + conn.commit() + + return len(rows) + + +def main() -> int: + parser = argparse.ArgumentParser(description="Import OSM buildings into PostGIS") + parser.add_argument("--pbf", type=Path, help="Path to OSM PBF file") + parser.add_argument("--bbox", nargs=4, type=float, metavar=("MIN_LAT", "MIN_LNG", "MAX_LAT", "MAX_LNG")) + parser.add_argument("--database-url", default="postgresql://radio:radio@localhost:5432/radio") + parser.add_argument("--default-height", type=float, default=10.0) + parser.add_argument("--truncate", action="store_true") + args = parser.parse_args() + + if args.pbf: + buildings = load_buildings_from_pbf(args.pbf) + elif args.bbox: + buildings = load_buildings_from_bbox(*args.bbox) + else: + parser.error("provide --pbf or --bbox") + + count = import_buildings(buildings, args.database_url, args.default_height, args.truncate) + print(f"imported {count} buildings") + return 0 + + +if __name__ == "__main__": + try: + raise SystemExit(main()) + except Exception as exc: + print(f"import failed: {exc}", file=sys.stderr) + raise SystemExit(1) from exc diff --git a/sql/init.sql b/sql/init.sql new file mode 100644 index 0000000..11d4b7c --- /dev/null +++ b/sql/init.sql @@ -0,0 +1,10 @@ +CREATE EXTENSION IF NOT EXISTS postgis; + +CREATE TABLE IF NOT EXISTS buildings ( + id BIGSERIAL PRIMARY KEY, + osm_id BIGINT, + height_m REAL NOT NULL DEFAULT 10, + geom GEOMETRY(Polygon, 4326) NOT NULL +); + +CREATE INDEX IF NOT EXISTS buildings_geom_gist ON buildings USING GIST (geom);