height correction, new defalt

flight_map/ahn.py

@@ -0,0 +1,108 @@
+"""PDOK AHN DTM coverage in RD New, elevations in metres relative to NAP."""
+
+import hashlib
+import logging
+import math
+import os
+import tempfile
+from pathlib import Path
+from urllib.parse import urlencode
+
+import rasterio
+from pyproj import Transformer
+from rasterio.errors import RasterioError
+
+from flight_map.geo import Bounds
+from flight_map.network import request_bytes
+
+LOGGER = logging.getLogger(__name__)
+WCS_URL = "https://service.pdok.nl/rws/ahn/wcs/v1_0"
+COVERAGE = "dtm_05m"
+AHN_EXTENT = Bounds(3.20030320333955, 50.728765211722, 7.273444085294, 53.555958411015)
+SIZE = 512
+
+
+def coverage_url(bounds: Bounds) -> str:
+    transform = Transformer.from_crs("EPSG:4326", "EPSG:28992", always_xy=True)
+    west, south, east, north = transform.transform_bounds(*bounds.values, densify_pts=21)
+    bbox = ",".join(
+        str(v) for v in (math.floor(west), math.floor(south), math.ceil(east), math.ceil(north))
+    )
+    return (
+        WCS_URL
+        + "?"
+        + urlencode(
+            {
+                "service": "WCS",
+                "version": "1.0.0",
+                "request": "GetCoverage",
+                "coverage": COVERAGE,
+                "crs": "EPSG:28992",
+                "response_crs": "EPSG:28992",
+                "bbox": bbox,
+                "width": SIZE,
+                "height": SIZE,
+                "format": "GEOTIFF",
+                "interpolation": "nearest neighbor",
+            }
+        )
+    )
+
+
+def validate_raster(path: Path) -> None:
+    with rasterio.open(path) as dataset:
+        if dataset.driver != "GTiff" or dataset.crs != rasterio.crs.CRS.from_epsg(28992):
+            raise ValueError("AHN response is not an RD New GeoTIFF.")
+        if dataset.count != 1 or dataset.width != SIZE or dataset.height != SIZE:
+            raise ValueError("AHN response has unexpected raster dimensions.")
+        dataset.read(1, masked=True)
+
+
+def terrain_data(bounds: Bounds, cache: Path, *, refresh: bool = False) -> Path | None:
+    if not (
+        AHN_EXTENT.contains(bounds.west, bounds.south)
+        and AHN_EXTENT.contains(bounds.east, bounds.north)
+    ):
+        LOGGER.warning("AHN height map unavailable: map bounds extend outside Dutch AHN coverage.")
+        return None
+    url = coverage_url(bounds)
+    key = hashlib.sha256(url.encode()).hexdigest()
+    path = cache / f"ahn-{key}.tif"
+    temporary: Path | None = None
+    try:
+        if path.exists() and not refresh:
+            validate_raster(path)
+            LOGGER.info("Reading cached AHN DTM (NAP): %s", path)
+            return path
+        LOGGER.info("Downloading AHN DTM height map from PDOK for %s", bounds.values)
+        content = request_bytes(url)
+        cache.mkdir(parents=True, exist_ok=True)
+        with tempfile.NamedTemporaryFile(dir=cache, suffix=".tif", delete=False) as stream:
+            temporary = Path(stream.name)
+            stream.write(content)
+            stream.flush()
+            os.fsync(stream.fileno())
+        validate_raster(temporary)
+        with rasterio.open(temporary, "r+") as dataset:
+            dataset.update_tags(
+                source=WCS_URL,
+                coverage=COVERAGE,
+                vertical_reference="NAP",
+                units="metres",
+                license="CC0",
+                request_url=url,
+            )
+        os.replace(temporary, path)
+        LOGGER.info("Cached AHN DTM (NAP): %s", path)
+        return path
+    except (RuntimeError, OSError, ValueError, RasterioError) as exc:
+        LOGGER.warning(
+            "AHN height map unavailable: %s. Land elevation will be N/A; "
+            "reported flight altitudes are unaffected. "
+            "Use --refresh-terrain to retry a damaged cache entry.",
+            exc,
+        )
+        return None
+    finally:
+        if temporary and temporary.exists():
+            temporary.unlink()