added land

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)