Add governed DHMV terrain analysis

backend/app/services/dhmv_acquisition_service.py

@@ -0,0 +1,501 @@
+from __future__ import annotations
+
+import hashlib
+import json
+import math
+from dataclasses import dataclass
+from datetime import UTC, datetime
+from email.parser import BytesParser
+from email.policy import default
+from pathlib import Path
+from typing import Any, Callable
+from urllib.error import HTTPError, URLError
+from urllib.parse import urlencode
+from urllib.request import Request, urlopen
+from uuid import UUID
+
+from geoalchemy2.shape import to_shape
+from pyproj import Transformer
+from shapely.geometry import box, mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.config import Settings, get_settings
+from app.core.errors import AppError
+from app.models import Area, Dataset, Project
+from app.schemas.dhmv import DhmvAcquireRequest, DhmvAcquisitionResult, DhmvProductRead
+from app.services.dataset_service import DatasetService
+
+
+@dataclass(frozen=True)
+class DhmvProduct:
+    key: str
+    display_name: str
+    surface_model: str
+    coverage_id: str
+    native_resolution_m: float
+    catalog_url: str
+    limitation_message: str
+
+
+class DhmvAcquisitionService:
+    PROVIDER = "digitaal_vlaanderen_dhmv"
+    SOURCE_CRS = "EPSG:31370"
+    VERTICAL_REFERENCE = "TAW (Tweede Algemene Waterpassing)"
+    ACQUISITION_PERIOD = "2013-2015"
+    SOURCE_VERSION = "DHMV II 2014.01"
+    NODATA = -9999.0
+    ATTRIBUTION = "Bron: Digitaal Vlaanderen, Digitaal Hoogtemodel Vlaanderen II"
+    LICENSE_NOTE = "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen."
+    DTM_CATALOG_URL = (
+        "https://www.vlaanderen.be/datavindplaats/catalogus/"
+        "digitaal-hoogtemodel-vlaanderen-ii-dtm-raster-1-m"
+    )
+    DSM_CATALOG_URL = (
+        "https://www.vlaanderen.be/datavindplaats/catalogus/"
+        "digitaal-hoogtemodel-vlaanderen-ii-dsm-raster-1-m"
+    )
+
+    @staticmethod
+    def _products() -> dict[str, DhmvProduct]:
+        products = (
+            DhmvProduct(
+                key="dtm_1m",
+                display_name="DHMV II terreinmodel (DTM)",
+                surface_model="terrain",
+                coverage_id="DHMVII_DTM_1m",
+                native_resolution_m=1.0,
+                catalog_url=DhmvAcquisitionService.DTM_CATALOG_URL,
+                limitation_message=(
+                    "Maaiveldhoogte uit de opnameperiode 2013-2015. Gebouwen en andere objecten zijn verwijderd. "
+                    "Afstroming is een afgeleide interpretatie; dit product bevat geen waterdiepte."
+                ),
+            ),
+            DhmvProduct(
+                key="dsm_1m",
+                display_name="DHMV II oppervlaktemodel (DSM)",
+                surface_model="surface",
+                coverage_id="DHMVII_DSM_1m",
+                native_resolution_m=1.0,
+                catalog_url=DhmvAcquisitionService.DSM_CATALOG_URL,
+                limitation_message=(
+                    "Oppervlaktehoogte uit de opnameperiode 2013-2015, inclusief gebouwen en vegetatie. "
+                    "Dit is geen maaiveldmodel, waterdiepte of rechtstreeks gebouwhoogteproduct."
+                ),
+            ),
+        )
+        return {product.key: product for product in products}
+
+    @staticmethod
+    def list_products() -> list[dict[str, Any]]:
+        return [
+            DhmvProductRead(
+                key=product.key,
+                display_name=product.display_name,
+                surface_model=product.surface_model,
+                coverage_id=product.coverage_id,
+                native_resolution_m=product.native_resolution_m,
+                source_crs=DhmvAcquisitionService.SOURCE_CRS,
+                vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
+                acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
+                catalog_url=product.catalog_url,
+                attribution=DhmvAcquisitionService.ATTRIBUTION,
+                limitation_message=product.limitation_message,
+            ).model_dump()
+            for product in DhmvAcquisitionService._products().values()
+        ]
+
+    @staticmethod
+    def _product(product_key: str) -> DhmvProduct:
+        product = DhmvAcquisitionService._products().get(product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="DHMV_PRODUCT_NOT_SUPPORTED",
+                message="Select DTM or DSM from the governed DHMV II product registry",
+                details={"product_key": product_key},
+                status_code=422,
+            )
+        return product
+
+    @staticmethod
+    def _prepared_request(payload: DhmvAcquireRequest, settings: Settings) -> dict[str, Any]:
+        if not settings.dhmv_enabled:
+            raise AppError(code="DHMV_NOT_CONFIGURED", message="DHMV acquisition is disabled", status_code=503)
+        product = DhmvAcquisitionService._product(payload.product_key)
+        resolution_m = float(payload.resolution_m or settings.dhmv_resolution_m)
+        if resolution_m < product.native_resolution_m or resolution_m > 10.0:
+            raise AppError(
+                code="DHMV_RESOLUTION_NOT_SUPPORTED",
+                message="DHMV analysis resolution must be between the native 1 metre and 10 metres",
+                status_code=422,
+            )
+        values = (payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+        if not all(math.isfinite(value) for value in values) or payload.bbox.min_x >= payload.bbox.max_x or payload.bbox.min_y >= payload.bbox.max_y:
+            raise AppError(code="INVALID_BBOX", message="DHMV selection must be a finite non-empty rectangle", status_code=400)
+
+        transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.SOURCE_CRS, always_xy=True)
+        lambert_bounds = transformer.transform_bounds(*values, densify_pts=21)
+        width_m = float(lambert_bounds[2] - lambert_bounds[0])
+        height_m = float(lambert_bounds[3] - lambert_bounds[1])
+        if width_m < settings.dhmv_min_side_m or height_m < settings.dhmv_min_side_m:
+            raise AppError(
+                code="DHMV_SELECTION_TOO_SMALL",
+                message=f"Select an area of at least {settings.dhmv_min_side_m:g} by {settings.dhmv_min_side_m:g} metres",
+                status_code=422,
+            )
+        if width_m > settings.dhmv_max_side_m or height_m > settings.dhmv_max_side_m:
+            raise AppError(
+                code="DHMV_SELECTION_TOO_LARGE",
+                message=f"Select an area no larger than {settings.dhmv_max_side_m:g} by {settings.dhmv_max_side_m:g} metres",
+                details={"width_m": width_m, "height_m": height_m},
+                status_code=422,
+            )
+        width = max(1, math.ceil(width_m / resolution_m))
+        height = max(1, math.ceil(height_m / resolution_m))
+        if width * height > settings.dhmv_max_pixels:
+            raise AppError(
+                code="DHMV_SELECTION_TOO_LARGE",
+                message="DHMV selection exceeds the configured raster cell limit",
+                details={"pixel_count": width * height, "max_pixels": settings.dhmv_max_pixels},
+                status_code=422,
+            )
+
+        bbox_4326 = [float(value) for value in values]
+        bbox_31370 = [float(value) for value in lambert_bounds]
+        request_identity = {
+            "provider": DhmvAcquisitionService.PROVIDER,
+            "coverage_id": product.coverage_id,
+            "bbox_epsg4326": [round(value, 8) for value in bbox_4326],
+            "bbox_epsg31370": [round(value, 3) for value in bbox_31370],
+            "resolution_m": resolution_m,
+            "area_id": str(payload.area_id) if payload.area_id else None,
+        }
+        request_hash = hashlib.sha256(json.dumps(request_identity, sort_keys=True).encode()).hexdigest()
+        params = {
+            "SERVICE": "WCS",
+            "VERSION": "2.0.1",
+            "REQUEST": "GetCoverage",
+            "COVERAGEID": product.coverage_id,
+            "FORMAT": "image/tiff",
+            "SUBSET": [
+                f"x({bbox_31370[0]:.3f},{bbox_31370[2]:.3f})",
+                f"y({bbox_31370[1]:.3f},{bbox_31370[3]:.3f})",
+            ],
+            "SCALEFACTOR": f"{resolution_m / product.native_resolution_m:g}",
+        }
+        query = [
+            ("SERVICE", params["SERVICE"]),
+            ("VERSION", params["VERSION"]),
+            ("REQUEST", params["REQUEST"]),
+            ("COVERAGEID", params["COVERAGEID"]),
+            ("FORMAT", params["FORMAT"]),
+            ("SUBSET", params["SUBSET"][0]),
+            ("SUBSET", params["SUBSET"][1]),
+            ("SCALEFACTOR", params["SCALEFACTOR"]),
+        ]
+        return {
+            **request_identity,
+            "product": product,
+            "request_hash": request_hash,
+            "request_url": f"{settings.dhmv_wcs_url}?{urlencode(query)}",
+            "params": params,
+            "bbox_epsg4326": bbox_4326,
+            "bbox_epsg31370": bbox_31370,
+            "width": width,
+            "height": height,
+        }
+
+    @staticmethod
+    def _scope_geometry(db, project_id: UUID, area_id: UUID | None, bbox_epsg4326: list[float]):
+        if not db.get(Project, project_id):
+            raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
+        selection = box(*bbox_epsg4326)
+        if area_id is None:
+            return selection
+        area = db.get(Area, area_id)
+        if not area:
+            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+        if area.project_id != project_id:
+            raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400)
+        intersection = to_shape(area.geometry).intersection(selection)
+        if intersection.is_empty or intersection.area <= 0:
+            raise AppError(
+                code="DHMV_SELECTION_OUTSIDE_AREA",
+                message="The DHMV selection does not overlap the selected work area",
+                status_code=422,
+            )
+        return intersection
+
+    @staticmethod
+    def _cached_dataset(db, project_id: UUID, filename: str) -> Dataset | None:
+        candidate = (
+            db.query(Dataset)
+            .filter(
+                Dataset.project_id == project_id,
+                Dataset.name == filename,
+                Dataset.source_name == DhmvAcquisitionService.PROVIDER,
+                Dataset.status == "ready",
+            )
+            .order_by(Dataset.imported_at.desc())
+            .first()
+        )
+        if candidate and candidate.storage_path and Path(candidate.storage_path).is_file():
+            return candidate
+        return None
+
+    @staticmethod
+    def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]:
+        request = Request(request_url, headers={"User-Agent": "GeoIntel/0.1 bounded-dhmv-acquisition"})
+        max_bytes = settings.dhmv_max_response_mb * 1024 * 1024
+        try:
+            with (opener or urlopen)(request, timeout=settings.dhmv_timeout_seconds) as response:
+                content_type = str(response.headers.get("Content-Type", ""))
+                content_length = response.headers.get("Content-Length")
+                if content_length and int(content_length) > max_bytes:
+                    raise AppError(code="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502)
+                content = response.read(max_bytes + 1)
+        except AppError:
+            raise
+        except (HTTPError, URLError, TimeoutError, OSError) as exc:
+            raise AppError(
+                code="DHMV_PROVIDER_UNAVAILABLE",
+                message="The official DHMV WCS could not complete the bounded request",
+                details={"reason": str(exc)},
+                status_code=502,
+            ) from exc
+        if len(content) > max_bytes:
+            raise AppError(code="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502)
+        return content, content_type
+
+    @staticmethod
+    def _extract_geotiff(content: bytes, content_type: str) -> bytes:
+        if content.startswith((b"II*\x00", b"MM\x00*")):
+            return content
+        if "multipart" not in content_type.lower():
+            preview = content[:300].decode("utf-8", errors="replace")
+            raise AppError(
+                code="DHMV_PROVIDER_INVALID_RESPONSE",
+                message="The official DHMV service did not return a GeoTIFF coverage",
+                details={"content_type": content_type, "response_preview": preview},
+                status_code=502,
+            )
+        message = BytesParser(policy=default).parsebytes(
+            f"Content-Type: {content_type}\r\nMIME-Version: 1.0\r\n\r\n".encode() + content
+        )
+        for part in message.iter_parts():
+            payload = part.get_payload(decode=True) or b""
+            if part.get_content_type() == "image/tiff" and payload.startswith((b"II*\x00", b"MM\x00*")):
+                return payload
+        raise AppError(
+            code="DHMV_PROVIDER_INVALID_RESPONSE",
+            message="The official DHMV multipart response contains no valid GeoTIFF coverage",
+            status_code=502,
+        )
+
+    @staticmethod
+    def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]:
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.io import MemoryFile
+            from rasterio.mask import mask
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for DHMV validation", status_code=503) from exc
+
+        try:
+            with MemoryFile(content) as source_memory, source_memory.open() as source:
+                if source.crs is None or source.crs.to_epsg() != 31370:
+                    raise AppError(code="DHMV_INVALID_CRS", message="DHMV coverage must use EPSG:31370", status_code=502)
+                if source.count != 1:
+                    raise AppError(code="DHMV_INVALID_BANDS", message="DHMV coverage must contain exactly one elevation band", status_code=502)
+                resolution = max(abs(float(source.res[0])), abs(float(source.res[1])))
+                if not math.isclose(resolution, prepared["resolution_m"], rel_tol=0.02, abs_tol=0.05):
+                    raise AppError(
+                        code="DHMV_INVALID_RESOLUTION",
+                        message="DHMV coverage resolution differs from the governed request",
+                        details={"expected_m": prepared["resolution_m"], "actual_m": resolution},
+                        status_code=502,
+                    )
+                transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.SOURCE_CRS, always_xy=True)
+                scope_metric = shapely_transform(transformer.transform, scope_geometry_4326)
+                clipped, transform = mask(
+                    source,
+                    [mapping(scope_metric)],
+                    crop=True,
+                    filled=False,
+                    indexes=[1],
+                )
+                band = np.ma.asarray(clipped[0], dtype="float32")
+                nodata = float(source.nodata if source.nodata is not None else DhmvAcquisitionService.NODATA)
+                invalid = ~np.isfinite(np.asarray(band.filled(np.nan), dtype="float64"))
+                combined_mask = np.ma.getmaskarray(band) | invalid | (np.asarray(band) == nodata)
+                normalized = np.ma.array(np.asarray(band, dtype="float32"), mask=combined_mask)
+                valid_pixel_count = int(normalized.count())
+                if valid_pixel_count == 0:
+                    raise AppError(code="DHMV_NO_VALID_DATA", message="DHMV coverage contains no valid elevation cells in this selection", status_code=422)
+                profile = source.profile.copy()
+                profile.pop("blockxsize", None)
+                profile.pop("blockysize", None)
+                profile.update(
+                    driver="GTiff",
+                    width=int(normalized.shape[1]),
+                    height=int(normalized.shape[0]),
+                    count=1,
+                    dtype="float32",
+                    crs=DhmvAcquisitionService.SOURCE_CRS,
+                    transform=transform,
+                    nodata=DhmvAcquisitionService.NODATA,
+                    compress="deflate",
+                    predictor=3,
+                )
+                with MemoryFile() as output_memory:
+                    with output_memory.open(**profile) as output:
+                        output.write(normalized.filled(DhmvAcquisitionService.NODATA), 1)
+                    normalized_content = output_memory.read()
+                valid_values = normalized.compressed().astype("float64")
+                return normalized_content, {
+                    "width": int(normalized.shape[1]),
+                    "height": int(normalized.shape[0]),
+                    "valid_pixel_count": valid_pixel_count,
+                    "nodata_value": DhmvAcquisitionService.NODATA,
+                    "resolution_m": resolution,
+                    "minimum_m_taw": float(valid_values.min()),
+                    "maximum_m_taw": float(valid_values.max()),
+                }
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="DHMV_RASTER_INVALID",
+                message="The official DHMV response could not be validated as a georeferenced elevation raster",
+                details={"reason": str(exc)},
+                status_code=502,
+            ) from exc
+
+    @staticmethod
+    def acquire(
+        db,
+        project_id: UUID,
+        payload: DhmvAcquireRequest,
+        *,
+        settings: Settings | None = None,
+        opener: Callable[..., Any] | None = None,
+    ) -> dict[str, Any]:
+        resolved_settings = settings or get_settings()
+        prepared = DhmvAcquisitionService._prepared_request(payload, resolved_settings)
+        product: DhmvProduct = prepared["product"]
+        scope_geometry = DhmvAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
+        resolution_token = f"{prepared['resolution_m']:g}".replace(".", "p")
+        filename = f"dhmvii_{product.surface_model}_{resolution_token}m_{prepared['request_hash'][:12]}.tif"
+        if not payload.force_refresh:
+            cached = DhmvAcquisitionService._cached_dataset(db, project_id, filename)
+            if cached is not None:
+                source_metadata = cached.source_metadata or {}
+                raster_metadata = cached.metadata_json or {}
+                return DhmvAcquisitionResult(
+                    output_dataset_id=cached.id,
+                    reused=True,
+                    provider=DhmvAcquisitionService.PROVIDER,
+                    product_key=product.key,
+                    display_name=product.display_name,
+                    surface_model=product.surface_model,
+                    coverage_id=product.coverage_id,
+                    native_resolution_m=product.native_resolution_m,
+                    resolution_m=float(source_metadata.get("analysis_resolution_m", prepared["resolution_m"])),
+                    width=int(raster_metadata.get("width", prepared["width"])),
+                    height=int(raster_metadata.get("height", prepared["height"])),
+                    valid_pixel_count=int(source_metadata.get("valid_pixel_count", 0)),
+                    nodata_value=float(raster_metadata.get("nodata", DhmvAcquisitionService.NODATA)),
+                    bbox_epsg4326=prepared["bbox_epsg4326"],
+                    bbox_epsg31370=prepared["bbox_epsg31370"],
+                    vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
+                    acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
+                    attribution=DhmvAcquisitionService.ATTRIBUTION,
+                    limitation_message=product.limitation_message,
+                ).model_dump(mode="json")
+
+        raw_content, content_type = DhmvAcquisitionService._fetch(prepared["request_url"], resolved_settings, opener)
+        coverage_content = DhmvAcquisitionService._extract_geotiff(raw_content, content_type)
+        normalized_content, validation = DhmvAcquisitionService._normalize_raster(coverage_content, scope_geometry, prepared)
+        acquired_at = datetime.now(UTC)
+        dataset = DatasetService.import_raster_bytes(
+            db,
+            project_id=project_id,
+            area_id=payload.area_id,
+            filename=filename,
+            content=normalized_content,
+            source=f"Digitaal Vlaanderen WCS {product.coverage_id}",
+            source_name=DhmvAcquisitionService.PROVIDER,
+            temporal_series_key=f"digitaal-vlaanderen:dhmvii:{product.key}:{prepared['request_hash'][:24]}",
+            observed_at=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC),
+            valid_from=datetime(2013, 1, 1, tzinfo=UTC),
+            valid_to=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC),
+            temporal_granularity="period",
+            source_version=DhmvAcquisitionService.SOURCE_VERSION,
+            content_type="image/tiff",
+            source_metadata={
+                "provider": DhmvAcquisitionService.PROVIDER,
+                "service": "WCS",
+                "service_version": "2.0.1",
+                "product_key": product.key,
+                "product_display_name": product.display_name,
+                "surface_model": product.surface_model,
+                "coverage_id": product.coverage_id,
+                "native_resolution_m": product.native_resolution_m,
+                "analysis_resolution_m": validation["resolution_m"],
+                "source_crs": DhmvAcquisitionService.SOURCE_CRS,
+                "vertical_reference": DhmvAcquisitionService.VERTICAL_REFERENCE,
+                "vertical_unit": "m",
+                "acquisition_period": DhmvAcquisitionService.ACQUISITION_PERIOD,
+                "observation_date_precision": "period",
+                "nodata_value": validation["nodata_value"],
+                "valid_pixel_count": validation["valid_pixel_count"],
+                "minimum_m_taw": validation["minimum_m_taw"],
+                "maximum_m_taw": validation["maximum_m_taw"],
+                "bbox_epsg4326": prepared["bbox_epsg4326"],
+                "bbox_epsg31370": prepared["bbox_epsg31370"],
+                "catalog_url": product.catalog_url,
+                "attribution": DhmvAcquisitionService.ATTRIBUTION,
+                "license_note": DhmvAcquisitionService.LICENSE_NOTE,
+                "theme": "elevation",
+                "coverage_scope": "municipality" if payload.area_id else "bounded_selection",
+            },
+            provenance_metadata={
+                "acquisition": "explicit_bounded_wcs_coverage",
+                "acquired_at": acquired_at.isoformat(),
+                "request_hash": prepared["request_hash"],
+                "request_url": prepared["request_url"],
+                "response_content_type": content_type,
+                "response_sha256": hashlib.sha256(raw_content).hexdigest(),
+                "coverage_sha256": hashlib.sha256(coverage_content).hexdigest(),
+                "normalized_sha256": hashlib.sha256(normalized_content).hexdigest(),
+                "bbox_epsg4326": prepared["bbox_epsg4326"],
+                "bbox_epsg31370": prepared["bbox_epsg31370"],
+                "requested_resolution_m": prepared["resolution_m"],
+                "clipped_to_area_id": str(payload.area_id) if payload.area_id else None,
+                "validation": validation,
+                "limitation_message": product.limitation_message,
+                "water_depth_available": False,
+                "water_volume_available": False,
+            },
+        )
+        return DhmvAcquisitionResult(
+            output_dataset_id=dataset.id,
+            reused=False,
+            provider=DhmvAcquisitionService.PROVIDER,
+            product_key=product.key,
+            display_name=product.display_name,
+            surface_model=product.surface_model,
+            coverage_id=product.coverage_id,
+            native_resolution_m=product.native_resolution_m,
+            resolution_m=validation["resolution_m"],
+            width=validation["width"],
+            height=validation["height"],
+            valid_pixel_count=validation["valid_pixel_count"],
+            nodata_value=validation["nodata_value"],
+            bbox_epsg4326=prepared["bbox_epsg4326"],
+            bbox_epsg31370=prepared["bbox_epsg31370"],
+            vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
+            acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
+            attribution=DhmvAcquisitionService.ATTRIBUTION,
+            limitation_message=product.limitation_message,
+        ).model_dump(mode="json")

backend/app/services/terrain_analysis_service.py

@@ -0,0 +1,254 @@
+from __future__ import annotations
+
+import io
+import math
+from datetime import UTC, datetime
+from pathlib import Path
+from uuid import UUID
+
+from geoalchemy2.shape import to_shape
+from pyproj import Transformer
+from shapely.geometry import box, mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.config import Settings, get_settings
+from app.core.errors import AppError
+from app.models import Area, Dataset
+from app.schemas.dhmv import TerrainMetric, TerrainSelectionRequest, TerrainSelectionResponse, TerrainSelectionSummary
+from app.services.dhmv_acquisition_service import DhmvAcquisitionService
+
+
+class TerrainAnalysisService:
+    UNSUPPORTED_METRICS = ["water_depth_m", "water_volume_m3"]
+    LIMITATION = (
+        "Hoogte, reliëf en helling zijn afgeleid uit DHMV II. Afstroming vraagt bijkomende hydrologische modellering. "
+        "Waterdiepte en watervolume zijn niet beschikbaar uit DTM/DSM alleen."
+    )
+
+    @staticmethod
+    def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> Dataset:
+        dataset = db.get(Dataset, dataset_id)
+        if not dataset or dataset.project_id != project_id:
+            raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
+        if dataset.dataset_type != "raster" or dataset.source_name != DhmvAcquisitionService.PROVIDER:
+            raise AppError(
+                code="INVALID_TERRAIN_DATASET",
+                message="Terrain analysis requires a governed DHMV raster dataset",
+                status_code=400,
+            )
+        if dataset.status != "ready" or not dataset.storage_path or not Path(dataset.storage_path).is_file():
+            raise AppError(code="DATASET_FILE_MISSING", message="Persisted DHMV raster file is unavailable", status_code=404)
+        return dataset
+
+    @staticmethod
+    def _selection_geometry(db, project_id: UUID, payload: TerrainSelectionRequest):
+        selection = box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+        if payload.area_id is None:
+            return selection
+        area = db.get(Area, payload.area_id)
+        if not area:
+            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+        if area.project_id != project_id:
+            raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400)
+        selection = selection.intersection(to_shape(area.geometry))
+        if selection.is_empty or selection.area <= 0:
+            raise AppError(code="TERRAIN_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+        return selection
+
+    @staticmethod
+    def analyze(
+        db,
+        project_id: UUID,
+        dataset_id: UUID,
+        payload: TerrainSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        dataset = TerrainAnalysisService._load_dataset(db, project_id, dataset_id)
+        selection_4326 = TerrainAnalysisService._selection_geometry(db, project_id, payload)
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.features import geometry_mask
+            from rasterio.mask import mask
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for terrain analysis", status_code=503) from exc
+
+        source_metadata = dataset.source_metadata or {}
+        product_key = str(source_metadata.get("product_key") or "")
+        surface_model = str(source_metadata.get("surface_model") or "")
+        if product_key not in DhmvAcquisitionService._products() or surface_model not in {"terrain", "surface"}:
+            raise AppError(code="INVALID_TERRAIN_METADATA", message="DHMV product provenance is incomplete", status_code=409)
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                if source.crs is None:
+                    raise AppError(code="INVALID_DATASET_CRS", message="DHMV raster CRS is missing", status_code=409)
+                transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True)
+                selection_metric = shapely_transform(transformer.transform, selection_4326)
+                source_extent = box(*source.bounds)
+                analysis_geometry = selection_metric.intersection(source_extent)
+                if analysis_geometry.is_empty or analysis_geometry.area <= 0:
+                    raise AppError(
+                        code="TERRAIN_SELECTION_OUTSIDE_DATASET",
+                        message="Selection does not overlap the persisted DHMV raster",
+                        status_code=422,
+                    )
+                min_x, min_y, max_x, max_y = analysis_geometry.bounds
+                expected_cells = math.ceil((max_x - min_x) / abs(source.res[0])) * math.ceil((max_y - min_y) / abs(source.res[1]))
+                if expected_cells > resolved_settings.dhmv_max_pixels:
+                    raise AppError(
+                        code="TERRAIN_SELECTION_TOO_LARGE",
+                        message="Terrain analysis exceeds the configured raster cell limit",
+                        details={"pixel_count": expected_cells, "max_pixels": resolved_settings.dhmv_max_pixels},
+                        status_code=422,
+                    )
+                clipped, clipped_transform = mask(
+                    source,
+                    [mapping(analysis_geometry)],
+                    crop=True,
+                    filled=False,
+                    indexes=[1],
+                )
+                elevation = np.ma.asarray(clipped[0], dtype="float64")
+                raw = elevation.filled(np.nan)
+                nodata = source.nodata
+                invalid = ~np.isfinite(raw)
+                if nodata is not None:
+                    invalid |= raw == float(nodata)
+                selected_cells = geometry_mask(
+                    [mapping(analysis_geometry)],
+                    out_shape=elevation.shape,
+                    transform=clipped_transform,
+                    invert=True,
+                )
+                valid_mask = selected_cells & ~np.ma.getmaskarray(elevation) & ~invalid
+                values = raw[valid_mask]
+                if values.size == 0:
+                    raise AppError(code="TERRAIN_NO_VALID_DATA", message="No valid DHMV height cells occur in this selection", status_code=422)
+
+                resolution_x = abs(float(source.res[0]))
+                resolution_y = abs(float(source.res[1]))
+                slope_values = np.asarray([], dtype="float64")
+                if raw.shape[0] >= 2 and raw.shape[1] >= 2:
+                    surface = np.where(valid_mask, raw, np.nan)
+                    gradient_y, gradient_x = np.gradient(surface, resolution_y, resolution_x)
+                    slope = np.degrees(np.arctan(np.hypot(gradient_x, gradient_y)))
+                    slope_values = slope[np.isfinite(slope) & valid_mask]
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="TERRAIN_ANALYSIS_FAILED",
+                message="The persisted DHMV raster could not be analysed",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc
+
+        def metric(key: str, label: str, value: float, unit: str, method: str) -> TerrainMetric:
+            return TerrainMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        prefix = "terrain" if surface_model == "terrain" else "surface"
+        elevation_label = "Gemiddelde maaiveldhoogte" if surface_model == "terrain" else "Gemiddelde oppervlaktehoogte"
+        metrics = [
+            metric(f"{prefix}_elevation_mean_m", elevation_label, values.mean(), "m TAW", "mean_valid_cells"),
+            metric(f"{prefix}_elevation_min_m", "Laagste hoogte", values.min(), "m TAW", "minimum_valid_cells"),
+            metric(f"{prefix}_elevation_max_m", "Hoogste hoogte", values.max(), "m TAW", "maximum_valid_cells"),
+            metric(f"{prefix}_elevation_p10_m", "10e percentiel hoogte", np.percentile(values, 10), "m TAW", "percentile_10_valid_cells"),
+            metric(f"{prefix}_elevation_p90_m", "90e percentiel hoogte", np.percentile(values, 90), "m TAW", "percentile_90_valid_cells"),
+            metric("relief_m", "Reliëfverschil", values.max() - values.min(), "m", "maximum_minus_minimum"),
+        ]
+        if slope_values.size:
+            metrics.extend(
+                [
+                    metric("slope_mean_deg", "Gemiddelde helling", slope_values.mean(), "°", "mean_finite_gradient"),
+                    metric("slope_p90_deg", "90e percentiel helling", np.percentile(slope_values, 90), "°", "percentile_90_finite_gradient"),
+                    metric("slope_max_deg", "Steilste helling", slope_values.max(), "°", "maximum_finite_gradient"),
+                ]
+            )
+        primary = metrics[0]
+        selected_cell_count = int(selected_cells.sum())
+        response = TerrainSelectionResponse(
+            dataset_id=dataset.id,
+            product_key=product_key,
+            surface_model=surface_model,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            sample_count=int(values.size),
+            slope_sample_count=int(slope_values.size),
+            coverage_ratio=round(float(values.size / max(1, selected_cell_count)), 6),
+            resolution_m=round(max(resolution_x, resolution_y), 4),
+            vertical_reference=str(source_metadata.get("vertical_reference") or DhmvAcquisitionService.VERTICAL_REFERENCE),
+            summary=TerrainSelectionSummary(
+                metric_label=primary.metric_label,
+                metric_value=primary.metric_value,
+                metric_unit=primary.metric_unit,
+                aggregation_method=primary.aggregation_method,
+                primary_metric_key=primary.metric_key,
+                metrics=metrics,
+            ),
+            unsupported_metrics=TerrainAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=TerrainAnalysisService.LIMITATION,
+            generated_at=datetime.now(UTC).isoformat(),
+        )
+        return response.model_dump(mode="json")
+
+    @staticmethod
+    def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
+        dataset = TerrainAnalysisService._load_dataset(db, project_id, dataset_id)
+        try:
+            import numpy as np
+            import rasterio
+            from PIL import Image
+            from rasterio.enums import Resampling
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio, numpy and Pillow are required for terrain rendering", status_code=503) from exc
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                scale = min(1.0, max_dimension / max(source.width, source.height))
+                width = max(1, round(source.width * scale))
+                height = max(1, round(source.height * scale))
+                data = source.read(1, out_shape=(height, width), masked=True, resampling=Resampling.bilinear)
+                values = np.asarray(data.filled(np.nan), dtype="float64")
+                valid = np.isfinite(values) & ~np.ma.getmaskarray(data)
+                if not valid.any():
+                    raise AppError(code="TERRAIN_NO_VALID_DATA", message="DHMV raster contains no renderable cells", status_code=422)
+                low, high = np.percentile(values[valid], [2, 98])
+                if high <= low:
+                    high = low + 1.0
+                normalized = np.clip((values - low) / (high - low), 0.0, 1.0)
+                stops = np.asarray([0.0, 0.25, 0.5, 0.75, 1.0])
+                colors = np.asarray(
+                    [
+                        [30, 94, 91],
+                        [79, 139, 102],
+                        [194, 183, 105],
+                        [173, 121, 79],
+                        [105, 94, 108],
+                    ],
+                    dtype="float64",
+                )
+                rgba = np.zeros((height, width, 4), dtype="uint8")
+                for channel in range(3):
+                    rgba[:, :, channel] = np.interp(normalized, stops, colors[:, channel]).astype("uint8")
+                rgba[:, :, 3] = np.where(valid, 225, 0).astype("uint8")
+                output = io.BytesIO()
+                Image.fromarray(rgba).save(output, format="PNG", optimize=True)
+                return output.getvalue()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="TERRAIN_PREVIEW_FAILED",
+                message="The persisted DHMV raster could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc