diff --git a/API.md b/API.md index b33057c..54c46b4 100644 --- a/API.md +++ b/API.md @@ -7,9 +7,9 @@ ## Состояние Реализации -- Работает: healthcheck, OpenAPI, DEM elevation/profile, OSM buildings query, terrain LOS/Fresnel/diffraction с учётом DEM и buildings, manual link budget, antenna sector/beam helpers. +- Работает: healthcheck, OpenAPI, DEM elevation/profile, OSM buildings query, landcover/canopy path sampling, terrain LOS/Fresnel/diffraction с учётом DEM и buildings, manual link budget, antenna sector/beam helpers. - Частично работает: coverage/viewshed создают job-stub, но тяжёлый расчёт ещё не реализован. -- Пока не реализовано: landcover/canopy sampling, vegetation attenuation по реальным данным, ITM/P.1812/P.452/P.676, реальные async results для coverage/viewshed. +- Пока не реализовано: vegetation attenuation по реальным данным, ITM/P.1812/P.452/P.676, реальные async results для coverage/viewshed. ## Общие Правила @@ -309,18 +309,49 @@ curl -s -X POST http://localhost:5603/api/v1/antenna/beam \ ### `POST /api/v1/landcover/path` -Контракт есть, но raster sampling landcover/canopy пока не реализован. Сейчас endpoint возвращает `501 Not Implemented`. +Сэмплит landcover raster вдоль трассы и сворачивает точки в сегменты. +Canopy raster используется опционально: если файлов нет, endpoint всё равно +работает, но `canopy_height_m` будет `null`. + +Ожидаемые данные: + +- `LANDCOVER_PATH`: GeoTIFF/COG с кодами ESA WorldCover. +- `CANOPY_PATH`: GeoTIFF/COG с высотой кроны в метрах. + +Поддерживаемые классы ESA WorldCover: + +- `10`: `tree_cover` +- `20`: `shrubland` +- `30`: `grassland` +- `40`: `cropland` +- `50`: `built_up` +- `60`: `bare_sparse_vegetation` +- `70`: `snow_ice` +- `80`: `water` +- `90`: `herbaceous_wetland` +- `95`: `mangroves` +- `100`: `moss_lichen` ```bash -curl -i -X POST http://localhost:5603/api/v1/landcover/path \ +curl -s -X POST http://localhost:5603/api/v1/landcover/path \ -H 'Content-Type: application/json' \ -d '{ "start": {"lat": 59.935, "lon": 30.305}, "end": {"lat": 59.945, "lon": 30.325}, "samples": 256 - }' + }' | jq ``` +Ключевые поля ответа: + +- `segments`: участки с одинаковым landcover class. +- `segments[].class`: класс покрова. +- `segments[].forest_type`: пока `unknown` для `tree_cover`/`mangroves`, `null` для остальных классов. +- `segments[].canopy_height_m`: средняя высота кроны на сегменте, если canopy raster доступен. +- `vegetation_depth_m`: суммарная длина участков `tree_cover`/`mangroves` по трассе. + +Если landcover raster отсутствует или не покрывает трассу, endpoint возвращает `501 Not Implemented`. + ## Viewshed ### `POST /api/v1/viewshed` diff --git a/README.md b/README.md index 8a8c74a..e8491dc 100644 --- a/README.md +++ b/README.md @@ -68,6 +68,19 @@ curl -s -X POST http://localhost:8000/api/v1/buildings/query \ -d '{"bbox":[30.30,59.93,30.33,59.95]}' | jq '.features | length' ``` +## Landcover And Canopy + +Place ESA WorldCover GeoTIFF/COG files under `data/landcover` and optional canopy +height GeoTIFF/COG files under `data/canopy`, then restart the API: + +```bash +docker compose restart api worker +``` + +The `/api/v1/landcover/path` endpoint samples all `.tif`/`.tiff` files under +`LANDCOVER_PATH` recursively. Canopy data is optional; when it is missing, +`canopy_height_m` is returned as `null`. + For local Python development: ```bash diff --git a/api/app/core/landcover.py b/api/app/core/landcover.py new file mode 100644 index 0000000..7aa1298 --- /dev/null +++ b/api/app/core/landcover.py @@ -0,0 +1,158 @@ +from __future__ import annotations + +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + +from app.core.geo import PathPoint +from app.core.raster_sampling import RasterNotConfiguredError, sample_along + + +class LandcoverNotConfiguredError(RasterNotConfiguredError): + """Raised when landcover rasters are missing or incomplete.""" + + +ESA_WORLDCOVER_CLASSES: dict[int, str] = { + 10: "tree_cover", + 20: "shrubland", + 30: "grassland", + 40: "cropland", + 50: "built_up", + 60: "bare_sparse_vegetation", + 70: "snow_ice", + 80: "water", + 90: "herbaceous_wetland", + 95: "mangroves", + 100: "moss_lichen", +} + + +@dataclass(frozen=True) +class LandcoverPoint: + distance_m: float + class_name: str + forest_type: str | None + canopy_height_m: float | None + + +@dataclass(frozen=True) +class LandcoverPathSegment: + from_m: float + to_m: float + class_name: str + forest_type: str | None + canopy_height_m: float | None + + +@dataclass(frozen=True) +class LandcoverPath: + segments: list[LandcoverPathSegment] + vegetation_depth_m: float + + +def worldcover_class(value: float) -> str: + if np.isnan(value): + return "unknown" + return ESA_WORLDCOVER_CLASSES.get(int(round(value)), "unknown") + + +def forest_type_for_class(class_name: str) -> str | None: + if class_name in {"tree_cover", "mangroves"}: + return "unknown" + return None + + +def landcover_path( + points: list[PathPoint], + landcover_path_dir: str | Path, + canopy_path_dir: str | Path | None = None, +) -> LandcoverPath: + if len(points) < 2: + raise ValueError("points must contain at least two samples") + + try: + landcover_values = sample_along(points, landcover_path_dir, "landcover") + except RasterNotConfiguredError as exc: + raise LandcoverNotConfiguredError(str(exc)) from exc + + canopy_values = _canopy_values(points, canopy_path_dir) + classified_points = [ + LandcoverPoint( + distance_m=point.distance_m, + class_name=worldcover_class(landcover_values[index]), + forest_type=forest_type_for_class(worldcover_class(landcover_values[index])), + canopy_height_m=_canopy_height(canopy_values[index]), + ) + for index, point in enumerate(points) + ] + return _segments_from_points(classified_points) + + +def _canopy_values(points: list[PathPoint], canopy_path_dir: str | Path | None) -> np.ndarray: + if canopy_path_dir is None: + return np.full(len(points), np.nan, dtype=float) + try: + return sample_along(points, canopy_path_dir, "canopy", require_all=False) + except RasterNotConfiguredError: + return np.full(len(points), np.nan, dtype=float) + + +def _canopy_height(value: float) -> float | None: + if np.isnan(value) or value < 0: + return None + return float(value) + + +def _segments_from_points(points: list[LandcoverPoint]) -> LandcoverPath: + segments: list[LandcoverPathSegment] = [] + vegetation_depth_m = 0.0 + start = points[0] + canopy_values: list[float] = [] + + for index in range(len(points) - 1): + current = points[index] + next_point = points[index + 1] + if current.canopy_height_m is not None: + canopy_values.append(current.canopy_height_m) + + if current.class_name in {"tree_cover", "mangroves"}: + vegetation_depth_m += next_point.distance_m - current.distance_m + + same_segment = ( + next_point.class_name == start.class_name + and next_point.forest_type == start.forest_type + ) + if not same_segment: + segments.append( + LandcoverPathSegment( + from_m=start.distance_m, + to_m=next_point.distance_m, + class_name=start.class_name, + forest_type=start.forest_type, + canopy_height_m=_mean_canopy(canopy_values), + ) + ) + start = next_point + canopy_values = [] + + last = points[-1] + if last.canopy_height_m is not None: + canopy_values.append(last.canopy_height_m) + segments.append( + LandcoverPathSegment( + from_m=start.distance_m, + to_m=last.distance_m, + class_name=start.class_name, + forest_type=start.forest_type, + canopy_height_m=_mean_canopy(canopy_values), + ) + ) + + return LandcoverPath(segments=segments, vegetation_depth_m=vegetation_depth_m) + + +def _mean_canopy(values: list[float]) -> float | None: + if not values: + return None + return float(np.mean(values)) diff --git a/api/app/core/raster_sampling.py b/api/app/core/raster_sampling.py new file mode 100644 index 0000000..687242e --- /dev/null +++ b/api/app/core/raster_sampling.py @@ -0,0 +1,104 @@ +from __future__ import annotations + +from collections.abc import Sequence +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.crs import CRS +from rasterio.warp import transform + +from app.core.geo import PathPoint + + +class RasterNotConfiguredError(NotImplementedError): + """Raised when raster files are missing or do not cover requested points.""" + + +def raster_files(path: str | Path) -> list[Path]: + raster_path = Path(path) + if not raster_path.exists(): + return [] + return sorted( + file_path + for pattern in ("*.tif", "*.tiff", "*.TIF", "*.TIFF") + for file_path in raster_path.rglob(pattern) + if file_path.is_file() + ) + + +def point_in_bounds(x: float, y: float, bounds: object) -> bool: + return bounds.left <= x <= bounds.right and bounds.bottom <= y <= bounds.top + + +def to_dataset_crs(lat: float, lon: float, dst_crs: CRS | None) -> tuple[float, float]: + if dst_crs is None or dst_crs == CRS.from_epsg(4326): + return lon, lat + xs, ys = transform(CRS.from_epsg(4326), dst_crs, [lon], [lat]) + return xs[0], ys[0] + + +def sample_dataset(dataset: rasterio.io.DatasetReader, lat: float, lon: float) -> float | None: + x, y = to_dataset_crs(lat, lon, dataset.crs) + if not point_in_bounds(x, y, dataset.bounds): + return None + + value = next(dataset.sample([(x, y)], masked=True))[0] + if np.ma.is_masked(value): + return None + if dataset.nodata is not None and float(value) == float(dataset.nodata): + return None + if not np.isfinite(value): + return None + return float(value) + + +def sample_at(lat: float, lon: float, path: str | Path, label: str) -> float: + files = raster_files(path) + if not files: + raise RasterNotConfiguredError(f"{label} raster files are not found in {path}") + + for file_path in files: + with rasterio.open(file_path) as dataset: + value = sample_dataset(dataset, lat, lon) + if value is not None: + return value + + raise RasterNotConfiguredError(f"No {label} raster tile covers lat={lat}, lon={lon}") + + +def sample_along( + points: Sequence[PathPoint], + path: str | Path, + label: str, + require_all: bool = True, +) -> np.ndarray: + if not points: + return np.array([], dtype=float) + + files = raster_files(path) + if not files: + raise RasterNotConfiguredError(f"{label} raster files are not found in {path}") + + values: list[float | None] = [None] * len(points) + remaining = set(range(len(points))) + + for file_path in files: + if not remaining: + break + with rasterio.open(file_path) as dataset: + for index in list(remaining): + point = points[index] + value = sample_dataset(dataset, point.lat, point.lon) + if value is not None: + values[index] = value + remaining.remove(index) + + if remaining and require_all: + missing = ", ".join(str(index) for index in sorted(remaining)[:10]) + raise RasterNotConfiguredError( + f"No {label} raster tile covers {len(remaining)} point(s), " + f"first missing indices: {missing}" + ) + + return np.array([np.nan if value is None else float(value) for value in values], dtype=float) diff --git a/api/app/services/__pycache__/landcover.cpython-313.pyc b/api/app/services/__pycache__/landcover.cpython-313.pyc index 25dbfa5..8750063 100644 Binary files a/api/app/services/__pycache__/landcover.cpython-313.pyc and b/api/app/services/__pycache__/landcover.cpython-313.pyc differ diff --git a/api/app/services/landcover.py b/api/app/services/landcover.py index 8934e66..802655b 100644 --- a/api/app/services/landcover.py +++ b/api/app/services/landcover.py @@ -1,7 +1,35 @@ -from app.models.landcover import LandcoverPathRequest, LandcoverPathResponse +from app.config import get_settings +from app.core.geo import GeoPoint, sample_path +from app.core.landcover import LandcoverNotConfiguredError, landcover_path +from app.models.landcover import LandcoverPathRequest, LandcoverPathResponse, LandcoverSegment def path_landcover(request: LandcoverPathRequest) -> LandcoverPathResponse: - raise NotImplementedError( - "Landcover and canopy raster sampling is implemented in a later stage" + settings = get_settings() + points = sample_path( + GeoPoint(lat=request.start.lat, lon=request.start.lon), + GeoPoint(lat=request.end.lat, lon=request.end.lon), + request.samples, + ) + try: + result = landcover_path( + points, + landcover_path_dir=settings.landcover_path, + canopy_path_dir=settings.canopy_path, + ) + except LandcoverNotConfiguredError as exc: + raise NotImplementedError(str(exc)) from exc + + return LandcoverPathResponse( + segments=[ + LandcoverSegment( + from_m=segment.from_m, + to_m=segment.to_m, + class_name=segment.class_name, + forest_type=segment.forest_type, + canopy_height_m=segment.canopy_height_m, + ) + for segment in result.segments + ], + vegetation_depth_m=result.vegetation_depth_m, ) diff --git a/api/tests/test_landcover.py b/api/tests/test_landcover.py new file mode 100644 index 0000000..3c5d91d --- /dev/null +++ b/api/tests/test_landcover.py @@ -0,0 +1,62 @@ +from pathlib import Path + +import numpy as np +import rasterio +from rasterio.transform import from_origin + +from app.core.geo import PathPoint +from app.core.landcover import landcover_path, worldcover_class + + +def write_raster(path: Path, data: np.ndarray, nodata: float = -9999.0) -> None: + transform = from_origin(30.0, 60.0, 0.01, 0.01) + with rasterio.open( + path, + "w", + driver="GTiff", + height=data.shape[0], + width=data.shape[1], + count=1, + dtype=data.dtype, + crs="EPSG:4326", + transform=transform, + nodata=nodata, + ) as dataset: + dataset.write(data, 1) + + +def test_worldcover_class_mapping() -> None: + assert worldcover_class(10) == "tree_cover" + assert worldcover_class(50) == "built_up" + assert worldcover_class(255) == "unknown" + + +def test_landcover_path_segments_and_vegetation_depth(tmp_path: Path) -> None: + landcover_dir = tmp_path / "landcover" + canopy_dir = tmp_path / "canopy" + landcover_dir.mkdir() + canopy_dir.mkdir() + write_raster( + landcover_dir / "worldcover.tif", + np.array([[10, 10, 50, 50]], dtype="uint8"), + nodata=0, + ) + write_raster( + canopy_dir / "canopy.tif", + np.array([[12.0, 14.0, 0.0, 0.0]], dtype="float32"), + ) + points = [ + PathPoint(lat=59.995, lon=30.005, distance_m=0), + PathPoint(lat=59.995, lon=30.015, distance_m=10), + PathPoint(lat=59.995, lon=30.025, distance_m=20), + PathPoint(lat=59.995, lon=30.035, distance_m=30), + ] + + result = landcover_path(points, landcover_dir, canopy_dir) + + assert result.vegetation_depth_m == 20 + assert len(result.segments) == 2 + assert result.segments[0].class_name == "tree_cover" + assert result.segments[0].forest_type == "unknown" + assert result.segments[0].canopy_height_m == 13.0 + assert result.segments[1].class_name == "built_up"