Add governed SPW bathymetry analysis

backend/README.md

@@ -1806,7 +1806,8 @@ matched name with `--show-names`.
 performs a bounded official VHA ArcGIS query, exact persisted-Area clipping,
 watercourse-name normalization and ordinary Dataset/VectorFeature persistence.
 `GET /api/v1/projects/{project_id}/datasets/bathymetry/sources` reports VHA as
-operational and the audited MDK/SPW candidates as unavailable for acquisition.
+operational, MDK as probe-only and the pinned SPW raster operator as
+operational.
 
 Runtime controls are `BATHYMETRY_PROFILES_ENABLED`,
 `BATHYMETRY_PROFILES_LAYER_URL`, `BATHYMETRY_WATERCOURSE_LAYER_URL`,
@@ -1854,6 +1855,29 @@ readiness state such as TLS or endpoint failure. Runtime controls are
 `MDK_BATHYMETRY_PROBE_TIMEOUT_SECONDS` and
 `MDK_BATHYMETRY_PROBE_MAX_RESPONSE_MB`. TLS verification cannot be disabled.
 
+### SPW waterbed raster
+
+The official 2023-05-23 SPW bathymetry ZIP is integrated only through the
+bounded operator. Stage the immutable ZIP under persistent storage and run:
+
+```bash
+docker exec geointel python /app/scripts/import_spw_bathymetry.py \
+  --base-url http://127.0.0.1:8000 \
+  --project-name "Belgium and North Sea Workbench" \
+  --area "RC Golden - Wallonia urban-rural" \
+  --bbox 4.85,50.45,4.87,50.47 \
+  --raw-zip /app/storage/operator-evidence/spw-bathymetry/2023-05-23/raw/BATHY_50CM_ALTITUDE_DNG_GEOTIFF_3812.zip \
+  --output-dir /app/storage/operator-evidence/spw-bathymetry/2023-05-23/derived
+```
+
+The script validates the pinned official checksum, safe archive members,
+EPSG:3812, one Float32 band, approximately 0.5 m cells and nodata `-9999`.
+It then creates a bounded COG and uploads it through `/datasets/upload`.
+`POST .../raster/bathymetry/select` returns waterbed elevation in mDNG,
+surveyed surface and coverage. Current depth, volume and datum conversion stay
+unavailable without a compatible water-surface source. Selection analysis is
+bounded by `BATHYMETRY_RASTER_MAX_PIXELS` (30 million by default).
+
 ## Governed regional official-vector acquisition
 
 The thematic raster registry includes forest and agricultural land-use masks

backend/app/api/routes/datasets.py

@@ -48,6 +48,8 @@ from app.schemas import (
     FloodHazardSelectionRequest,
     BathymetryPartitionFinalizeRequest,
     BathymetryProfileAcquireRequest,
+    BathymetryRasterSelectionRequest,
+    BathymetryRasterSelectionResponse,
     ThematicRasterAcquireRequest,
     ThematicRasterProductRead,
     ThematicRasterSelectionResponse,
@@ -89,6 +91,7 @@ from app.services.terrain_analysis_service import TerrainAnalysisService
 from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
 from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService
 from app.services.bathymetry_profile_acquisition_service import BathymetryProfileAcquisitionService
+from app.services.bathymetry_raster_analysis_service import BathymetryRasterAnalysisService
 from app.services.mdk_bathymetry_probe_service import MdkBathymetryProbeService
 from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService
 from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService
@@ -888,6 +891,33 @@ def raster_terrain_image(
     )
 
 
+@router.post(
+    "/datasets/{dataset_id}/raster/bathymetry/select",
+    response_model=Envelope[BathymetryRasterSelectionResponse],
+)
+def raster_bathymetry_selection(
+    project_id: UUID,
+    dataset_id: UUID,
+    payload: BathymetryRasterSelectionRequest,
+    db: Session = Depends(get_db),
+):
+    return envelope(BathymetryRasterAnalysisService.analyze(db, project_id, dataset_id, payload))
+
+
+@router.get("/datasets/{dataset_id}/raster/bathymetry/image")
+def raster_bathymetry_image(
+    project_id: UUID,
+    dataset_id: UUID,
+    db: Session = Depends(get_db),
+):
+    content = BathymetryRasterAnalysisService.render_png(db, project_id, dataset_id)
+    return Response(
+        content=content,
+        media_type="image/png",
+        headers={"Cache-Control": "private, max-age=86400"},
+    )
+
+
 @router.post(
     "/datasets/{dataset_id}/raster/flood-hazard/select",
     response_model=Envelope[FloodHazardSelectionResponse],

backend/app/core/config.py

@@ -220,6 +220,11 @@ class Settings(BaseSettings):
         le=256,
         validation_alias="BATHYMETRY_PROFILES_MAX_RESPONSE_MB",
     )
+    bathymetry_raster_max_pixels: int = Field(
+        default=30_000_000,
+        ge=1,
+        validation_alias="BATHYMETRY_RASTER_MAX_PIXELS",
+    )
     mdk_bathymetry_probe_enabled: bool = Field(default=True, validation_alias="MDK_BATHYMETRY_PROBE_ENABLED")
     mdk_bathymetry_wcs_url: str = Field(
         default="https://bathy.agentschapmdk.be/spatialfusionserver/services/ows/wcs/EL_wcs",

backend/app/schemas/__init__.py

@@ -93,6 +93,10 @@ from .bathymetry import (
     BathymetryPartitionFinalizationResult,
     BathymetryProfileAcquireRequest,
     BathymetryProfileAcquisitionResult,
+    BathymetryRasterMetric,
+    BathymetryRasterSelectionRequest,
+    BathymetryRasterSelectionResponse,
+    BathymetryRasterSelectionSummary,
     BathymetrySourceProbeRead,
     BathymetrySourceRead,
 )

backend/app/schemas/bathymetry.py

@@ -109,3 +109,43 @@ class BathymetrySourceProbeRead(BaseModel):
     checked_at: datetime
     message: str
     limitation_message: str
+
+
+class BathymetryRasterSelectionRequest(BaseModel):
+    bbox: VectorSelectionBBox
+    area_id: UUID | None = None
+
+
+class BathymetryRasterMetric(BaseModel):
+    metric_key: str
+    metric_label: str
+    metric_value: float
+    metric_unit: str
+    aggregation_method: str
+    is_estimate: bool = False
+
+
+class BathymetryRasterSelectionSummary(BaseModel):
+    metric_label: str
+    metric_value: float
+    metric_unit: str
+    aggregation_method: str
+    primary_metric_key: str
+    metrics: list[BathymetryRasterMetric]
+
+
+class BathymetryRasterSelectionResponse(BaseModel):
+    dataset_id: UUID
+    product_key: str
+    selection_bbox: VectorSelectionBBox
+    selection_area_id: UUID | None = None
+    selected_cell_count: int = Field(ge=1)
+    valid_cell_count: int = Field(ge=1)
+    coverage_ratio: float = Field(ge=0, le=1)
+    resolution_m: float = Field(gt=0)
+    vertical_reference: str
+    survey_period: str
+    summary: BathymetryRasterSelectionSummary
+    unsupported_metrics: list[str]
+    limitation_message: str
+    generated_at: str

backend/app/services/bathymetry_profile_acquisition_service.py

@@ -91,20 +91,21 @@ class BathymetryProfileAcquisitionService:
             "authority_level": "authoritative",
             "geographic_coverage": "Waalse bevaarbare waterwegen en stuwmeren met uitgevoerde opmetingen",
             "data_kind": "bodemhoogteraster en XYZ-puntenwolk",
-            "query_modes": ["download", "arcgis_map_service"],
+            "query_modes": ["operator_archive", "bounded_raster", "arcgis_map_service"],
             "vertical_reference": "mDNG",
             "horizontal_crs": "EPSG:3812; visualisatieservice kan EPSG:31370 aanbieden",
             "native_resolution": "0,5 m",
-            "integration_status": "available_not_integrated",
-            "acquisition_supported": False,
-            "configured": False,
+            "integration_status": "operational",
+            "acquisition_supported": True,
+            "configured": True,
             "service_url": "https://geoservices.wallonie.be/arcgis/rest/services/EAU/BATHY/MapServer",
-            "catalog_url": "https://geoportail.wallonie.be/catalogue/c450c28f-d357-48af-8423-62d524632cf9.html",
+            "catalog_url": "https://geoportail.wallonie.be/catalogue/0a544b42-0b30-4c8e-85e7-38149b99eae0.html",
             "attribution": "Service public de Wallonie",
             "license_note": "CC BY 4.0 volgens de officiële Geoportail-metadata.",
             "limitation_message": (
-                "Dekking en meetjaar verschillen per vaarweg of reservoir. Integratie vereist een beheerde "
-                "download- en mosaïekstroom plus expliciete omzetting van mDNG."
+                "De gepinde officiële release kan begrensd als raster worden geïmporteerd via de operator. "
+                "Dekking verschilt per vaarweg; de waarden zijn bodemhoogtes in mDNG uit 2019-2022, "
+                "zonder stilzwijgende datumconversie of afleiding van actuele waterdiepte."
             ),
         },
         {

backend/app/services/bathymetry_raster_analysis_service.py

@@ -0,0 +1,359 @@
+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.bathymetry import (
+    BathymetryRasterMetric,
+    BathymetryRasterSelectionRequest,
+    BathymetryRasterSelectionResponse,
+    BathymetryRasterSelectionSummary,
+)
+
+
+class BathymetryRasterAnalysisService:
+    SOURCE_NAME = "spw_bathymetry"
+    PRODUCT_KEY = "spw_bathymetry_50cm_mdng"
+    UNSUPPORTED_METRICS = [
+        "current_water_depth_m",
+        "water_volume_m3",
+        "vertical_datum_conversion",
+    ]
+    LIMITATION = (
+        "De rasterwaarden zijn waterbodemhoogtes in mDNG uit een samengestelde SPW-opmeting "
+        "(2019-2022). Zonder een gelijktijdig waterpeil zijn actuele waterdiepte en watervolume "
+        "niet berekenbaar. mDNG wordt niet stilzwijgend naar TAW, LAT of een ander verticaal datum omgezet."
+    )
+
+    @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 != BathymetryRasterAnalysisService.SOURCE_NAME:
+            raise AppError(
+                code="INVALID_BATHYMETRY_RASTER_DATASET",
+                message="Bathymetry analysis requires a governed SPW bathymetry 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 bathymetry raster file is unavailable",
+                status_code=404,
+            )
+        return dataset
+
+    @staticmethod
+    def _metadata(dataset: Dataset) -> dict:
+        metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
+        if (
+            metadata.get("product_key") != BathymetryRasterAnalysisService.PRODUCT_KEY
+            or metadata.get("theme") != "bathymetry"
+            or metadata.get("value_semantics") != "bed_elevation"
+            or metadata.get("vertical_reference") != "mDNG"
+            or metadata.get("source_crs") != "EPSG:3812"
+        ):
+            raise AppError(
+                code="INVALID_BATHYMETRY_RASTER_METADATA",
+                message="Bathymetry raster provenance or value semantics are incomplete",
+                status_code=409,
+            )
+        return metadata
+
+    @staticmethod
+    def _selection_geometry(db, project_id: UUID, payload: BathymetryRasterSelectionRequest):
+        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="BATHYMETRY_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: BathymetryRasterSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        dataset = BathymetryRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        source_metadata = BathymetryRasterAnalysisService._metadata(dataset)
+        selection_4326 = BathymetryRasterAnalysisService._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 bathymetry analysis",
+                status_code=503,
+            ) from exc
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                if source.crs is None or source.crs.to_epsg() != 3812:
+                    raise AppError(
+                        code="INVALID_DATASET_CRS",
+                        message="SPW bathymetry raster CRS must be EPSG:3812",
+                        status_code=409,
+                    )
+                if source.count != 1:
+                    raise AppError(
+                        code="INVALID_BATHYMETRY_RASTER_BANDS",
+                        message="SPW bathymetry requires one bed-elevation band",
+                        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="BATHYMETRY_SELECTION_OUTSIDE_DATASET",
+                        message="Selection does not overlap the persisted bathymetry 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.bathymetry_raster_max_pixels:
+                    raise AppError(
+                        code="BATHYMETRY_SELECTION_TOO_LARGE",
+                        message="Bathymetry analysis exceeds the configured raster cell limit",
+                        details={
+                            "pixel_count": expected_cells,
+                            "max_pixels": resolved_settings.bathymetry_raster_max_pixels,
+                        },
+                        status_code=422,
+                    )
+                clipped, clipped_transform = mask(
+                    source,
+                    [mapping(analysis_geometry)],
+                    crop=True,
+                    filled=False,
+                    indexes=[1],
+                )
+                band = np.ma.asarray(clipped[0], dtype="float64")
+                raw = band.filled(np.nan)
+                selected_cells = geometry_mask(
+                    [mapping(analysis_geometry)],
+                    out_shape=band.shape,
+                    transform=clipped_transform,
+                    invert=True,
+                )
+                valid_cells = selected_cells & ~np.ma.getmaskarray(band) & np.isfinite(raw)
+                if source.nodata is not None:
+                    valid_cells &= ~np.isclose(raw, float(source.nodata))
+                values = raw[valid_cells]
+                if values.size == 0:
+                    raise AppError(
+                        code="BATHYMETRY_NO_VALID_DATA",
+                        message="No surveyed waterbed cells occur in this selection",
+                        status_code=422,
+                    )
+                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="BATHYMETRY_ANALYSIS_FAILED",
+                message="The persisted bathymetry 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) -> BathymetryRasterMetric:
+            return BathymetryRasterMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        selected_cell_count = int(selected_cells.sum())
+        valid_cell_count = int(values.size)
+        vertical_unit = str(source_metadata["vertical_reference"])
+        coverage_ratio = valid_cell_count / max(1, selected_cell_count)
+        metrics = [
+            metric(
+                "bed_elevation_mean_m",
+                "Gemiddelde waterbodemhoogte",
+                values.mean(),
+                f"m {vertical_unit}",
+                "mean_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_min_m",
+                "Laagste waterbodemhoogte",
+                values.min(),
+                f"m {vertical_unit}",
+                "minimum_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_max_m",
+                "Hoogste waterbodemhoogte",
+                values.max(),
+                f"m {vertical_unit}",
+                "maximum_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_p10_m",
+                "10e percentiel waterbodemhoogte",
+                np.percentile(values, 10),
+                f"m {vertical_unit}",
+                "percentile_10_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_p90_m",
+                "90e percentiel waterbodemhoogte",
+                np.percentile(values, 90),
+                f"m {vertical_unit}",
+                "percentile_90_valid_source_cells",
+            ),
+            metric(
+                "surveyed_bed_surface_ha",
+                "Oppervlakte met gemeten waterbodem",
+                valid_cell_count * cell_area_m2 / 10_000.0,
+                "ha",
+                "valid_source_cells_times_cell_area",
+            ),
+            metric(
+                "bathymetry_coverage_pct",
+                "Dekking waterbodemmeting",
+                coverage_ratio * 100.0,
+                "%",
+                "valid_source_cells_divided_by_selected_cells",
+            ),
+        ]
+        primary = metrics[0]
+        response = BathymetryRasterSelectionResponse(
+            dataset_id=dataset.id,
+            product_key=BathymetryRasterAnalysisService.PRODUCT_KEY,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            selected_cell_count=selected_cell_count,
+            valid_cell_count=valid_cell_count,
+            coverage_ratio=round(coverage_ratio, 6),
+            resolution_m=round(max(resolution_x, resolution_y), 4),
+            vertical_reference=vertical_unit,
+            survey_period=str(source_metadata.get("survey_period") or "2019-2022"),
+            summary=BathymetryRasterSelectionSummary(
+                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=BathymetryRasterAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=BathymetryRasterAnalysisService.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 = BathymetryRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        BathymetryRasterAnalysisService._metadata(dataset)
+        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 bathymetry 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 source.nodata is not None:
+                    valid &= ~np.isclose(values, float(source.nodata))
+                if not valid.any():
+                    raise AppError(
+                        code="BATHYMETRY_NO_VALID_DATA",
+                        message="Bathymetry 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)
+                normalized = np.where(valid, normalized, 0.0)
+                stops = np.asarray([0.0, 0.35, 0.7, 1.0])
+                colors = np.asarray(
+                    [
+                        [8, 47, 73],
+                        [15, 118, 140],
+                        [103, 190, 170],
+                        [236, 224, 163],
+                    ],
+                    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, 220, 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="BATHYMETRY_PREVIEW_FAILED",
+                message="The persisted bathymetry raster could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc

backend/app/services/coverage_registry_service.py

@@ -233,11 +233,11 @@ SOURCE_DEFINITIONS = (
         attribution="Service public de Wallonie",
         license_note="Consult the license of each Geoportail Wallonie product.",
         limitation_message=(
-            "Bounded PICC buildings, road axes and hydrography are operational; "
+            "Bounded PICC buildings, road axes, hydrography and operator-imported SPW bathymetry are operational; "
             "other Walloon themes remain unavailable until separately governed."
         ),
-        materialized_source_names=("spw_picc",),
-        operational_themes=("buildings", "roads", "surface_water"),
+        materialized_source_names=("spw_picc", "spw_bathymetry"),
+        operational_themes=("buildings", "roads", "surface_water", "bathymetry"),
     ),
     _contract(
         source_name="urbis",
@@ -348,15 +348,16 @@ FLANDERS_THEME_DATASETS: dict[str, dict[str, tuple[str, ...]]] = {
     "flood_climate": {"vmm_flood_hazard": ()},
 }
 
-REGIONAL_THEME_DATASETS: dict[str, dict[str, tuple[str, ...]]] = {
+REGIONAL_THEME_DATASETS: dict[str, dict[str, dict[str, tuple[str, ...]]]] = {
     "spw_geoportail": {
-        "buildings": ("buildings",),
-        "roads": ("roads",),
-        "surface_water": ("water",),
+        "buildings": {"spw_picc": ("buildings",)},
+        "roads": {"spw_picc": ("roads",)},
+        "surface_water": {"spw_picc": ("water",)},
+        "bathymetry": {"spw_bathymetry": ()},
     },
     "urbis": {
-        "buildings": ("buildings",),
-        "parcels": ("parcels",),
+        "buildings": {"urbis": ("buildings",)},
+        "parcels": {"urbis": ("parcels",)},
     },
 }
 
@@ -437,10 +438,10 @@ class CoverageRegistryService:
                     continue
                 layer_names = theme_sources[dataset.source_name]
             elif definition.contract.source_name in REGIONAL_THEME_DATASETS:
-                theme_layers = REGIONAL_THEME_DATASETS[definition.contract.source_name].get(theme)
-                if theme_layers is None:
+                theme_sources = REGIONAL_THEME_DATASETS[definition.contract.source_name].get(theme, {})
+                if dataset.source_name not in theme_sources:
                     continue
-                layer_names = theme_layers
+                layer_names = theme_sources[dataset.source_name]
             metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
             coverage_zones = metadata.get("coverage_zones") or metadata.get("coverage_zone") or []
             if isinstance(coverage_zones, str):

backend/tests/test_rc4_national_coverage.py

@@ -221,6 +221,53 @@ def test_bounded_api_materialization_only_covers_its_persisted_bbox() -> None:
     assert outside.items[0].materialized_dataset_ids == []
 
 
+def test_spw_bathymetry_materialization_is_source_specific() -> None:
+    project_id = uuid4()
+    scope = [
+        scope_area("Belgium land", box(2.5, 49.5, 6.4, 51.5)),
+        scope_area("Wallonia", box(2.5, 49.5, 6.4, 50.8)),
+    ]
+    selection = CoverageBBox(minx=4.851, miny=50.451, maxx=4.869, maxy=50.469)
+    spw_picc = SimpleNamespace(
+        id=uuid4(),
+        status="ready",
+        source_name="spw_picc",
+        reference_layer_name="buildings",
+        source_metadata={
+            "coverage_zones": ["wallonia"],
+            "bbox_epsg4326": [4.85, 50.45, 4.87, 50.47],
+        },
+    )
+    without_bathymetry = CoverageRegistryService.resolve(
+        FakeSession(project=SimpleNamespace(id=project_id), areas=scope, datasets=[spw_picc]),
+        project_id,
+        selection,
+        ["bathymetry"],
+    )
+    assert without_bathymetry.items[0].status == "partial"
+    assert without_bathymetry.items[0].materialized_dataset_ids == []
+
+    bathymetry_id = uuid4()
+    bathymetry = SimpleNamespace(
+        id=bathymetry_id,
+        status="ready",
+        source_name="spw_bathymetry",
+        reference_layer_name=None,
+        source_metadata={
+            "coverage_zones": ["wallonia"],
+            "bbox_epsg4326": [4.85, 50.45, 4.87, 50.47],
+        },
+    )
+    with_bathymetry = CoverageRegistryService.resolve(
+        FakeSession(project=SimpleNamespace(id=project_id), areas=scope, datasets=[spw_picc, bathymetry]),
+        project_id,
+        selection,
+        ["bathymetry"],
+    )
+    assert with_bathymetry.items[0].status == "operational"
+    assert with_bathymetry.items[0].materialized_dataset_ids == [bathymetry_id]
+
+
 def test_mixed_land_and_north_sea_selection_remains_split() -> None:
     project_id = uuid4()
     project = SimpleNamespace(id=project_id)

backend/tests/test_sprint241_spw_bathymetry_raster.py

@@ -0,0 +1,270 @@
+from __future__ import annotations
+
+import importlib.util
+import io
+from pathlib import Path
+import sys
+import zipfile
+from uuid import uuid4
+
+import numpy as np
+import pytest
+import rasterio
+from fastapi.testclient import TestClient
+from pyproj import Transformer
+from rasterio.io import MemoryFile
+from rasterio.transform import from_origin
+from shapely.geometry import shape
+
+from app.core.config import Settings
+from app.core.errors import AppError
+from app.db.session import get_db
+from app.main import app
+from app.models import Dataset
+from app.schemas.bathymetry import BathymetryRasterSelectionRequest
+from app.services.bathymetry_raster_analysis_service import BathymetryRasterAnalysisService
+
+
+ROOT = Path(__file__).resolve().parents[2]
+SCRIPT_PATH = ROOT / "scripts" / "import_spw_bathymetry.py"
+
+
+def load_operator():
+    name = "test_import_spw_bathymetry_sprint241"
+    spec = importlib.util.spec_from_file_location(name, SCRIPT_PATH)
+    assert spec is not None
+    assert spec.loader is not None
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+OPERATOR = load_operator()
+
+
+class FakeSession:
+    def __init__(self, rows):
+        self.rows = rows
+
+    def get(self, model, row_id):
+        return self.rows.get((model, row_id))
+
+
+def bathymetry_tiff(*, nodata_only: bool = False) -> bytes:
+    values = np.linspace(72.0, 80.0, 400, dtype="float32").reshape(20, 20)
+    values[:, :5] = -9999.0
+    if nodata_only:
+        values[:] = -9999.0
+    with MemoryFile() as memory:
+        with memory.open(
+            driver="GTiff",
+            width=20,
+            height=20,
+            count=1,
+            dtype="float32",
+            crs="EPSG:3812",
+            transform=from_origin(684_000, 629_000, 0.5, 0.5),
+            nodata=-9999.0,
+        ) as output:
+            output.write(values, 1)
+        return memory.read()
+
+
+def selection_payload() -> BathymetryRasterSelectionRequest:
+    transformer = Transformer.from_crs("EPSG:3812", "EPSG:4326", always_xy=True)
+    min_x, min_y = transformer.transform(684_000, 628_990)
+    max_x, max_y = transformer.transform(684_010, 629_000)
+    return BathymetryRasterSelectionRequest(
+        bbox={
+            "min_x": min(min_x, max_x),
+            "min_y": min(min_y, max_y),
+            "max_x": max(min_x, max_x),
+            "max_y": max(min_y, max_y),
+            "crs": "EPSG:4326",
+        }
+    )
+
+
+def persisted_dataset(path: Path, *, metadata: dict | None = None) -> Dataset:
+    path.write_bytes(bathymetry_tiff())
+    return Dataset(
+        id=uuid4(),
+        project_id=uuid4(),
+        name="spw_bathymetry_test_3812.tif",
+        dataset_type="raster",
+        source="SPW official operator archive",
+        source_name="spw_bathymetry",
+        source_metadata=metadata
+        or {
+            "product_key": "spw_bathymetry_50cm_mdng",
+            "theme": "bathymetry",
+            "value_semantics": "bed_elevation",
+            "vertical_reference": "mDNG",
+            "source_crs": "EPSG:3812",
+            "survey_period": "2019-2022",
+        },
+        storage_path=str(path),
+        status="ready",
+    )
+
+
+def test_bathymetry_analysis_returns_real_bed_elevation_and_surface_metrics(tmp_path: Path) -> None:
+    dataset = persisted_dataset(tmp_path / "bathymetry.tif")
+    db = FakeSession({(Dataset, dataset.id): dataset})
+
+    result = BathymetryRasterAnalysisService.analyze(
+        db,
+        dataset.project_id,
+        dataset.id,
+        selection_payload(),
+        settings=Settings(_env_file=None, bathymetry_raster_max_pixels=10_000),
+    )
+
+    metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]}
+    assert result["product_key"] == "spw_bathymetry_50cm_mdng"
+    assert result["vertical_reference"] == "mDNG"
+    assert result["survey_period"] == "2019-2022"
+    assert result["selected_cell_count"] == 400
+    assert result["valid_cell_count"] == 300
+    assert result["coverage_ratio"] == pytest.approx(0.75)
+    assert metrics["bed_elevation_mean_m"]["metric_unit"] == "m mDNG"
+    assert metrics["surveyed_bed_surface_ha"]["metric_value"] == pytest.approx(0.0075)
+    assert metrics["bathymetry_coverage_pct"]["metric_value"] == pytest.approx(75.0)
+    assert result["unsupported_metrics"] == [
+        "current_water_depth_m",
+        "water_volume_m3",
+        "vertical_datum_conversion",
+    ]
+
+
+def test_bathymetry_analysis_fails_closed_for_metadata_size_and_empty_cells(tmp_path: Path) -> None:
+    dataset = persisted_dataset(tmp_path / "bathymetry.tif")
+    db = FakeSession({(Dataset, dataset.id): dataset})
+
+    with pytest.raises(AppError) as size_error:
+        BathymetryRasterAnalysisService.analyze(
+            db,
+            dataset.project_id,
+            dataset.id,
+            selection_payload(),
+            settings=Settings(_env_file=None, bathymetry_raster_max_pixels=100),
+        )
+    assert size_error.value.code == "BATHYMETRY_SELECTION_TOO_LARGE"
+
+    dataset.source_metadata = {"theme": "bathymetry"}
+    with pytest.raises(AppError) as metadata_error:
+        BathymetryRasterAnalysisService.analyze(
+            db,
+            dataset.project_id,
+            dataset.id,
+            selection_payload(),
+        )
+    assert metadata_error.value.code == "INVALID_BATHYMETRY_RASTER_METADATA"
+
+    dataset.source_metadata = {
+        "product_key": "spw_bathymetry_50cm_mdng",
+        "theme": "bathymetry",
+        "value_semantics": "bed_elevation",
+        "vertical_reference": "mDNG",
+        "source_crs": "EPSG:3812",
+    }
+    Path(dataset.storage_path).write_bytes(bathymetry_tiff(nodata_only=True))
+    with pytest.raises(AppError) as empty_error:
+        BathymetryRasterAnalysisService.analyze(
+            db,
+            dataset.project_id,
+            dataset.id,
+            selection_payload(),
+        )
+    assert empty_error.value.code == "BATHYMETRY_NO_VALID_DATA"
+
+
+def test_bathymetry_image_and_route_use_persisted_raster_and_canonical_envelope(tmp_path: Path) -> None:
+    dataset = persisted_dataset(tmp_path / "bathymetry.tif")
+    db = FakeSession({(Dataset, dataset.id): dataset})
+    image = BathymetryRasterAnalysisService.render_png(db, dataset.project_id, dataset.id)
+
+    assert image.startswith(b"\x89PNG\r\n\x1a\n")
+
+    app.dependency_overrides[get_db] = lambda: db
+    try:
+        response = TestClient(app).post(
+            f"/api/v1/projects/{dataset.project_id}/datasets/{dataset.id}/raster/bathymetry/select",
+            json=selection_payload().model_dump(mode="json"),
+        )
+    finally:
+        app.dependency_overrides.clear()
+    assert response.status_code == 200
+    payload = response.json()
+    assert set(payload) == {"data"}
+    assert payload["data"]["dataset_id"] == str(dataset.id)
+    assert payload["data"]["summary"]["primary_metric_key"] == "bed_elevation_mean_m"
+
+
+def test_operator_validates_pinned_archive_and_rejects_unsafe_members(
+    tmp_path: Path,
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
+    safe_path = tmp_path / "safe.zip"
+    with zipfile.ZipFile(safe_path, "w", compression=zipfile.ZIP_DEFLATED) as archive:
+        archive.writestr(OPERATOR.SOURCE_MEMBER, bathymetry_tiff())
+    monkeypatch.setattr(OPERATOR, "SOURCE_SHA256", OPERATOR.sha256_file(safe_path))
+
+    member = OPERATOR.validate_archive(safe_path)
+
+    assert member.filename == OPERATOR.SOURCE_MEMBER
+
+    unsafe_path = tmp_path / "unsafe.zip"
+    with zipfile.ZipFile(unsafe_path, "w", compression=zipfile.ZIP_DEFLATED) as archive:
+        archive.writestr("../escape.txt", "unsafe")
+        archive.writestr(OPERATOR.SOURCE_MEMBER, bathymetry_tiff())
+    monkeypatch.setattr(OPERATOR, "SOURCE_SHA256", OPERATOR.sha256_file(unsafe_path))
+    with pytest.raises(OPERATOR.SpwBathymetryImportError, match="unsafe member"):
+        OPERATOR.validate_archive(unsafe_path)
+
+
+def test_operator_crops_zip_member_to_valid_cog(tmp_path: Path, monkeypatch: pytest.MonkeyPatch) -> None:
+    archive_path = tmp_path / "source.zip"
+    with zipfile.ZipFile(archive_path, "w", compression=zipfile.ZIP_DEFLATED) as archive:
+        archive.writestr(OPERATOR.SOURCE_MEMBER, bathymetry_tiff())
+    monkeypatch.setattr(OPERATOR, "SOURCE_SHA256", OPERATOR.sha256_file(archive_path))
+    member = OPERATOR.validate_archive(archive_path)
+    output_path = tmp_path / "bounded.tif"
+
+    diagnostics = OPERATOR.crop_source(
+        archive_path,
+        member,
+        shape(
+            {
+                "type": "Polygon",
+                "coordinates": [[
+                    [selection_payload().bbox.min_x, selection_payload().bbox.min_y],
+                    [selection_payload().bbox.max_x, selection_payload().bbox.min_y],
+                    [selection_payload().bbox.max_x, selection_payload().bbox.max_y],
+                    [selection_payload().bbox.min_x, selection_payload().bbox.max_y],
+                    [selection_payload().bbox.min_x, selection_payload().bbox.min_y],
+                ]],
+            }
+        ),
+        output_path,
+        max_pixels=10_000,
+    )
+
+    with rasterio.open(output_path) as output:
+        assert output.crs.to_epsg() == 3812
+        assert output.driver == "GTiff"
+        assert output.nodata == -9999.0
+        assert output.profile["tiled"] is True
+    assert diagnostics["valid_cell_count"] == 300
+    assert len(diagnostics["output_sha256"]) == 64
+
+
+def test_operator_is_api_only_and_does_not_claim_depth_or_volume() -> None:
+    source = SCRIPT_PATH.read_text(encoding="utf-8")
+
+    assert "/datasets/upload" in source
+    assert "water_depth_available" in source
+    assert '"water_volume_available": False' in source
+    assert "SessionLocal" not in source
+    assert "db.add(" not in source