feat: complete Wallonia land cover and terrain sources

backend/app/services/coverage_registry_service.py

@@ -253,10 +253,10 @@ SOURCE_DEFINITIONS = (
         license_note="Consult the license of each Geoportail Wallonie product.",
         limitation_message=(
             "Bounded PICC buildings, road axes and hydrography, the legally current flood-hazard polygons, "
-            "and operator-imported SPW bathymetry are operational; other Walloon themes remain separately governed."
+            "operator-imported SPW bathymetry and bounded SPW MNT terrain are operational; other Walloon themes remain separately governed."
         ),
-        materialized_source_names=("spw_picc", "spw_flood_hazard", "spw_walous_land_cover", "spw_bathymetry"),
-        operational_themes=("buildings", "roads", "surface_water", "land_cover_use", "flood_climate", "bathymetry"),
+        materialized_source_names=("spw_picc", "spw_flood_hazard", "spw_walous_land_cover", "spw_bathymetry", "spw_terrain"),
+        operational_themes=("buildings", "roads", "surface_water", "land_cover_use", "elevation", "flood_climate", "bathymetry"),
     ),
     _contract(
         source_name="urbis",
@@ -376,6 +376,7 @@ REGIONAL_THEME_DATASETS: dict[str, dict[str, dict[str, tuple[str, ...]]]] = {
         "roads": {"spw_picc": ("roads",)},
         "surface_water": {"spw_picc": ("water",)},
         "land_cover_use": {"spw_walous_land_cover": ()},
+        "elevation": {"spw_terrain": ()},
         "flood_climate": {"spw_flood_hazard": ("flood_hazard",)},
         "bathymetry": {"spw_bathymetry": ()},
     },

backend/app/services/model_registry_service.py

@@ -5,7 +5,10 @@ from typing import Type
 
 from app.core.config import Settings, get_settings
 from app.schemas.detection import DetectionModelCapability
-from app.services.segmentation_adapter import SamSegmentationAdapter, YoloSegmentationAdapter
+from app.services.segmentation_adapter import (
+    SamSegmentationAdapter,
+    YoloSegmentationAdapter,
+)
 from app.services.yolo_adapter import YoloDetectionAdapter
 
 
@@ -39,7 +42,9 @@ class ModelRegistryService:
                 limitation_message="YOLO/PyTorch inference is not configured in Sprint 8; no model is downloaded or executed.",
                 version=None,
             ),
-            ModelRegistryService._configured_yolo_capability(resolved_settings, yolo_adapter_class),
+            ModelRegistryService._configured_yolo_capability(
+                resolved_settings, yolo_adapter_class
+            ),
             DetectionModelCapability(
                 model_id="manual-fixture-detector",
                 display_name="Manual fixture detector",
@@ -104,8 +109,12 @@ class ModelRegistryService:
                 limitation_message="Fixture segmenter is for explicit tests/demo fixtures only and is not production inference.",
                 version="fixture-v1",
             ),
-            ModelRegistryService._configured_yolo_seg_capability(resolved_settings, yolo_seg_adapter_class),
-            ModelRegistryService._configured_sam_capability(resolved_settings, sam_adapter_class),
+            ModelRegistryService._configured_yolo_seg_capability(
+                resolved_settings, yolo_seg_adapter_class
+            ),
+            ModelRegistryService._configured_sam_capability(
+                resolved_settings, sam_adapter_class
+            ),
         ]
 
     @staticmethod
@@ -119,7 +128,11 @@ class ModelRegistryService:
             "YOLO segmentation is disabled. Set YOLO_SEG_ENABLED=true and YOLO_SEG_MODEL_PATH to a local "
             "segmentation model file to enable inference. GeoIntel never downloads model weights automatically."
         )
-        model_path = Path(settings.yolo_seg_model_path).expanduser() if settings.yolo_seg_model_path else None
+        model_path = (
+            Path(settings.yolo_seg_model_path).expanduser()
+            if settings.yolo_seg_model_path
+            else None
+        )
 
         if settings.yolo_seg_enabled:
             if not adapter_class.dependencies_available():
@@ -157,7 +170,11 @@ class ModelRegistryService:
             "SAM is disabled. Set SAM_ENABLED=true and SAM_MODEL_PATH to a local SAM-compatible model file to "
             "enable class-agnostic segmentation. GeoIntel never downloads model weights automatically."
         )
-        model_path = Path(settings.sam_model_path).expanduser() if settings.sam_model_path else None
+        model_path = (
+            Path(settings.sam_model_path).expanduser()
+            if settings.sam_model_path
+            else None
+        )
 
         if settings.sam_enabled:
             if not adapter_class.dependencies_available():
@@ -192,7 +209,11 @@ class ModelRegistryService:
         configured = False
         status = "not_configured"
         limitation = "YOLO is disabled. Set YOLO_ENABLED=true and YOLO_MODEL_PATH to a local model file to enable inference."
-        model_path = Path(settings.yolo_model_path).expanduser() if settings.yolo_model_path else None
+        model_path = (
+            Path(settings.yolo_model_path).expanduser()
+            if settings.yolo_model_path
+            else None
+        )
 
         if settings.yolo_enabled:
             if not yolo_adapter_class.dependencies_available():
@@ -217,4 +238,11 @@ class ModelRegistryService:
             status=status,
             limitation_message=limitation,
             version=settings.yolo_model_version,
+            training_scope=(
+                "Operator-managed local weights; the runtime has no nationally governed training-corpus evidence."
+            ),
+            validation_scope="Mol and the Kempen operator evidence; no Belgian national validation matrix is bound.",
+            validated_regions=["flanders_mol_kempen"],
+            nationally_validated=False,
+            operator_review_required=True,
         )

backend/app/services/spw_terrain_service.py

@@ -0,0 +1,467 @@
+from __future__ import annotations
+
+from datetime import UTC, datetime
+import hashlib
+import json
+import math
+from pathlib import Path
+from typing import Any
+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.spw_terrain import (
+    SpwTerrainAcquireRequest,
+    SpwTerrainAcquisitionResult,
+    SpwTerrainProductRead,
+)
+from app.services.dataset_service import DatasetService
+
+
+class SpwTerrainService:
+    PROVIDER = "spw_terrain"
+    PRODUCT_KEY = "spw_mnt_1m_2021_2022"
+    DISPLAY_NAME = "SPW terreinmodel (MNT) 2021-2022"
+    SOURCE_FILENAME = "spw_mnt_1m_2021_2022_3812.tif"
+    SOURCE_SHA256_FILENAME = "spw_mnt_1m_2021_2022_3812.sha256"
+    SOURCE_CRS = "EPSG:3812"
+    SOURCE_RESOLUTION_M = 1.0
+    SURFACE_MODEL = "terrain"
+    VERTICAL_REFERENCE = "DNG / Deuxieme Nivellement General (EPSG:5710)"
+    VERTICAL_UNIT_LABEL = "m DNG"
+    ACQUISITION_PERIOD = "2021-02-19/2022-03-05"
+    CATALOG_URL = "https://geoportail.wallonie.be/catalogue/fe13bc84-e371-46ca-9632-8ad4139f1ee5.html"
+    DOWNLOAD_URL = (
+        "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/"
+        "fe13bc84-e371-46ca-9632-8ad4139f1ee5/RELIEF_WALLONIE_MNT_1M_2021_2022_GEOTIFF_3812.zip"
+    )
+    ATTRIBUTION = (
+        "Service public de Wallonie (SPW) - Relief de la Wallonie MNT 2021-2022"
+    )
+    LICENSE_NOTE = "CC BY 4.0; cite SPW and identify modifications."
+    NODATA = -9999.0
+    LIMITATION = (
+        "GeoIntel leest uitsluitend een begrensd venster uit het checksum-gevalideerde officiele 1 m MNT en "
+        "bewaart een analyse-afgeleide op de gekozen resolutie. Het MNT beschrijft maaiveldhoogte in DNG, niet "
+        "oppervlaktehoogte, afstroming, waterdiepte of watervolume. Kleine bronzones zijn door SPW geinterpoleerd."
+    )
+
+    @staticmethod
+    def _source_path(settings: Settings) -> Path:
+        return Path(settings.spw_terrain_source_dir) / SpwTerrainService.SOURCE_FILENAME
+
+    @staticmethod
+    def list_products(*, settings: Settings | None = None) -> list[dict[str, Any]]:
+        resolved = settings or get_settings()
+        configured = (
+            resolved.spw_terrain_enabled
+            and SpwTerrainService._source_path(resolved).is_file()
+        )
+        product = SpwTerrainProductRead(
+            key=SpwTerrainService.PRODUCT_KEY,
+            display_name=SpwTerrainService.DISPLAY_NAME,
+            surface_model=SpwTerrainService.SURFACE_MODEL,
+            source_filename=SpwTerrainService.SOURCE_FILENAME,
+            native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M,
+            analysis_resolution_m=resolved.spw_terrain_analysis_resolution_m,
+            source_crs=SpwTerrainService.SOURCE_CRS,
+            vertical_reference=SpwTerrainService.VERTICAL_REFERENCE,
+            acquisition_period=SpwTerrainService.ACQUISITION_PERIOD,
+            catalog_url=SpwTerrainService.CATALOG_URL,
+            attribution=SpwTerrainService.ATTRIBUTION,
+            license_note=SpwTerrainService.LICENSE_NOTE,
+            limitation_message=SpwTerrainService.LIMITATION,
+            coverage_zones=["wallonia"],
+            configured=configured,
+            status="configured" if configured else "source_not_provisioned",
+        )
+        return [product.model_dump()]
+
+    @staticmethod
+    def _scope_geometry(db, project_id: UUID, payload: SpwTerrainAcquireRequest):
+        if not db.get(Project, project_id):
+            raise AppError(
+                code="PROJECT_NOT_FOUND", message="Project not found", status_code=404
+            )
+        if payload.product_key.strip().lower() != SpwTerrainService.PRODUCT_KEY:
+            raise AppError(
+                code="SPW_TERRAIN_PRODUCT_NOT_SUPPORTED",
+                message="Select the governed SPW MNT 2021-2022 product",
+                details={"product_key": payload.product_key},
+                status_code=422,
+            )
+        if payload.bbox.crs.upper() != "EPSG:4326":
+            raise AppError(
+                code="INVALID_BBOX_CRS",
+                message="SPW terrain acquisition requires EPSG:4326",
+                status_code=400,
+            )
+        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 values[0] >= values[2]
+            or values[1] >= values[3]
+        ):
+            raise AppError(
+                code="INVALID_BBOX",
+                message="SPW terrain selection must be a finite non-empty rectangle",
+                status_code=400,
+            )
+        selection = box(*values)
+        if payload.area_id is None:
+            return selection, values
+        area = db.get(Area, payload.area_id)
+        if area is None or area.project_id != project_id:
+            raise AppError(
+                code="AREA_NOT_FOUND", message="Area not found", status_code=404
+            )
+        selection = selection.intersection(to_shape(area.geometry))
+        if selection.is_empty or selection.area <= 0:
+            raise AppError(
+                code="SPW_TERRAIN_SELECTION_OUTSIDE_AREA",
+                message="Selection does not overlap the selected work area",
+                status_code=422,
+            )
+        return selection, values
+
+    @staticmethod
+    def _read_source_window(
+        source_path: Path, scope_4326, resolution: float, settings: Settings
+    ) -> tuple[bytes, dict[str, Any]]:
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.enums import Resampling
+            from rasterio.features import geometry_mask
+            from rasterio.io import MemoryFile
+            from rasterio.transform import from_bounds
+            from rasterio.windows import from_bounds as window_from_bounds
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio and numpy are required for SPW terrain",
+                status_code=503,
+            ) from exc
+
+        scope_metric = shapely_transform(
+            Transformer.from_crs(
+                "EPSG:4326", SpwTerrainService.SOURCE_CRS, always_xy=True
+            ).transform,
+            scope_4326,
+        )
+        try:
+            with rasterio.open(source_path) as source:
+                if (
+                    source.crs is None
+                    or source.crs.to_epsg() != 3812
+                    or source.count != 1
+                ):
+                    raise AppError(
+                        code="SPW_TERRAIN_SOURCE_INVALID",
+                        message="SPW MNT must be a one-band EPSG:3812 raster",
+                        status_code=409,
+                    )
+                if not all(
+                    math.isclose(abs(float(value)), 1.0, abs_tol=0.05)
+                    for value in source.res
+                ):
+                    raise AppError(
+                        code="SPW_TERRAIN_SOURCE_INVALID",
+                        message="SPW MNT must retain the official 1 m resolution",
+                        status_code=409,
+                    )
+                clipped_geometry = scope_metric.intersection(box(*source.bounds))
+                if clipped_geometry.is_empty or clipped_geometry.area <= 0:
+                    raise AppError(
+                        code="SPW_TERRAIN_SELECTION_OUTSIDE_COVERAGE",
+                        message="Selection does not overlap SPW MNT coverage",
+                        status_code=422,
+                    )
+                min_x, min_y, max_x, max_y = clipped_geometry.bounds
+                bounds = (
+                    math.floor(min_x / resolution) * resolution,
+                    math.floor(min_y / resolution) * resolution,
+                    math.ceil(max_x / resolution) * resolution,
+                    math.ceil(max_y / resolution) * resolution,
+                )
+                width_m, height_m = bounds[2] - bounds[0], bounds[3] - bounds[1]
+                if (
+                    width_m > settings.spw_terrain_max_side_m
+                    or height_m > settings.spw_terrain_max_side_m
+                ):
+                    raise AppError(
+                        code="SPW_TERRAIN_SELECTION_TOO_LARGE",
+                        message="SPW terrain selection exceeds the configured side limit",
+                        status_code=422,
+                    )
+                width, height = (
+                    max(1, round(width_m / resolution)),
+                    max(1, round(height_m / resolution)),
+                )
+                if width * height > settings.spw_terrain_max_pixels:
+                    raise AppError(
+                        code="SPW_TERRAIN_SELECTION_TOO_LARGE",
+                        message="SPW terrain selection exceeds the configured cell limit",
+                        details={
+                            "pixel_count": width * height,
+                            "max_pixels": settings.spw_terrain_max_pixels,
+                        },
+                        status_code=422,
+                    )
+                window = window_from_bounds(*bounds, transform=source.transform)
+                band = source.read(
+                    1,
+                    window=window,
+                    out_shape=(height, width),
+                    masked=True,
+                    resampling=Resampling.bilinear,
+                )
+                output_transform = from_bounds(*bounds, width, height)
+                outside_scope = geometry_mask(
+                    [mapping(clipped_geometry)],
+                    out_shape=(height, width),
+                    transform=output_transform,
+                    invert=False,
+                )
+                values = np.asarray(np.ma.getdata(band), dtype="float32")
+                invalid = (
+                    np.ma.getmaskarray(band) | outside_scope | ~np.isfinite(values)
+                )
+                if source.nodata is not None:
+                    invalid |= np.isclose(
+                        values.astype("float64"), float(source.nodata)
+                    )
+                values[invalid] = SpwTerrainService.NODATA
+                valid = values[~invalid]
+                if valid.size == 0:
+                    raise AppError(
+                        code="SPW_TERRAIN_NO_VALID_DATA",
+                        message="SPW MNT contains no valid cells in this selection",
+                        status_code=422,
+                    )
+                if float(valid.min()) < -100.0 or float(valid.max()) > 1000.0:
+                    raise AppError(
+                        code="SPW_TERRAIN_SOURCE_INVALID_VALUES",
+                        message="SPW MNT contains implausible elevations for Wallonia",
+                        details={
+                            "minimum": float(valid.min()),
+                            "maximum": float(valid.max()),
+                        },
+                        status_code=409,
+                    )
+                profile = {
+                    "driver": "GTiff",
+                    "width": width,
+                    "height": height,
+                    "count": 1,
+                    "dtype": "float32",
+                    "crs": SpwTerrainService.SOURCE_CRS,
+                    "transform": output_transform,
+                    "nodata": SpwTerrainService.NODATA,
+                    "compress": "deflate",
+                    "predictor": 3,
+                }
+                with MemoryFile() as memory:
+                    with memory.open(**profile) as output:
+                        output.write(values, 1)
+                    content = memory.read()
+                return content, {
+                    "width": width,
+                    "height": height,
+                    "valid_pixel_count": int(valid.size),
+                    "bbox_epsg3812": list(bounds),
+                    "source_width": int(source.width),
+                    "source_height": int(source.height),
+                    "source_nodata": None
+                    if source.nodata is None
+                    else float(source.nodata),
+                    "source_resolution_m": 1.0,
+                    "analysis_resolution_m": resolution,
+                    "elevation_min_m": float(valid.min()),
+                    "elevation_max_m": float(valid.max()),
+                }
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="SPW_TERRAIN_SOURCE_READ_FAILED",
+                message="The provisioned SPW MNT could not be read",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) 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 == SpwTerrainService.PROVIDER,
+                Dataset.status == "ready",
+            )
+            .order_by(Dataset.imported_at.desc())
+            .first()
+        )
+        return (
+            candidate
+            if candidate
+            and candidate.storage_path
+            and Path(candidate.storage_path).is_file()
+            else None
+        )
+
+    @staticmethod
+    def acquire(
+        db,
+        project_id: UUID,
+        payload: SpwTerrainAcquireRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict[str, Any]:
+        resolved = settings or get_settings()
+        if not resolved.spw_terrain_enabled:
+            raise AppError(
+                code="SPW_TERRAIN_NOT_CONFIGURED",
+                message="SPW terrain bounded analysis is disabled",
+                status_code=503,
+            )
+        source_path = SpwTerrainService._source_path(resolved)
+        if not source_path.is_file():
+            raise AppError(
+                code="SPW_TERRAIN_SOURCE_NOT_PROVISIONED",
+                message="The official SPW MNT source archive has not been provisioned on this runtime",
+                details={
+                    "expected_path": str(source_path),
+                    "operator_command": "python scripts/provision_spw_terrain_source.py",
+                },
+                status_code=503,
+            )
+        scope, bbox_4326 = SpwTerrainService._scope_geometry(db, project_id, payload)
+        resolution = float(
+            payload.resolution_m or resolved.spw_terrain_analysis_resolution_m
+        )
+        identity = {
+            "product_key": SpwTerrainService.PRODUCT_KEY,
+            "bbox_epsg4326": [round(float(value), 8) for value in bbox_4326],
+            "area_id": str(payload.area_id) if payload.area_id else None,
+            "analysis_resolution_m": resolution,
+        }
+        request_hash = hashlib.sha256(
+            json.dumps(identity, sort_keys=True).encode()
+        ).hexdigest()
+        filename = f"spw_mnt_2021_2022_{request_hash[:12]}_3812.tif"
+        if not payload.force_refresh:
+            cached = SpwTerrainService._cached_dataset(db, project_id, filename)
+            if cached is not None:
+                metadata = cached.source_metadata or {}
+                return SpwTerrainAcquisitionResult(
+                    output_dataset_id=cached.id,
+                    reused=True,
+                    provider=SpwTerrainService.PROVIDER,
+                    product_key=SpwTerrainService.PRODUCT_KEY,
+                    display_name=SpwTerrainService.DISPLAY_NAME,
+                    surface_model=SpwTerrainService.SURFACE_MODEL,
+                    native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M,
+                    resolution_m=resolution,
+                    width=int((cached.metadata_json or {}).get("width", 0)),
+                    height=int((cached.metadata_json or {}).get("height", 0)),
+                    valid_pixel_count=int(metadata.get("valid_pixel_count", 0)),
+                    nodata_value=SpwTerrainService.NODATA,
+                    bbox_epsg4326=bbox_4326,
+                    bbox_epsg3812=list(metadata.get("bbox_epsg3812") or []),
+                    vertical_reference=SpwTerrainService.VERTICAL_REFERENCE,
+                    acquisition_period=SpwTerrainService.ACQUISITION_PERIOD,
+                    attribution=SpwTerrainService.ATTRIBUTION,
+                    limitation_message=SpwTerrainService.LIMITATION,
+                ).model_dump(mode="json")
+
+        content, validation = SpwTerrainService._read_source_window(
+            source_path, scope, resolution, resolved
+        )
+        source_sha256_path = source_path.with_name(
+            SpwTerrainService.SOURCE_SHA256_FILENAME
+        )
+        source_sha256 = (
+            source_sha256_path.read_text(encoding="ascii").strip().split()[0]
+            if source_sha256_path.is_file()
+            else None
+        )
+        acquired_at = datetime.now(UTC)
+        dataset = DatasetService.import_raster_bytes(
+            db,
+            project_id=project_id,
+            area_id=payload.area_id,
+            filename=filename,
+            content=content,
+            source="SPW Relief de la Wallonie MNT 2021-2022 operator-provisioned GeoTIFF",
+            source_name=SpwTerrainService.PROVIDER,
+            observed_at=datetime(2022, 3, 5, 23, 59, 59, tzinfo=UTC),
+            valid_from=datetime(2021, 2, 19, tzinfo=UTC),
+            valid_to=datetime(2022, 3, 5, 23, 59, 59, tzinfo=UTC),
+            temporal_granularity="acquisition_period",
+            source_version="RELIEF_WALLONIE_MNT_1M_2021_2022",
+            source_metadata={
+                "provider": SpwTerrainService.PROVIDER,
+                "product_key": SpwTerrainService.PRODUCT_KEY,
+                "product_display_name": SpwTerrainService.DISPLAY_NAME,
+                "surface_model": SpwTerrainService.SURFACE_MODEL,
+                "source_crs": SpwTerrainService.SOURCE_CRS,
+                "source_resolution_m": SpwTerrainService.SOURCE_RESOLUTION_M,
+                "analysis_resolution_m": validation["analysis_resolution_m"],
+                "valid_pixel_count": validation["valid_pixel_count"],
+                "bbox_epsg4326": bbox_4326,
+                "bbox_epsg3812": validation["bbox_epsg3812"],
+                "coverage_zones": ["wallonia"],
+                "vertical_reference": SpwTerrainService.VERTICAL_REFERENCE,
+                "vertical_unit": "m",
+                "vertical_unit_label": SpwTerrainService.VERTICAL_UNIT_LABEL,
+                "acquisition_period": SpwTerrainService.ACQUISITION_PERIOD,
+                "catalog_url": SpwTerrainService.CATALOG_URL,
+                "download_url": SpwTerrainService.DOWNLOAD_URL,
+                "attribution": SpwTerrainService.ATTRIBUTION,
+                "license_note": SpwTerrainService.LICENSE_NOTE,
+                "limitation_message": SpwTerrainService.LIMITATION,
+            },
+            provenance_metadata={
+                "acquisition": "operator_provisioned_official_archive_bounded_window",
+                "acquired_at": acquired_at.isoformat(),
+                "request_hash": request_hash,
+                "source_filename": SpwTerrainService.SOURCE_FILENAME,
+                "source_sha256": source_sha256,
+                "derived_sha256": hashlib.sha256(content).hexdigest(),
+                "resampling": "bilinear",
+                "validation": validation,
+            },
+        )
+        return SpwTerrainAcquisitionResult(
+            output_dataset_id=dataset.id,
+            reused=False,
+            provider=SpwTerrainService.PROVIDER,
+            product_key=SpwTerrainService.PRODUCT_KEY,
+            display_name=SpwTerrainService.DISPLAY_NAME,
+            surface_model=SpwTerrainService.SURFACE_MODEL,
+            native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M,
+            resolution_m=validation["analysis_resolution_m"],
+            width=validation["width"],
+            height=validation["height"],
+            valid_pixel_count=validation["valid_pixel_count"],
+            nodata_value=SpwTerrainService.NODATA,
+            bbox_epsg4326=bbox_4326,
+            bbox_epsg3812=validation["bbox_epsg3812"],
+            vertical_reference=SpwTerrainService.VERTICAL_REFERENCE,
+            acquisition_period=SpwTerrainService.ACQUISITION_PERIOD,
+            attribution=SpwTerrainService.ATTRIBUTION,
+            limitation_message=SpwTerrainService.LIMITATION,
+        ).model_dump(mode="json")

backend/app/services/terrain_analysis_service.py

@@ -22,10 +22,13 @@ from app.schemas.dhmv import (
     TerrainSelectionSummary,
 )
 from app.services.dhmv_acquisition_service import DhmvAcquisitionService
-from app.services.raster_partition_analysis_service import RasterPartitionAnalysisService
+from app.services.raster_partition_analysis_service import (
+    RasterPartitionAnalysisService,
+)
 
 
 class TerrainAnalysisService:
+    SUPPORTED_PROVIDERS = {DhmvAcquisitionService.PROVIDER, "spw_terrain"}
     UNSUPPORTED_METRICS = ["water_depth_m", "water_volume_m3"]
     LIMITATION = (
         "Hoogte, reliëf en helling zijn afgeleid uit DHMV II. Afstroming vraagt bijkomende hydrologische modellering. "
@@ -36,30 +39,58 @@ class TerrainAnalysisService:
     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="DATASET_NOT_FOUND", message="Dataset not found", status_code=404
+            )
+        if (
+            dataset.dataset_type != "raster"
+            or dataset.source_name not in TerrainAnalysisService.SUPPORTED_PROVIDERS
+        ):
             raise AppError(
                 code="INVALID_TERRAIN_DATASET",
-                message="Terrain analysis requires a governed DHMV raster dataset",
+                message="Terrain analysis requires a governed regional elevation 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 DHMV raster file is unavailable", status_code=404)
+        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 terrain 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)
+        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)
+            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)
+            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)
+            raise AppError(
+                code="TERRAIN_SELECTION_OUTSIDE_AREA",
+                message="Selection does not overlap the selected work area",
+                status_code=422,
+            )
         return selection
 
     @staticmethod
@@ -73,27 +104,51 @@ class TerrainAnalysisService:
     ) -> 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)
+        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
+            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)
+        product_is_governed = (
+            product_key in DhmvAcquisitionService._products()
+            if dataset.source_name == DhmvAcquisitionService.PROVIDER
+            else product_key == "spw_mnt_1m_2021_2022"
+        )
+        if not product_is_governed or surface_model not in {"terrain", "surface"}:
+            raise AppError(
+                code="INVALID_TERRAIN_METADATA",
+                message="Regional terrain product provenance is incomplete",
+                status_code=409,
+            )
+        vertical_unit_label = str(source_metadata.get("vertical_unit_label") or "m TAW")
 
         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)
+                    raise AppError(
+                        code="INVALID_DATASET_CRS",
+                        message="Terrain 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:
@@ -103,12 +158,17 @@ class TerrainAnalysisService:
                         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]))
+                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},
+                        details={
+                            "pixel_count": expected_cells,
+                            "max_pixels": resolved_settings.dhmv_max_pixels,
+                        },
                         status_code=422,
                     )
                 clipped, clipped_transform = mask(
@@ -133,14 +193,20 @@ class TerrainAnalysisService:
                 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)
+                    raise AppError(
+                        code="TERRAIN_NO_VALID_DATA",
+                        message="No valid terrain 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)
+                    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:
@@ -148,12 +214,14 @@ class TerrainAnalysisService:
         except Exception as exc:
             raise AppError(
                 code="TERRAIN_ANALYSIS_FAILED",
-                message="The persisted DHMV raster could not be analysed",
+                message="The persisted terrain 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:
+        def metric(
+            key: str, label: str, value: float, unit: str, method: str
+        ) -> TerrainMetric:
             return TerrainMetric(
                 metric_key=key,
                 metric_label=label,
@@ -163,21 +231,79 @@ class TerrainAnalysisService:
             )
 
         prefix = "terrain" if surface_model == "terrain" else "surface"
-        elevation_label = "Gemiddelde maaiveldhoogte" if surface_model == "terrain" else "Gemiddelde oppervlaktehoogte"
+        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"),
+            metric(
+                f"{prefix}_elevation_mean_m",
+                elevation_label,
+                values.mean(),
+                vertical_unit_label,
+                "mean_valid_cells",
+            ),
+            metric(
+                f"{prefix}_elevation_min_m",
+                "Laagste hoogte",
+                values.min(),
+                vertical_unit_label,
+                "minimum_valid_cells",
+            ),
+            metric(
+                f"{prefix}_elevation_max_m",
+                "Hoogste hoogte",
+                values.max(),
+                vertical_unit_label,
+                "maximum_valid_cells",
+            ),
+            metric(
+                f"{prefix}_elevation_p10_m",
+                "10e percentiel hoogte",
+                np.percentile(values, 10),
+                vertical_unit_label,
+                "percentile_10_valid_cells",
+            ),
+            metric(
+                f"{prefix}_elevation_p90_m",
+                "90e percentiel hoogte",
+                np.percentile(values, 90),
+                vertical_unit_label,
+                "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"),
+                    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]
@@ -194,7 +320,10 @@ class TerrainAnalysisService:
             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),
+            vertical_reference=str(
+                source_metadata.get("vertical_reference")
+                or DhmvAcquisitionService.VERTICAL_REFERENCE
+            ),
             summary=TerrainSelectionSummary(
                 metric_label=primary.metric_label,
                 metric_value=primary.metric_value,
@@ -204,7 +333,10 @@ class TerrainAnalysisService:
                 metrics=metrics,
             ),
             unsupported_metrics=TerrainAnalysisService.UNSUPPORTED_METRICS,
-            limitation_message=TerrainAnalysisService.LIMITATION,
+            limitation_message=str(
+                source_metadata.get("limitation_message")
+                or TerrainAnalysisService.LIMITATION
+            ),
             generated_at=datetime.now(UTC).isoformat(),
         )
         return response.model_dump(mode="json")
@@ -218,7 +350,9 @@ class TerrainAnalysisService:
         settings: Settings | None = None,
     ) -> dict:
         resolved_settings = settings or get_settings()
-        product = DhmvAcquisitionService._products().get(payload.product_key.strip().lower())
+        product = DhmvAcquisitionService._products().get(
+            payload.product_key.strip().lower()
+        )
         if product is None:
             raise AppError(
                 code="DHMV_PRODUCT_NOT_SUPPORTED",
@@ -226,7 +360,9 @@ class TerrainAnalysisService:
                 details={"product_key": payload.product_key},
                 status_code=422,
             )
-        selection_4326 = TerrainAnalysisService._selection_geometry(db, project_id, payload)
+        selection_4326 = TerrainAnalysisService._selection_geometry(
+            db, project_id, payload
+        )
         partition = RasterPartitionAnalysisService.select(
             db,
             project_id,
@@ -279,7 +415,9 @@ class TerrainAnalysisService:
             slope = np.degrees(np.arctan(np.hypot(gradient_x, gradient_y)))
             slope_values = slope[np.isfinite(slope) & valid_mask]
 
-        def metric(key: str, label: str, value: float, unit: str, method: str) -> TerrainMetric:
+        def metric(
+            key: str, label: str, value: float, unit: str, method: str
+        ) -> TerrainMetric:
             return TerrainMetric(
                 metric_key=key,
                 metric_label=label,
@@ -295,19 +433,73 @@ class TerrainAnalysisService:
             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"),
+            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"),
+                    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]
@@ -344,7 +536,9 @@ class TerrainAnalysisService:
         return response.model_dump(mode="json")
 
     @staticmethod
-    def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
+    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
@@ -352,18 +546,31 @@ class TerrainAnalysisService:
             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
+            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)
+                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)
+                    raise AppError(
+                        code="TERRAIN_NO_VALID_DATA",
+                        message="Terrain raster contains no renderable cells",
+                        status_code=422,
+                    )
                 low, high = np.percentile(values[valid], [2, 98])
                 if high <= low:
                     high = low + 1.0
@@ -381,7 +588,9 @@ class TerrainAnalysisService:
                 )
                 rgba = np.zeros((height, width, 4), dtype="uint8")
                 for channel in range(3):
-                    rgba[:, :, channel] = np.interp(normalized, stops, colors[:, channel]).astype("uint8")
+                    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)
@@ -391,7 +600,7 @@ class TerrainAnalysisService:
         except Exception as exc:
             raise AppError(
                 code="TERRAIN_PREVIEW_FAILED",
-                message="The persisted DHMV raster could not be rendered",
+                message="The persisted terrain raster could not be rendered",
                 details={"reason": str(exc)},
                 status_code=500,
             ) from exc

backend/app/services/walous_land_cover_service.py

@@ -40,7 +40,10 @@ class WalousProduct:
     catalog_url: str
     download_url: str
     source_sha256_filename: str
+    attribution: str
     accuracy_label: str
+    raw_class_crosswalk: dict[int, int] | None
+    comparability_note: str
     observation_start: datetime
     observation_end: datetime
 
@@ -54,7 +57,9 @@ class WalousLandCoverService:
     METRIC_KIND = "categorical_area"
     NODATA = 255
     ATTRIBUTION = "Service public de Wallonie (SPW), Aerospacelab S.A."
-    LICENSE_NOTE = "CC BY 4.0; cite the official SPW WALOUS edition and identify modifications."
+    LICENSE_NOTE = (
+        "CC BY 4.0; cite the official SPW WALOUS edition and identify modifications."
+    )
     LIMITATION = (
         "GeoIntel analyseert een nearest-neighbour afgeleide van het officiele 1 m WALOUS-raster op de "
         "geconfigureerde analyseresolutie. Oppervlakten zijn celgebaseerde schattingen; de kaart is landbedekking, "
@@ -88,10 +93,74 @@ class WalousLandCoverService:
         80: (78, 125, 70),
         90: (107, 164, 87),
     }
+    # The original 2018 product retains stacked two-digit codes. The official
+    # "Classe vue" legend resolves those codes to the visible top class. The
+    # only 2018-only visible class, greenhouses (62), is explicitly normalized
+    # to artificial constructions so the stable 11-class series can be used.
+    WALOUS_2018_CLASS_CROSSWALK = {
+        0: NODATA,
+        1: 1,
+        11: 1,
+        15: 1,
+        18: 1,
+        19: 1,
+        31: 1,
+        51: 1,
+        71: 1,
+        81: 1,
+        91: 1,
+        2: 2,
+        28: 2,
+        29: 2,
+        62: 2,
+        3: 3,
+        38: 3,
+        39: 3,
+        73: 3,
+        83: 3,
+        93: 3,
+        4: 4,
+        5: 5,
+        55: 5,
+        58: 5,
+        59: 5,
+        75: 5,
+        85: 5,
+        95: 5,
+        6: 6,
+        7: 7,
+        8: 8,
+        9: 9,
+        80: 80,
+        90: 90,
+    }
 
     @staticmethod
     def _products() -> dict[str, WalousProduct]:
         products = (
+            WalousProduct(
+                key="walous_land_cover_2018",
+                display_name="WALOUS landbedekking 2018",
+                observation_year=2018,
+                source_filename="walous_land_cover_2018_3812.tif",
+                source_version="WALOUS_OCS__2018",
+                catalog_url="https://geoportail.wallonie.be/catalogue/a0ad23a1-1845-4bd5-8c2f-0f62d3f1ec75.html",
+                download_url=(
+                    "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/"
+                    "a0ad23a1-1845-4bd5-8c2f-0f62d3f1ec75/WALOUS_OCS__2018_GEOTIFF_3812.zip"
+                ),
+                source_sha256_filename="walous_land_cover_2018_3812.sha256",
+                attribution="Service public de Wallonie (SPW), UCLouvain, ULB, ISSeP",
+                accuracy_label="Officiele globale nauwkeurigheid 91,5%",
+                raw_class_crosswalk=WalousLandCoverService.WALOUS_2018_CLASS_CROSSWALK,
+                comparability_note=(
+                    "De 2018-editie gebruikt een eerdere, deels handmatig geconsolideerde methode. GeoIntel past de "
+                    "officiele 'Classe vue'-crosswalk toe en groepeert de 2018-only serreklasse bij constructies; "
+                    "trends blijven methodologisch begrensde schattingen."
+                ),
+                observation_start=datetime(2018, 1, 1, tzinfo=UTC),
+                observation_end=datetime(2018, 12, 31, 23, 59, 59, tzinfo=UTC),
+            ),
             WalousProduct(
                 key="walous_land_cover_2020",
                 display_name="WALOUS landbedekking 2020",
@@ -104,7 +173,10 @@ class WalousLandCoverService:
                     "47b348f1-6e7a-4baa-963c-0232a43c0cff/WAL_OCS_IA__2020_GEOTIFF_3812.zip"
                 ),
                 source_sha256_filename="walous_land_cover_2020_3812.sha256",
+                attribution=WalousLandCoverService.ATTRIBUTION,
                 accuracy_label="Officiele globale nauwkeurigheid 83,30%",
+                raw_class_crosswalk=None,
+                comparability_note="",
                 observation_start=datetime(2020, 4, 1, tzinfo=UTC),
                 observation_end=datetime(2020, 4, 24, 23, 59, 59, tzinfo=UTC),
             ),
@@ -120,7 +192,10 @@ class WalousLandCoverService:
                     "4e780ba1-463c-478e-95df-d2f1963a150d/WAL_OCS_IA__2023_GEOTIFF_3812.zip"
                 ),
                 source_sha256_filename="walous_land_cover_2023_3812.sha256",
+                attribution=WalousLandCoverService.ATTRIBUTION,
                 accuracy_label="Officiele globale nauwkeurigheid 87,10%",
+                raw_class_crosswalk=None,
+                comparability_note="",
                 observation_start=datetime(2023, 5, 27, tzinfo=UTC),
                 observation_end=datetime(2023, 6, 25, 23, 59, 59, tzinfo=UTC),
             ),
@@ -136,7 +211,10 @@ class WalousLandCoverService:
         resolved = settings or get_settings()
         result: list[dict[str, Any]] = []
         for product in WalousLandCoverService._products().values():
-            configured = resolved.walous_enabled and WalousLandCoverService._source_path(resolved, product).is_file()
+            configured = (
+                resolved.walous_enabled
+                and WalousLandCoverService._source_path(resolved, product).is_file()
+            )
             result.append(
                 ThematicRasterProductRead(
                     key=product.key,
@@ -151,12 +229,20 @@ class WalousLandCoverService:
                     observation_year=product.observation_year,
                     source_version=product.source_version,
                     catalog_url=product.catalog_url,
-                    attribution=WalousLandCoverService.ATTRIBUTION,
+                    attribution=product.attribution,
                     license_note=WalousLandCoverService.LICENSE_NOTE,
                     legend_min_label="WALOUS klasse 1 (kunstmatige bodem)",
                     legend_max_label="WALOUS klasse 90 (loofbomen tot 3 m)",
                     included_source_values=list(WalousLandCoverService.CLASS_LABELS),
-                    limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
+                    limitation_message=" ".join(
+                        part
+                        for part in (
+                            WalousLandCoverService.LIMITATION,
+                            f"{product.accuracy_label}.",
+                            product.comparability_note,
+                        )
+                        if part
+                    ),
                     coverage_zones=["wallonia"],
                     configured=configured,
                     status="configured" if configured else "source_not_provisioned",
@@ -179,21 +265,46 @@ class WalousLandCoverService:
     @staticmethod
     def _scope_geometry(db, project_id: UUID, payload: ThematicRasterAcquireRequest):
         if not db.get(Project, project_id):
-            raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
+            raise AppError(
+                code="PROJECT_NOT_FOUND", message="Project not found", status_code=404
+            )
         if payload.bbox.crs.upper() != "EPSG:4326":
-            raise AppError(code="INVALID_BBOX_CRS", message="WALOUS acquisition requires EPSG:4326", status_code=400)
-        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 values[0] >= values[2] or values[1] >= values[3]:
-            raise AppError(code="INVALID_BBOX", message="WALOUS selection must be a finite non-empty rectangle", status_code=400)
+            raise AppError(
+                code="INVALID_BBOX_CRS",
+                message="WALOUS acquisition requires EPSG:4326",
+                status_code=400,
+            )
+        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 values[0] >= values[2]
+            or values[1] >= values[3]
+        ):
+            raise AppError(
+                code="INVALID_BBOX",
+                message="WALOUS selection must be a finite non-empty rectangle",
+                status_code=400,
+            )
         selection = box(*values)
         if payload.area_id is None:
             return selection, values
         area = db.get(Area, payload.area_id)
         if area is None or area.project_id != project_id:
-            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+            raise AppError(
+                code="AREA_NOT_FOUND", message="Area not found", status_code=404
+            )
         selection = selection.intersection(to_shape(area.geometry))
         if selection.is_empty or selection.area <= 0:
-            raise AppError(code="WALOUS_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+            raise AppError(
+                code="WALOUS_SELECTION_OUTSIDE_AREA",
+                message="Selection does not overlap the selected work area",
+                status_code=422,
+            )
         return selection, values
 
     @staticmethod
@@ -201,6 +312,7 @@ class WalousLandCoverService:
         source_path: Path,
         scope_4326,
         settings: Settings,
+        product: WalousProduct,
     ) -> tuple[bytes, dict[str, Any]]:
         try:
             import numpy as np
@@ -211,20 +323,45 @@ class WalousLandCoverService:
             from rasterio.transform import from_bounds
             from rasterio.windows import from_bounds as window_from_bounds
         except ImportError as exc:
-            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for WALOUS", status_code=503) from exc
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio and numpy are required for WALOUS",
+                status_code=503,
+            ) from exc
 
         resolution = float(settings.walous_analysis_resolution_m)
-        transformer = Transformer.from_crs("EPSG:4326", WalousLandCoverService.SOURCE_CRS, always_xy=True)
+        transformer = Transformer.from_crs(
+            "EPSG:4326", WalousLandCoverService.SOURCE_CRS, always_xy=True
+        )
         scope_metric = shapely_transform(transformer.transform, scope_4326)
         try:
             with rasterio.open(source_path) as source:
-                if source.crs is None or source.crs.to_epsg() != 3812 or source.count != 1:
-                    raise AppError(code="WALOUS_SOURCE_INVALID", message="WALOUS source must be a one-band EPSG:3812 raster", status_code=409)
-                if not all(math.isclose(abs(float(value)), 1.0, abs_tol=0.05) for value in source.res):
-                    raise AppError(code="WALOUS_SOURCE_INVALID", message="WALOUS source must retain the official 1 m resolution", status_code=409)
+                if (
+                    source.crs is None
+                    or source.crs.to_epsg() != 3812
+                    or source.count != 1
+                ):
+                    raise AppError(
+                        code="WALOUS_SOURCE_INVALID",
+                        message="WALOUS source must be a one-band EPSG:3812 raster",
+                        status_code=409,
+                    )
+                if not all(
+                    math.isclose(abs(float(value)), 1.0, abs_tol=0.05)
+                    for value in source.res
+                ):
+                    raise AppError(
+                        code="WALOUS_SOURCE_INVALID",
+                        message="WALOUS source must retain the official 1 m resolution",
+                        status_code=409,
+                    )
                 clipped_geometry = scope_metric.intersection(box(*source.bounds))
                 if clipped_geometry.is_empty or clipped_geometry.area <= 0:
-                    raise AppError(code="WALOUS_SELECTION_OUTSIDE_COVERAGE", message="Selection does not overlap WALOUS coverage", status_code=422)
+                    raise AppError(
+                        code="WALOUS_SELECTION_OUTSIDE_COVERAGE",
+                        message="Selection does not overlap WALOUS coverage",
+                        status_code=422,
+                    )
                 min_x, min_y, max_x, max_y = clipped_geometry.bounds
                 bounds = (
                     math.floor(min_x / resolution) * resolution,
@@ -233,20 +370,45 @@ class WalousLandCoverService:
                     math.ceil(max_y / resolution) * resolution,
                 )
                 width_m, height_m = bounds[2] - bounds[0], bounds[3] - bounds[1]
-                if width_m > settings.walous_max_side_m or height_m > settings.walous_max_side_m:
+                if (
+                    width_m > settings.walous_max_side_m
+                    or height_m > settings.walous_max_side_m
+                ):
                     raise AppError(
                         code="WALOUS_SELECTION_TOO_LARGE",
                         message=f"Select no more than {settings.walous_max_side_m:g} by {settings.walous_max_side_m:g} metres",
                         details={"width_m": width_m, "height_m": height_m},
                         status_code=422,
                     )
-                width, height = max(1, round(width_m / resolution)), max(1, round(height_m / resolution))
+                width, height = (
+                    max(1, round(width_m / resolution)),
+                    max(1, round(height_m / resolution)),
+                )
                 if width * height > settings.walous_max_pixels:
-                    raise AppError(code="WALOUS_SELECTION_TOO_LARGE", message="WALOUS selection exceeds the configured cell limit", details={"pixel_count": width * height, "max_pixels": settings.walous_max_pixels}, status_code=422)
+                    raise AppError(
+                        code="WALOUS_SELECTION_TOO_LARGE",
+                        message="WALOUS selection exceeds the configured cell limit",
+                        details={
+                            "pixel_count": width * height,
+                            "max_pixels": settings.walous_max_pixels,
+                        },
+                        status_code=422,
+                    )
                 window = window_from_bounds(*bounds, transform=source.transform)
-                band = source.read(1, window=window, out_shape=(height, width), masked=True, resampling=Resampling.nearest)
+                band = source.read(
+                    1,
+                    window=window,
+                    out_shape=(height, width),
+                    masked=True,
+                    resampling=Resampling.nearest,
+                )
                 output_transform = from_bounds(*bounds, width, height)
-                outside_scope = geometry_mask([mapping(clipped_geometry)], out_shape=(height, width), transform=output_transform, invert=False)
+                outside_scope = geometry_mask(
+                    [mapping(clipped_geometry)],
+                    out_shape=(height, width),
+                    transform=output_transform,
+                    invert=False,
+                )
                 # The official 2023 GeoTIFF is signed int8 while GDAL exposes
                 # its nodata sentinel as 255. Filling before widening would
                 # therefore reject the sentinel as out of range for int8.
@@ -255,11 +417,18 @@ class WalousLandCoverService:
                 if source.nodata is not None:
                     invalid |= np.isclose(raw.astype("float64"), float(source.nodata))
                 raw[invalid] = WalousLandCoverService.NODATA
-                valid = raw[raw != WalousLandCoverService.NODATA]
-                if valid.size == 0:
-                    raise AppError(code="WALOUS_NO_VALID_DATA", message="WALOUS contains no valid cells in this selection", status_code=422)
-                classes = set(np.unique(valid).astype(int).tolist())
-                unexpected = sorted(classes - set(WalousLandCoverService.CLASS_LABELS))
+                source_valid = raw[raw != WalousLandCoverService.NODATA]
+                if source_valid.size == 0:
+                    raise AppError(
+                        code="WALOUS_NO_VALID_DATA",
+                        message="WALOUS contains no valid cells in this selection",
+                        status_code=422,
+                    )
+                source_classes = set(np.unique(source_valid).astype(int).tolist())
+                governed_source_classes = set(
+                    product.raw_class_crosswalk or WalousLandCoverService.CLASS_LABELS
+                )
+                unexpected = sorted(source_classes - governed_source_classes)
                 if unexpected:
                     raise AppError(
                         code="WALOUS_SOURCE_INVALID_VALUES",
@@ -267,6 +436,18 @@ class WalousLandCoverService:
                         details={"unexpected_classes": unexpected},
                         status_code=409,
                     )
+                if product.raw_class_crosswalk:
+                    normalized = np.full(
+                        raw.shape, WalousLandCoverService.NODATA, dtype="uint8"
+                    )
+                    for (
+                        source_value,
+                        normalized_value,
+                    ) in product.raw_class_crosswalk.items():
+                        normalized[(raw == source_value) & ~invalid] = normalized_value
+                    raw = normalized
+                valid = raw[raw != WalousLandCoverService.NODATA]
+                classes = set(np.unique(valid).astype(int).tolist())
                 profile = {
                     "driver": "GTiff",
                     "width": width,
@@ -288,50 +469,87 @@ class WalousLandCoverService:
                     "height": height,
                     "valid_pixel_count": int(valid.size),
                     "classes_present": sorted(classes),
+                    "source_classes_present": sorted(source_classes),
+                    "class_crosswalk": product.raw_class_crosswalk,
                     "bbox_epsg3812": list(bounds),
                     "source_width": int(source.width),
                     "source_height": int(source.height),
-                    "source_nodata": None if source.nodata is None else float(source.nodata),
+                    "source_nodata": None
+                    if source.nodata is None
+                    else float(source.nodata),
                     "source_resolution_m": 1.0,
                     "analysis_resolution_m": resolution,
                 }
         except AppError:
             raise
         except Exception as exc:
-            raise AppError(code="WALOUS_SOURCE_READ_FAILED", message="The provisioned WALOUS source could not be read", details={"reason": str(exc)}, status_code=500) from exc
+            raise AppError(
+                code="WALOUS_SOURCE_READ_FAILED",
+                message="The provisioned WALOUS source could not be read",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) 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 == WalousLandCoverService.PROVIDER, Dataset.status == "ready")
+            .filter(
+                Dataset.project_id == project_id,
+                Dataset.name == filename,
+                Dataset.source_name == WalousLandCoverService.PROVIDER,
+                Dataset.status == "ready",
+            )
             .order_by(Dataset.imported_at.desc())
             .first()
         )
-        return candidate if candidate and candidate.storage_path and Path(candidate.storage_path).is_file() else None
+        return (
+            candidate
+            if candidate
+            and candidate.storage_path
+            and Path(candidate.storage_path).is_file()
+            else None
+        )
 
     @staticmethod
-    def acquire(db, project_id: UUID, payload: ThematicRasterAcquireRequest, *, settings: Settings | None = None) -> dict[str, Any]:
+    def acquire(
+        db,
+        project_id: UUID,
+        payload: ThematicRasterAcquireRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict[str, Any]:
         resolved = settings or get_settings()
         if not resolved.walous_enabled:
-            raise AppError(code="WALOUS_NOT_CONFIGURED", message="WALOUS bounded analysis is disabled", status_code=503)
+            raise AppError(
+                code="WALOUS_NOT_CONFIGURED",
+                message="WALOUS bounded analysis is disabled",
+                status_code=503,
+            )
         product = WalousLandCoverService._product(payload.product_key)
         source_path = WalousLandCoverService._source_path(resolved, product)
         if not source_path.is_file():
             raise AppError(
                 code="WALOUS_SOURCE_NOT_PROVISIONED",
                 message="The official WALOUS source archive has not been provisioned on this runtime",
-                details={"expected_path": str(source_path), "operator_command": "python scripts/provision_walous_sources.py --years 2020 2023"},
+                details={
+                    "expected_path": str(source_path),
+                    "operator_command": "python scripts/provision_walous_sources.py --years 2018 2020 2023",
+                },
                 status_code=503,
             )
-        scope, bbox_4326 = WalousLandCoverService._scope_geometry(db, project_id, payload)
+        scope, bbox_4326 = WalousLandCoverService._scope_geometry(
+            db, project_id, payload
+        )
         identity = {
             "product_key": product.key,
             "bbox_epsg4326": [round(float(value), 8) for value in bbox_4326],
             "area_id": str(payload.area_id) if payload.area_id else None,
             "analysis_resolution_m": resolved.walous_analysis_resolution_m,
         }
-        request_hash = hashlib.sha256(json.dumps(identity, sort_keys=True).encode()).hexdigest()
+        request_hash = hashlib.sha256(
+            json.dumps(identity, sort_keys=True).encode()
+        ).hexdigest()
         filename = f"walous_{product.observation_year}_{request_hash[:12]}_3812.tif"
         if not payload.force_refresh:
             cached = WalousLandCoverService._cached_dataset(db, project_id, filename)
@@ -345,7 +563,12 @@ class WalousLandCoverService:
                     display_name=product.display_name,
                     theme=WalousLandCoverService.THEME,
                     metric_kind=WalousLandCoverService.METRIC_KIND,
-                    resolution_m=float(metadata.get("analysis_resolution_m", resolved.walous_analysis_resolution_m)),
+                    resolution_m=float(
+                        metadata.get(
+                            "analysis_resolution_m",
+                            resolved.walous_analysis_resolution_m,
+                        )
+                    ),
                     width=int((cached.metadata_json or {}).get("width", 0)),
                     height=int((cached.metadata_json or {}).get("height", 0)),
                     valid_pixel_count=int(metadata.get("valid_pixel_count", 0)),
@@ -353,16 +576,39 @@ class WalousLandCoverService:
                     bbox_epsg3812=list(metadata.get("bbox_epsg3812") or []),
                     observation_year=product.observation_year,
                     source_value_unit=WalousLandCoverService.SOURCE_VALUE_UNIT,
-                    attribution=WalousLandCoverService.ATTRIBUTION,
-                    limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
+                    attribution=product.attribution,
+                    limitation_message=" ".join(
+                        part
+                        for part in (
+                            WalousLandCoverService.LIMITATION,
+                            f"{product.accuracy_label}.",
+                            product.comparability_note,
+                        )
+                        if part
+                    ),
                 ).model_dump(mode="json")
 
-        content, validation = WalousLandCoverService._read_source_window(source_path, scope, resolved)
+        content, validation = WalousLandCoverService._read_source_window(
+            source_path, scope, resolved, product
+        )
         source_sha256_path = source_path.with_name(product.source_sha256_filename)
-        source_sha256 = source_sha256_path.read_text(encoding="ascii").strip().split()[0] if source_sha256_path.is_file() else None
+        source_sha256 = (
+            source_sha256_path.read_text(encoding="ascii").strip().split()[0]
+            if source_sha256_path.is_file()
+            else None
+        )
         acquired_at = datetime.now(UTC)
         observed_at = product.observation_end
-        spatial_series_hash = hashlib.sha256(json.dumps({"bbox": identity["bbox_epsg4326"], "area_id": identity["area_id"], "resolution": identity["analysis_resolution_m"]}, sort_keys=True).encode()).hexdigest()[:24]
+        spatial_series_hash = hashlib.sha256(
+            json.dumps(
+                {
+                    "bbox": identity["bbox_epsg4326"],
+                    "area_id": identity["area_id"],
+                    "resolution": identity["analysis_resolution_m"],
+                },
+                sort_keys=True,
+            ).encode()
+        ).hexdigest()[:24]
         dataset = DatasetService.import_raster_bytes(
             db,
             project_id=project_id,
@@ -394,14 +640,24 @@ class WalousLandCoverService:
                 "observation_end": product.observation_end.isoformat(),
                 "valid_pixel_count": validation["valid_pixel_count"],
                 "classes_present": validation["classes_present"],
+                "source_classes_present": validation["source_classes_present"],
+                "class_crosswalk": validation["class_crosswalk"],
                 "bbox_epsg4326": bbox_4326,
                 "bbox_epsg3812": validation["bbox_epsg3812"],
                 "coverage_zones": ["wallonia"],
                 "catalog_url": product.catalog_url,
                 "download_url": product.download_url,
-                "attribution": WalousLandCoverService.ATTRIBUTION,
+                "attribution": product.attribution,
                 "license_note": WalousLandCoverService.LICENSE_NOTE,
-                "limitation_message": f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
+                "limitation_message": " ".join(
+                    part
+                    for part in (
+                        WalousLandCoverService.LIMITATION,
+                        f"{product.accuracy_label}.",
+                        product.comparability_note,
+                    )
+                    if part
+                ),
             },
             provenance_metadata={
                 "acquisition": "operator_provisioned_official_archive_bounded_window",
@@ -430,77 +686,169 @@ class WalousLandCoverService:
             bbox_epsg3812=validation["bbox_epsg3812"],
             observation_year=product.observation_year,
             source_value_unit=WalousLandCoverService.SOURCE_VALUE_UNIT,
-            attribution=WalousLandCoverService.ATTRIBUTION,
-            limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
+            attribution=product.attribution,
+            limitation_message=" ".join(
+                part
+                for part in (
+                    WalousLandCoverService.LIMITATION,
+                    f"{product.accuracy_label}.",
+                    product.comparability_note,
+                )
+                if part
+            ),
         ).model_dump(mode="json")
 
     @staticmethod
-    def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> tuple[Dataset, WalousProduct]:
+    def _load_dataset(
+        db, project_id: UUID, dataset_id: UUID
+    ) -> tuple[Dataset, WalousProduct]:
         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 != WalousLandCoverService.PROVIDER:
-            raise AppError(code="INVALID_WALOUS_DATASET", message="WALOUS analysis requires a governed WALOUS 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 WALOUS raster is unavailable", status_code=404)
-        product = WalousLandCoverService._product(str((dataset.source_metadata or {}).get("product_key") or ""))
+            raise AppError(
+                code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404
+            )
+        if (
+            dataset.dataset_type != "raster"
+            or dataset.source_name != WalousLandCoverService.PROVIDER
+        ):
+            raise AppError(
+                code="INVALID_WALOUS_DATASET",
+                message="WALOUS analysis requires a governed WALOUS 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 WALOUS raster is unavailable",
+                status_code=404,
+            )
+        product = WalousLandCoverService._product(
+            str((dataset.source_metadata or {}).get("product_key") or "")
+        )
         return dataset, product
 
     @staticmethod
-    def _analysis_geometry(db, project_id: UUID, payload: ThematicRasterSelectionRequest):
-        selection = box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+    def _analysis_geometry(
+        db, project_id: UUID, payload: ThematicRasterSelectionRequest
+    ):
+        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 area is None or area.project_id != project_id:
-            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+            raise AppError(
+                code="AREA_NOT_FOUND", message="Area not found", status_code=404
+            )
         selection = selection.intersection(to_shape(area.geometry))
         if selection.is_empty or selection.area <= 0:
-            raise AppError(code="WALOUS_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+            raise AppError(
+                code="WALOUS_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: ThematicRasterSelectionRequest) -> dict[str, Any]:
-        dataset, product = WalousLandCoverService._load_dataset(db, project_id, dataset_id)
-        selection_4326 = WalousLandCoverService._analysis_geometry(db, project_id, payload)
+    def analyze(
+        db, project_id: UUID, dataset_id: UUID, payload: ThematicRasterSelectionRequest
+    ) -> dict[str, Any]:
+        dataset, product = WalousLandCoverService._load_dataset(
+            db, project_id, dataset_id
+        )
+        selection_4326 = WalousLandCoverService._analysis_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 WALOUS analysis", status_code=503) from exc
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio and numpy are required for WALOUS analysis",
+                status_code=503,
+            ) from exc
         try:
             with rasterio.open(dataset.storage_path) as source:
-                transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True)
-                selection_metric = shapely_transform(transformer.transform, selection_4326)
+                transformer = Transformer.from_crs(
+                    "EPSG:4326", source.crs, always_xy=True
+                )
+                selection_metric = shapely_transform(
+                    transformer.transform, selection_4326
+                )
                 geometry = selection_metric.intersection(box(*source.bounds))
                 if geometry.is_empty or geometry.area <= 0:
-                    raise AppError(code="WALOUS_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted WALOUS raster", status_code=422)
-                clipped, transform = mask(source, [mapping(geometry)], crop=True, filled=False, indexes=[1])
+                    raise AppError(
+                        code="WALOUS_SELECTION_OUTSIDE_DATASET",
+                        message="Selection does not overlap the persisted WALOUS raster",
+                        status_code=422,
+                    )
+                clipped, transform = mask(
+                    source, [mapping(geometry)], crop=True, filled=False, indexes=[1]
+                )
                 band = np.ma.asarray(clipped[0])
                 raw = np.asarray(np.ma.getdata(band), dtype="uint8")
-                selected = geometry_mask([mapping(geometry)], out_shape=raw.shape, transform=transform, invert=True)
-                valid = selected & ~np.ma.getmaskarray(band) & (raw != WalousLandCoverService.NODATA)
+                selected = geometry_mask(
+                    [mapping(geometry)],
+                    out_shape=raw.shape,
+                    transform=transform,
+                    invert=True,
+                )
+                valid = (
+                    selected
+                    & ~np.ma.getmaskarray(band)
+                    & (raw != WalousLandCoverService.NODATA)
+                )
                 values = raw[valid]
                 selected_count = int(selected.sum())
                 valid_count = int(values.size)
                 if not valid_count:
-                    raise AppError(code="WALOUS_NO_VALID_DATA", message="WALOUS contains no valid cells in this selection", status_code=422)
+                    raise AppError(
+                        code="WALOUS_NO_VALID_DATA",
+                        message="WALOUS contains no valid cells in this selection",
+                        status_code=422,
+                    )
                 cell_area_m2 = abs(float(source.res[0]) * float(source.res[1]))
         except AppError:
             raise
         except Exception as exc:
-            raise AppError(code="WALOUS_ANALYSIS_FAILED", message="The persisted WALOUS raster could not be analysed", details={"reason": str(exc)}, status_code=500) from exc
+            raise AppError(
+                code="WALOUS_ANALYSIS_FAILED",
+                message="The persisted WALOUS raster could not be analysed",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc
 
         def area_for(classes: set[int]) -> float:
-            return float(np.count_nonzero(np.isin(values, list(classes))) * cell_area_m2 / 10_000.0)
+            return float(
+                np.count_nonzero(np.isin(values, list(classes)))
+                * cell_area_m2
+                / 10_000.0
+            )
 
         metric_specs = [
-            ("land_cover_observed_area_ha", "Gekarteerde landbedekking", set(WalousLandCoverService.CLASS_LABELS)),
+            (
+                "land_cover_observed_area_ha",
+                "Gekarteerde landbedekking",
+                set(WalousLandCoverService.CLASS_LABELS),
+            ),
             ("forest_cover_area_ha", "Boom- en bosbedekking", {8, 9, 80, 90}),
             ("surface_water_area_ha", "Oppervlaktewater", {5}),
-            ("artificial_cover_area_ha", "Kunstmatige bedekking en constructies", {1, 2, 3}),
+            (
+                "artificial_cover_area_ha",
+                "Kunstmatige bedekking en constructies",
+                {1, 2, 3},
+            ),
             ("annual_herbaceous_cover_area_ha", "Jaarlijks wisselende kruidlaag", {6}),
             ("permanent_herbaceous_cover_area_ha", "Jaarronde kruidlaag", {7}),
             ("bare_soil_area_ha", "Kale bodem", {4}),
@@ -537,25 +885,52 @@ class WalousLandCoverService:
                 primary_metric_key=primary.metric_key,
                 metrics=metrics,
             ),
-            unsupported_metrics=["legal_land_use", "ownership", "tree_count", "timber_volume", "water_volume"],
-            limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
+            unsupported_metrics=[
+                "legal_land_use",
+                "ownership",
+                "tree_count",
+                "timber_volume",
+                "water_volume",
+            ],
+            limitation_message=" ".join(
+                part
+                for part in (
+                    WalousLandCoverService.LIMITATION,
+                    f"{product.accuracy_label}.",
+                    product.comparability_note,
+                )
+                if part
+            ),
             generated_at=datetime.now(UTC).isoformat(),
         ).model_dump(mode="json")
 
     @staticmethod
-    def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
-        dataset, _product = WalousLandCoverService._load_dataset(db, project_id, dataset_id)
+    def render_png(
+        db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800
+    ) -> bytes:
+        dataset, _product = WalousLandCoverService._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 WALOUS rendering", status_code=503) from exc
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio, numpy and Pillow are required for WALOUS rendering",
+                status_code=503,
+            ) from exc
         with rasterio.open(dataset.storage_path) as source:
             scale = min(1.0, max_dimension / max(source.width, source.height))
-            width, height = max(1, round(source.width * scale)), max(1, round(source.height * scale))
-            values = source.read(1, out_shape=(height, width), masked=True, resampling=Resampling.nearest)
+            width, height = (
+                max(1, round(source.width * scale)),
+                max(1, round(source.height * scale)),
+            )
+            values = source.read(
+                1, out_shape=(height, width), masked=True, resampling=Resampling.nearest
+            )
             raw = np.asarray(np.ma.getdata(values), dtype="uint8")
             rgba = np.zeros((height, width, 4), dtype="uint8")
             for value, color in WalousLandCoverService.CLASS_COLORS.items():