Add governed VMM flood hazard scenarios

backend/app/services/flood_hazard_acquisition_service.py

@@ -0,0 +1,551 @@
+from __future__ import annotations
+
+import hashlib
+import json
+import math
+import time
+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.flood_hazard import FloodHazardAcquireRequest, FloodHazardAcquisitionResult, FloodHazardProductRead
+from app.services.dataset_service import DatasetService
+
+
+@dataclass(frozen=True)
+class FloodHazardProduct:
+    key: str
+    display_name: str
+    mechanism: str
+    climate_context: str
+    probability_class: str
+    return_period_years: int
+    coverage_id: str
+    published_on: str
+    catalog_url: str
+
+
+class FloodHazardAcquisitionService:
+    PROVIDER = "vmm_flood_hazard"
+    SOURCE_CRS = "EPSG:31370"
+    NATIVE_RESOLUTION_M = 2.0
+    SOURCE_VALUE_UNIT = "cm"
+    NORMALIZED_VALUE_UNIT = "m"
+    NODATA = -9999.0
+    SOURCE_VERSION = "VMM OGRK flood hazard maps"
+    ATTRIBUTION = "Bron: VMM"
+    LICENSE_NOTE = "Publieke toegang; gebruik en bronvermelding volgens de metadata van VMM/GDI-Vlaanderen."
+    WCS_TILE_SIDE_M = 10_000.0
+    WCS_REQUEST_INTERVAL_SECONDS = 1.0
+    WCS_RETRY_DELAY_SECONDS = 3.0
+    WCS_TRANSIENT_STATUS_CODES = frozenset({400, 429, 502, 503, 504})
+    SERVICE_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/publieke-inspire-coverage-service-van-ogrk"
+    LIMITATION = (
+        "Gemodelleerde maximale overstromingsdiepte voor een vast kans- en klimaatscenario. "
+        "Dit is geen actuele waterstand, geen bathymetrie en geen permanente diepte of inhoud van een waterlichaam."
+    )
+
+    @staticmethod
+    def _products() -> dict[str, FloodHazardProduct]:
+        products: list[FloodHazardProduct] = []
+        probability = {
+            10: ("grote kans", "grote-kans"),
+            100: ("middelgrote kans", "middelgrote-kans"),
+            1000: ("kleine kans", "kleine-kans"),
+        }
+        for mechanism, code in (("pluviaal", "PLU"), ("fluviaal", "FLU")):
+            for climate_key, climate_code, climate_label, published_on in (
+                ("current", "noCC", "huidig klimaat", "2021-08-31"),
+                ("future_2050", "hCC", "klimaatprojectie 2050", "2021-08-31" if mechanism == "fluviaal" else "2019-12-22"),
+            ):
+                for period, (probability_label, probability_slug) in probability.items():
+                    climate_slug = (
+                        "huidig-klimaat"
+                        if climate_key == "current"
+                        else "toekomstig-klimaat-met-klimaatprojectie-2050"
+                    )
+                    catalog_url = (
+                        "https://www.vlaanderen.be/datavindplaats/catalogus/"
+                        f"overstromingsgevaarkaart-waterdiepte-{mechanism}-{climate_slug}-{probability_slug}"
+                    )
+                    products.append(
+                        FloodHazardProduct(
+                            key=f"{mechanism}_{climate_key}_t{period}",
+                            display_name=(
+                                f"{mechanism.capitalize()} - {climate_label} - {probability_label} (T{period})"
+                            ),
+                            mechanism=mechanism,
+                            climate_context=climate_label,
+                            probability_class=probability_label,
+                            return_period_years=period,
+                            coverage_id=(
+                                f"Overstromingsgevaarkaarten-{code.replace('PLU', 'PLUVIAAL').replace('FLU', 'FLUVIAAL')}:"
+                                f"waterdiepte_{code}_{climate_code}_T{period}"
+                            ),
+                            published_on=published_on,
+                            catalog_url=catalog_url,
+                        )
+                    )
+        return {product.key: product for product in products}
+
+    @staticmethod
+    def list_products() -> list[dict[str, Any]]:
+        return [
+            FloodHazardProductRead(
+                key=product.key,
+                display_name=product.display_name,
+                mechanism=product.mechanism,
+                climate_context=product.climate_context,
+                probability_class=product.probability_class,
+                return_period_years=product.return_period_years,
+                coverage_id=product.coverage_id,
+                native_resolution_m=FloodHazardAcquisitionService.NATIVE_RESOLUTION_M,
+                source_crs=FloodHazardAcquisitionService.SOURCE_CRS,
+                source_value_unit=FloodHazardAcquisitionService.SOURCE_VALUE_UNIT,
+                normalized_value_unit=FloodHazardAcquisitionService.NORMALIZED_VALUE_UNIT,
+                published_on=product.published_on,
+                catalog_url=product.catalog_url,
+                attribution=FloodHazardAcquisitionService.ATTRIBUTION,
+                limitation_message=FloodHazardAcquisitionService.LIMITATION,
+            ).model_dump()
+            for product in FloodHazardAcquisitionService._products().values()
+        ]
+
+    @staticmethod
+    def _product(product_key: str) -> FloodHazardProduct:
+        product = FloodHazardAcquisitionService._products().get(product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="FLOOD_HAZARD_PRODUCT_NOT_SUPPORTED",
+                message="Select a governed VMM fluvial or pluvial flood-depth scenario",
+                details={"product_key": product_key},
+                status_code=422,
+            )
+        return product
+
+    @staticmethod
+    def _prepared_request(payload: FloodHazardAcquireRequest, settings: Settings) -> dict[str, Any]:
+        if not settings.flood_hazard_enabled:
+            raise AppError(code="FLOOD_HAZARD_NOT_CONFIGURED", message="VMM flood-hazard acquisition is disabled", status_code=503)
+        product = FloodHazardAcquisitionService._product(payload.product_key)
+        resolution_m = float(payload.resolution_m or settings.flood_hazard_resolution_m)
+        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="Flood-hazard selection must be a finite non-empty rectangle", status_code=400)
+        transformer = Transformer.from_crs("EPSG:4326", FloodHazardAcquisitionService.SOURCE_CRS, always_xy=True)
+        metric_bounds = transformer.transform_bounds(*values, densify_pts=21)
+        width_m = float(metric_bounds[2] - metric_bounds[0])
+        height_m = float(metric_bounds[3] - metric_bounds[1])
+        if width_m < settings.flood_hazard_min_side_m or height_m < settings.flood_hazard_min_side_m:
+            raise AppError(
+                code="FLOOD_HAZARD_SELECTION_TOO_SMALL",
+                message=f"Select an area of at least {settings.flood_hazard_min_side_m:g} by {settings.flood_hazard_min_side_m:g} metres",
+                status_code=422,
+            )
+        if width_m > settings.flood_hazard_max_side_m or height_m > settings.flood_hazard_max_side_m:
+            raise AppError(
+                code="FLOOD_HAZARD_SELECTION_TOO_LARGE",
+                message=f"Select an area no larger than {settings.flood_hazard_max_side_m:g} by {settings.flood_hazard_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.flood_hazard_max_pixels:
+            raise AppError(
+                code="FLOOD_HAZARD_SELECTION_TOO_LARGE",
+                message="Flood-hazard selection exceeds the configured raster cell limit",
+                details={"pixel_count": width * height, "max_pixels": settings.flood_hazard_max_pixels},
+                status_code=422,
+            )
+        bbox_4326 = [float(value) for value in values]
+        bbox_31370 = [float(value) for value in metric_bounds]
+        identity = {
+            "provider": FloodHazardAcquisitionService.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(identity, sort_keys=True).encode()).hexdigest()
+        return {
+            **identity,
+            "product": product,
+            "request_hash": request_hash,
+            "bbox_epsg4326": bbox_4326,
+            "bbox_epsg31370": bbox_31370,
+            "width": width,
+            "height": height,
+        }
+
+    @staticmethod
+    def _wcs_request_url(settings: Settings, product: FloodHazardProduct, bounds: tuple[float, float, float, float], resolution_m: float) -> str:
+        crs = "urn:ogc:def:crs:EPSG::31370"
+        query = [
+            ("SERVICE", "WCS"),
+            ("VERSION", "1.1.0"),
+            ("REQUEST", "GetCoverage"),
+            ("IDENTIFIER", product.coverage_id),
+            ("BOUNDINGBOX", f"{bounds[0]:.3f},{bounds[1]:.3f},{bounds[2]:.3f},{bounds[3]:.3f},{crs}"),
+            ("FORMAT", "image/tiff"),
+            ("GRIDBASECRS", crs),
+            ("GRIDCS", "urn:ogc:def:cs:OGC:0.0:Grid2dSquareCS"),
+            ("GRIDTYPE", "urn:ogc:def:method:WCS:1.1:2dSimpleGrid"),
+            ("GRIDORIGIN", f"{bounds[0]:.3f},{bounds[3]:.3f}"),
+            ("GRIDOFFSETS", f"{resolution_m:g},-{resolution_m:g}"),
+        ]
+        return f"{settings.flood_hazard_wcs_url}?{urlencode(query)}"
+
+    @staticmethod
+    def _tile_bounds(prepared: dict[str, Any]) -> list[tuple[float, float, float, float]]:
+        min_x, min_y, max_x, max_y = prepared["bbox_epsg31370"]
+        tiles: list[tuple[float, float, float, float]] = []
+        y = min_y
+        while y < max_y:
+            tile_max_y = min(y + FloodHazardAcquisitionService.WCS_TILE_SIDE_M, max_y)
+            x = min_x
+            while x < max_x:
+                tile_max_x = min(x + FloodHazardAcquisitionService.WCS_TILE_SIDE_M, max_x)
+                tiles.append((x, y, tile_max_x, tile_max_y))
+                x = tile_max_x
+            y = tile_max_y
+        return tiles
+
+    @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="FLOOD_HAZARD_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+        return intersection
+
+    @staticmethod
+    def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]:
+        request = Request(request_url, headers={"Accept": "*/*", "User-Agent": "GeoIntel/0.1 bounded-vmm-flood-hazard-acquisition"})
+        max_bytes = settings.flood_hazard_max_response_mb * 1024 * 1024
+        try:
+            with (opener or urlopen)(request, timeout=settings.flood_hazard_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="FLOOD_HAZARD_RESPONSE_TOO_LARGE", message="Official VMM response exceeds the configured size limit", status_code=502)
+                content = response.read(max_bytes + 1)
+        except AppError:
+            raise
+        except HTTPError as exc:
+            preview = exc.read(300).decode("utf-8", errors="replace")
+            raise AppError(
+                code="FLOOD_HAZARD_PROVIDER_UNAVAILABLE",
+                message="The official VMM WCS could not complete the bounded request",
+                details={"reason": str(exc), "provider_status_code": int(exc.code), "response_preview": preview},
+                status_code=502,
+            ) from exc
+        except (URLError, TimeoutError, OSError) as exc:
+            raise AppError(
+                code="FLOOD_HAZARD_PROVIDER_UNAVAILABLE",
+                message="The official VMM WCS could not complete the bounded request",
+                details={"reason": str(exc)},
+                status_code=502,
+            ) from exc
+        if len(content) > max_bytes:
+            raise AppError(code="FLOOD_HAZARD_RESPONSE_TOO_LARGE", message="Official VMM 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():
+            raise AppError(
+                code="FLOOD_HAZARD_PROVIDER_INVALID_RESPONSE",
+                message="The official VMM service did not return a GeoTIFF coverage",
+                details={"content_type": content_type, "response_preview": content[:300].decode("utf-8", errors="replace")},
+                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.walk():
+            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="FLOOD_HAZARD_PROVIDER_INVALID_RESPONSE",
+            message="The official VMM multipart response contains no valid GeoTIFF coverage",
+            status_code=502,
+        )
+
+    @staticmethod
+    def _mosaic_geotiffs(coverages: list[bytes], expected_resolution_m: float) -> bytes:
+        if len(coverages) == 1:
+            return coverages[0]
+        try:
+            from rasterio.io import MemoryFile
+            from rasterio.merge import merge
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio is required to assemble VMM flood-hazard tiles", status_code=503) from exc
+        memories = [MemoryFile(content) for content in coverages]
+        sources = []
+        try:
+            sources = [memory.open() for memory in memories]
+            for source in sources:
+                if source.crs is None or source.crs.to_epsg() != 31370 or source.count != 1:
+                    raise AppError(code="FLOOD_HAZARD_TILE_MISMATCH", message="VMM coverage tiles do not share the governed CRS and band layout", status_code=502)
+                if not all(math.isclose(abs(float(value)), expected_resolution_m, rel_tol=0.02, abs_tol=0.05) for value in source.res):
+                    raise AppError(code="FLOOD_HAZARD_TILE_MISMATCH", message="VMM coverage tile resolution differs from the governed request", status_code=502)
+            mosaic, transform = merge(sources, res=(expected_resolution_m, expected_resolution_m), nodata=0.0, dtype="float32")
+            profile = sources[0].profile.copy()
+            profile.pop("blockxsize", None)
+            profile.pop("blockysize", None)
+            profile.update(driver="GTiff", width=mosaic.shape[2], height=mosaic.shape[1], count=1, dtype="float32", crs=FloodHazardAcquisitionService.SOURCE_CRS, transform=transform, nodata=0.0, compress="deflate", predictor=3)
+            with MemoryFile() as output_memory:
+                with output_memory.open(**profile) as output:
+                    output.write(mosaic)
+                return output_memory.read()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(code="FLOOD_HAZARD_TILE_MOSAIC_FAILED", message="VMM flood-hazard tiles could not be assembled", details={"reason": str(exc)}, status_code=502) from exc
+        finally:
+            for source in sources:
+                source.close()
+            for memory in memories:
+                memory.close()
+
+    @staticmethod
+    def _fetch_coverage(prepared: dict[str, Any], settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, dict[str, Any]]:
+        product: FloodHazardProduct = prepared["product"]
+        request_urls = [
+            FloodHazardAcquisitionService._wcs_request_url(settings, product, bounds, prepared["resolution_m"])
+            for bounds in FloodHazardAcquisitionService._tile_bounds(prepared)
+        ]
+        raw_hash = hashlib.sha256()
+        coverage_hash = hashlib.sha256()
+        content_types: list[str] = []
+        coverages: list[bytes] = []
+        for index, request_url in enumerate(request_urls):
+            if index > 0 and opener is None:
+                time.sleep(FloodHazardAcquisitionService.WCS_REQUEST_INTERVAL_SECONDS)
+            try:
+                raw_content, content_type = FloodHazardAcquisitionService._fetch(request_url, settings, opener)
+            except AppError as exc:
+                provider_status = (exc.details or {}).get("provider_status_code")
+                if opener is not None or provider_status not in FloodHazardAcquisitionService.WCS_TRANSIENT_STATUS_CODES:
+                    raise
+                time.sleep(FloodHazardAcquisitionService.WCS_RETRY_DELAY_SECONDS)
+                raw_content, content_type = FloodHazardAcquisitionService._fetch(request_url, settings, opener)
+            coverage = FloodHazardAcquisitionService._extract_geotiff(raw_content, content_type)
+            raw_hash.update(len(raw_content).to_bytes(8, "big")); raw_hash.update(raw_content)
+            coverage_hash.update(len(coverage).to_bytes(8, "big")); coverage_hash.update(coverage)
+            content_types.append(content_type)
+            coverages.append(coverage)
+        return FloodHazardAcquisitionService._mosaic_geotiffs(coverages, prepared["resolution_m"]), {
+            "tile_count": len(request_urls),
+            "request_urls": request_urls,
+            "response_content_types": content_types,
+            "response_sha256": raw_hash.hexdigest(),
+            "coverage_sha256": coverage_hash.hexdigest(),
+        }
+
+    @staticmethod
+    def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]:
+        try:
+            import numpy as np
+            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 flood-hazard 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="FLOOD_HAZARD_INVALID_CRS", message="VMM flood-hazard coverage must use EPSG:31370", status_code=502)
+                if source.count != 1:
+                    raise AppError(code="FLOOD_HAZARD_INVALID_BANDS", message="VMM flood-hazard coverage must contain one depth 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="FLOOD_HAZARD_INVALID_RESOLUTION", message="VMM coverage resolution differs from the governed request", status_code=502)
+                transformer = Transformer.from_crs("EPSG:4326", FloodHazardAcquisitionService.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])
+                source_values = np.ma.asarray(clipped[0], dtype="float32")
+                raw_cm = np.asarray(source_values.filled(0.0), dtype="float32")
+                positive = (~np.ma.getmaskarray(source_values)) & np.isfinite(raw_cm) & (raw_cm > 0.0)
+                normalized_m = np.full(raw_cm.shape, FloodHazardAcquisitionService.NODATA, dtype="float32")
+                normalized_m[positive] = raw_cm[positive] / 100.0
+                valid_values = normalized_m[positive].astype("float64")
+                profile = source.profile.copy()
+                profile.pop("blockxsize", None); profile.pop("blockysize", None)
+                profile.update(driver="GTiff", width=normalized_m.shape[1], height=normalized_m.shape[0], count=1, dtype="float32", crs=FloodHazardAcquisitionService.SOURCE_CRS, transform=transform, nodata=FloodHazardAcquisitionService.NODATA, compress="deflate", predictor=3)
+                with MemoryFile() as output_memory:
+                    with output_memory.open(**profile) as output:
+                        output.write(normalized_m, 1)
+                    normalized_content = output_memory.read()
+                return normalized_content, {
+                    "width": int(normalized_m.shape[1]),
+                    "height": int(normalized_m.shape[0]),
+                    "inundated_pixel_count": int(positive.sum()),
+                    "nodata_value": FloodHazardAcquisitionService.NODATA,
+                    "resolution_m": resolution,
+                    "minimum_depth_m": float(valid_values.min()) if valid_values.size else None,
+                    "maximum_depth_m": float(valid_values.max()) if valid_values.size else None,
+                    "source_value_unit": FloodHazardAcquisitionService.SOURCE_VALUE_UNIT,
+                    "normalized_value_unit": FloodHazardAcquisitionService.NORMALIZED_VALUE_UNIT,
+                }
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(code="FLOOD_HAZARD_RASTER_INVALID", message="The official VMM response is not a valid georeferenced flood-depth raster", details={"reason": str(exc)}, status_code=502) from exc
+
+    @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 == FloodHazardAcquisitionService.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 acquire(db, project_id: UUID, payload: FloodHazardAcquireRequest, *, settings: Settings | None = None, opener: Callable[..., Any] | None = None) -> dict[str, Any]:
+        resolved_settings = settings or get_settings()
+        prepared = FloodHazardAcquisitionService._prepared_request(payload, resolved_settings)
+        product: FloodHazardProduct = prepared["product"]
+        scope_geometry = FloodHazardAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
+        resolution_token = f"{prepared['resolution_m']:g}".replace(".", "p")
+        filename = f"vmm_flood_depth_{product.key}_{resolution_token}m_{prepared['request_hash'][:12]}.tif"
+        if not payload.force_refresh:
+            cached = FloodHazardAcquisitionService._cached_dataset(db, project_id, filename)
+            if cached is not None:
+                metadata = cached.source_metadata or {}
+                raster = cached.metadata_json or {}
+                return FloodHazardAcquisitionResult(
+                    output_dataset_id=cached.id,
+                    reused=True,
+                    provider=FloodHazardAcquisitionService.PROVIDER,
+                    product_key=product.key,
+                    display_name=product.display_name,
+                    mechanism=product.mechanism,
+                    climate_context=product.climate_context,
+                    probability_class=product.probability_class,
+                    return_period_years=product.return_period_years,
+                    coverage_id=product.coverage_id,
+                    resolution_m=float(metadata.get("analysis_resolution_m", prepared["resolution_m"])),
+                    width=int(raster.get("width", prepared["width"])),
+                    height=int(raster.get("height", prepared["height"])),
+                    inundated_pixel_count=int(metadata.get("inundated_pixel_count", 0)),
+                    bbox_epsg4326=prepared["bbox_epsg4326"],
+                    bbox_epsg31370=prepared["bbox_epsg31370"],
+                    attribution=FloodHazardAcquisitionService.ATTRIBUTION,
+                    limitation_message=FloodHazardAcquisitionService.LIMITATION,
+                ).model_dump(mode="json")
+        content, transfer = FloodHazardAcquisitionService._fetch_coverage(prepared, resolved_settings, opener)
+        normalized, validation = FloodHazardAcquisitionService._normalize_raster(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,
+            source=f"VMM OGRK WCS {product.coverage_id}",
+            source_name=FloodHazardAcquisitionService.PROVIDER,
+            source_version=FloodHazardAcquisitionService.SOURCE_VERSION,
+            content_type="image/tiff",
+            source_metadata={
+                "provider": FloodHazardAcquisitionService.PROVIDER,
+                "service": "WCS",
+                "service_version": "1.1.0",
+                "product_key": product.key,
+                "product_display_name": product.display_name,
+                "mechanism": product.mechanism,
+                "climate_context": product.climate_context,
+                "probability_class": product.probability_class,
+                "return_period_years": product.return_period_years,
+                "coverage_id": product.coverage_id,
+                "native_resolution_m": FloodHazardAcquisitionService.NATIVE_RESOLUTION_M,
+                "analysis_resolution_m": validation["resolution_m"],
+                "source_crs": FloodHazardAcquisitionService.SOURCE_CRS,
+                "source_value_unit": FloodHazardAcquisitionService.SOURCE_VALUE_UNIT,
+                "normalized_value_unit": FloodHazardAcquisitionService.NORMALIZED_VALUE_UNIT,
+                "inundated_pixel_count": validation["inundated_pixel_count"],
+                "minimum_depth_m": validation["minimum_depth_m"],
+                "maximum_depth_m": validation["maximum_depth_m"],
+                "bbox_epsg4326": prepared["bbox_epsg4326"],
+                "bbox_epsg31370": prepared["bbox_epsg31370"],
+                "published_on": product.published_on,
+                "catalog_url": product.catalog_url,
+                "service_catalog_url": FloodHazardAcquisitionService.SERVICE_CATALOG_URL,
+                "attribution": FloodHazardAcquisitionService.ATTRIBUTION,
+                "license_note": FloodHazardAcquisitionService.LICENSE_NOTE,
+                "theme": "flood_hazard",
+                "layer_name": "modelled_flood_depth",
+                "coverage_scope": "municipality" if payload.area_id else "bounded_selection",
+            },
+            provenance_metadata={
+                "acquisition": "explicit_bounded_tiled_wcs_coverage",
+                "acquired_at": acquired_at.isoformat(),
+                "request_hash": prepared["request_hash"],
+                "tile_count": transfer["tile_count"],
+                "tile_request_urls": transfer["request_urls"],
+                "response_content_types": transfer["response_content_types"],
+                "response_sha256": transfer["response_sha256"],
+                "coverage_sha256": transfer["coverage_sha256"],
+                "normalized_sha256": hashlib.sha256(normalized).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": FloodHazardAcquisitionService.LIMITATION,
+                "bathymetry_available": False,
+                "permanent_water_depth_available": False,
+                "permanent_water_volume_available": False,
+                "concurrent_flood_volume_available": False,
+            },
+        )
+        return FloodHazardAcquisitionResult(
+            output_dataset_id=dataset.id,
+            reused=False,
+            provider=FloodHazardAcquisitionService.PROVIDER,
+            product_key=product.key,
+            display_name=product.display_name,
+            mechanism=product.mechanism,
+            climate_context=product.climate_context,
+            probability_class=product.probability_class,
+            return_period_years=product.return_period_years,
+            coverage_id=product.coverage_id,
+            resolution_m=validation["resolution_m"],
+            width=validation["width"],
+            height=validation["height"],
+            inundated_pixel_count=validation["inundated_pixel_count"],
+            bbox_epsg4326=prepared["bbox_epsg4326"],
+            bbox_epsg31370=prepared["bbox_epsg31370"],
+            attribution=FloodHazardAcquisitionService.ATTRIBUTION,
+            limitation_message=FloodHazardAcquisitionService.LIMITATION,
+        ).model_dump(mode="json")

backend/app/services/flood_hazard_analysis_service.py

@@ -0,0 +1,238 @@
+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.flood_hazard import (
+    FloodHazardMetric,
+    FloodHazardSelectionRequest,
+    FloodHazardSelectionResponse,
+    FloodHazardSelectionSummary,
+)
+from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
+
+
+class FloodHazardAnalysisService:
+    UNSUPPORTED_METRICS = [
+        "bathymetry_depth_m",
+        "permanent_water_volume_m3",
+        "concurrent_flood_volume_m3",
+    ]
+    LIMITATION = (
+        "Alle waarden horen bij het gekozen VMM-overstromingsscenario. De diepte-oppervlakte-integraal telt lokale "
+        "gemodelleerde maxima op en is geen gelijktijdig opgeslagen watervolume, actuele waterstand of bathymetrie."
+    )
+
+    @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 != FloodHazardAcquisitionService.PROVIDER:
+            raise AppError(
+                code="INVALID_FLOOD_HAZARD_DATASET",
+                message="Flood-hazard analysis requires a governed VMM flood-depth raster",
+                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 VMM flood-depth raster is unavailable", status_code=404)
+        return dataset
+
+    @staticmethod
+    def _selection_geometry(db, project_id: UUID, payload: FloodHazardSelectionRequest):
+        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)
+        intersection = selection.intersection(to_shape(area.geometry))
+        if intersection.is_empty or intersection.area <= 0:
+            raise AppError(code="FLOOD_HAZARD_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+        return intersection
+
+    @staticmethod
+    def analyze(
+        db,
+        project_id: UUID,
+        dataset_id: UUID,
+        payload: FloodHazardSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        dataset = FloodHazardAnalysisService._load_dataset(db, project_id, dataset_id)
+        selection_4326 = FloodHazardAnalysisService._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 flood-hazard analysis", status_code=503) from exc
+
+        source_metadata = dataset.source_metadata or {}
+        product_key = str(source_metadata.get("product_key") or "")
+        product = FloodHazardAcquisitionService._products().get(product_key)
+        if product is None or str(source_metadata.get("normalized_value_unit") or "") != "m":
+            raise AppError(code="INVALID_FLOOD_HAZARD_METADATA", message="VMM flood-hazard 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="VMM flood-depth 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)
+                analysis_geometry = selection_metric.intersection(box(*source.bounds))
+                if analysis_geometry.is_empty or analysis_geometry.area <= 0:
+                    raise AppError(code="FLOOD_HAZARD_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted flood-depth 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.flood_hazard_max_pixels:
+                    raise AppError(
+                        code="FLOOD_HAZARD_SELECTION_TOO_LARGE",
+                        message="Flood-hazard analysis exceeds the configured raster cell limit",
+                        details={"pixel_count": expected_cells, "max_pixels": resolved_settings.flood_hazard_max_pixels},
+                        status_code=422,
+                    )
+                clipped, clipped_transform = mask(source, [mapping(analysis_geometry)], crop=True, filled=False, indexes=[1])
+                depth = np.ma.asarray(clipped[0], dtype="float64")
+                raw = depth.filled(np.nan)
+                selected_cells = geometry_mask([mapping(analysis_geometry)], out_shape=depth.shape, transform=clipped_transform, invert=True)
+                nodata = source.nodata
+                valid = selected_cells & ~np.ma.getmaskarray(depth) & np.isfinite(raw) & (raw > 0.0)
+                if nodata is not None:
+                    valid &= raw != float(nodata)
+                values = raw[valid]
+                selected_cell_count = int(selected_cells.sum())
+                inundated_cell_count = int(values.size)
+                resolution_x = abs(float(source.res[0]))
+                resolution_y = abs(float(source.res[1]))
+                cell_area_m2 = resolution_x * resolution_y
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="FLOOD_HAZARD_ANALYSIS_FAILED",
+                message="The persisted VMM flood-depth 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) -> FloodHazardMetric:
+            return FloodHazardMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        inundated_area_ha = inundated_cell_count * cell_area_m2 / 10_000.0
+        metrics = [
+            metric("modelled_inundated_area_ha", "Gemodelleerd overstroomd oppervlak", inundated_area_ha, "ha", "positive_depth_cells_times_cell_area"),
+            metric(
+                "modelled_inundated_share_pct",
+                "Aandeel selectie met gemodelleerde diepte",
+                inundated_cell_count / max(1, selected_cell_count) * 100.0,
+                "%",
+                "positive_depth_cells_divided_by_selected_cells",
+            ),
+        ]
+        if inundated_cell_count:
+            metrics.extend(
+                [
+                    metric("modelled_depth_mean_m", "Gemiddelde gemodelleerde maximumdiepte", values.mean(), "m", "mean_positive_depth_cells"),
+                    metric("modelled_depth_p90_m", "90e percentiel gemodelleerde maximumdiepte", np.percentile(values, 90), "m", "percentile_90_positive_depth_cells"),
+                    metric("modelled_depth_max_m", "Hoogste gemodelleerde maximumdiepte", values.max(), "m", "maximum_positive_depth_cells"),
+                    metric(
+                        "modelled_max_depth_area_integral_m3",
+                        "Diepte-oppervlakte-integraal (geen gelijktijdig volume)",
+                        values.sum() * cell_area_m2,
+                        "m3",
+                        "sum_local_max_depth_times_cell_area",
+                    ),
+                ]
+            )
+        primary = metrics[0]
+        response = FloodHazardSelectionResponse(
+            dataset_id=dataset.id,
+            product_key=product.key,
+            mechanism=product.mechanism,
+            climate_context=product.climate_context,
+            probability_class=product.probability_class,
+            return_period_years=product.return_period_years,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            selected_cell_count=selected_cell_count,
+            inundated_cell_count=inundated_cell_count,
+            inundated_fraction=round(inundated_cell_count / max(1, selected_cell_count), 6),
+            resolution_m=round(max(resolution_x, resolution_y), 4),
+            summary=FloodHazardSelectionSummary(
+                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=FloodHazardAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=FloodHazardAnalysisService.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 = FloodHazardAnalysisService._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 flood-hazard 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) & (values > 0.0)
+                normalized = np.clip(values / 2.0, 0.0, 1.0)
+                normalized = np.where(valid, normalized, 0.0)
+                stops = np.asarray([0.0, 0.15, 0.35, 0.65, 1.0])
+                colors = np.asarray(
+                    [[190, 228, 255], [105, 184, 235], [42, 132, 201], [19, 83, 154], [8, 36, 92]],
+                    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, np.clip(150 + normalized * 90, 0, 235), 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="FLOOD_HAZARD_PREVIEW_FAILED",
+                message="The persisted VMM flood-depth raster could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc

backend/app/services/geo_assistant_service.py

@@ -21,6 +21,9 @@ from app.schemas.assistant import (
     AssistantStatus,
     AssistantTemporalSeries,
 )
+from app.schemas.flood_hazard import FloodHazardSelectionRequest
+from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
+from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService
 from app.services.vector_feature_service import VectorFeatureService
 
 
@@ -252,16 +255,22 @@ class GeoAssistantService:
             db.query(Dataset)
             .filter(Dataset.project_id == project_id)
             .filter(Dataset.status == "ready")
-            .filter(Dataset.dataset_type.in_(["vector", "geojson"]))
             .all()
         )
+        vector_datasets = [dataset for dataset in datasets if dataset.dataset_type in {"vector", "geojson"}]
+        flood_hazard_datasets = [
+            dataset
+            for dataset in datasets
+            if dataset.dataset_type == "raster" and dataset.source_name == FloodHazardAcquisitionService.PROVIDER
+            and (area is None or dataset.area_id is None or dataset.area_id == area.id)
+        ]
         warnings: list[str] = []
         context_metrics: list[AssistantContextMetric] = []
         source_dataset_ids: list[UUID] = []
         current_context: list[dict[str, Any]] = []
 
         if bbox is not None:
-            for dataset in self._current_datasets(datasets):
+            for dataset in self._current_datasets(vector_datasets):
                 kwargs: dict[str, Any] = {"dataset": dataset, "bbox": bbox}
                 if area is not None:
                     kwargs["selection_geometry"] = area.geometry
@@ -314,10 +323,60 @@ class GeoAssistantService:
                     }
                 )
 
+            for dataset in sorted(
+                flood_hazard_datasets,
+                key=lambda item: str((item.source_metadata or {}).get("product_key") or item.name),
+            ):
+                try:
+                    result = FloodHazardAnalysisService.analyze(
+                        db,
+                        project_id,
+                        dataset.id,
+                        FloodHazardSelectionRequest(bbox=bbox, area_id=area.id if area is not None else None),
+                        settings=self.settings,
+                    )
+                except AppError as exc:
+                    warnings.append(f"{dataset.name}: {exc.message}")
+                    continue
+                metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
+                scenario_label = str(metadata.get("product_display_name") or result["product_key"])
+                serialized_metrics: list[dict[str, Any]] = []
+                for metric in result["summary"]["metrics"]:
+                    item = AssistantContextMetric(
+                        theme="flood_hazard",
+                        label=f"{metric['metric_label']} - {scenario_label}",
+                        value=float(metric["metric_value"]),
+                        unit=str(metric["metric_unit"]),
+                        source=FloodHazardAcquisitionService.ATTRIBUTION,
+                        dataset_id=dataset.id,
+                        is_estimate=False,
+                    )
+                    context_metrics.append(item)
+                    serialized_metrics.append(item.model_dump(mode="json"))
+                    serialized_metrics[-1]["measurement_quality"] = "exacte_berekening_binnen_gemodelleerd_scenario"
+                source_dataset_ids.append(dataset.id)
+                current_context.append(
+                    {
+                        "dataset_name": dataset.name,
+                        "dataset_id": str(dataset.id),
+                        "theme": "flood_hazard",
+                        "source": FloodHazardAcquisitionService.ATTRIBUTION,
+                        "scenario": {
+                            "label": scenario_label,
+                            "mechanism": result["mechanism"],
+                            "climate_context": result["climate_context"],
+                            "probability_class": result["probability_class"],
+                            "return_period_years": result["return_period_years"],
+                        },
+                        "metrics": serialized_metrics,
+                        "warning": result["limitation_message"],
+                    }
+                )
+
         temporal_series: list[AssistantTemporalSeries] = []
         temporal_context: list[dict[str, Any]] = []
         include_history = self.history_requested(payload.question)
-        for key, observations in self._series(datasets):
+        for key, observations in self._series(vector_datasets):
             first = observations[0]
             last = observations[-1]
             source_metadata = last.source_metadata if isinstance(last.source_metadata, dict) else {}
@@ -364,7 +423,9 @@ class GeoAssistantService:
             "available_temporal_series": temporal_context,
             "rules": {
                 "water_volume_available": False,
-                "water_volume_reason": "Geen gebiedsdekkende waterdiepte of bathymetrie gekoppeld.",
+                "water_volume_reason": "Geen bathymetrie gekoppeld voor de permanente inhoud van waterlichamen.",
+                "flood_hazard_scenarios_available": bool(flood_hazard_datasets),
+                "flood_depth_area_integral_is_concurrent_volume": False,
                 "object_counts_are_supporting_metrics": True,
                 "causal_explanations_available": False,
                 "forecast_available": False,
@@ -402,6 +463,7 @@ class GeoAssistantService:
             "Neem waarden en jaren letterlijk over en bereken zelf geen gemiddelde, tempo, oorzaak of afgeleide trend. "
             "Gebruik platte tekst met korte alinea's en opsommingen, zonder Markdown-symbolen. "
             "Bereken of suggereer nooit watervolume zonder gekoppelde diepte of bathymetrie. "
+            "Noem de VMM-diepte-oppervlakte-integraal nooit een werkelijk, permanent of gelijktijdig watervolume. "
             "Als de gevraagde informatie niet in de context staat, zeg precies welke bron of meting ontbreekt. "
             "CONTEXT_JSON:\n" + json.dumps(context, ensure_ascii=False, separators=(",", ":"))
         )