feat: add governed hydrology and historical imagery

CHANGELOG.md

@@ -7,6 +7,23 @@
 
 # Changelog
 
+## Sprint 203 Governed hydrology and historical imagery (2026-07-15)
+
+- Added an explicit Waterinfo KiWIS operator for annual water-level and
+  discharge station histories. It filters against the persisted Area, retains
+  raw JSON/checksums and imports each real station/year through DatasetService.
+- Added numeric `mean` selection aggregation for station measurements while
+  keeping every station in an independent temporal series with an explicit
+  point-versus-area/volume limitation.
+- Added a fixed official orthophoto product registry covering current, annual
+  2012-2025, older winter periods, RGB 1979-1990 and panchromatic 1971.
+- Added historical raster temporal provenance, a constrained browser PNG
+  endpoint and a MapLibre image overlay/product selector.
+- Kept configured-YOLO/current-GRB QA exclusive to the most-recent product;
+  historical imagery never produces fake current-state quality metrics.
+- Converted BWK/Natura 2000, agricultural parcels, Buildings Register, DHMV
+  and bathymetry into an ordered acceptance-criteria backlog.
+
 ## Sprint 202 Source intelligence, full evolution metrics and local assistant (2026-07-15)
 
 - Extended temporal comparisons with exact persisted-Area filtering, every

backend/README.md

@@ -1121,16 +1121,40 @@ is never presented as a complete result.
 
 ## Bounded official orthophoto acquisition
 
-`POST /api/v1/projects/{project_id}/datasets/orthophoto/acquire` accepts an
-explicit EPSG:4326 map rectangle and stores the official Digitaal Vlaanderen
-`OMWRGBMRVL`/`Ortho` response as a canonical EPSG:31370 raster Dataset. The
+`GET /api/v1/projects/{project_id}/datasets/orthophoto/products` lists the
+governed product allowlist. `POST .../datasets/orthophoto/acquire` accepts an
+explicit EPSG:4326 map rectangle plus `product_key` and stores the official
+Digitaal Vlaanderen WMS response as a canonical EPSG:31370 raster Dataset. The
 default safety envelope is 128-1,024 m per side, 1 m/pixel, 32 MiB and a
 24-hour exact-request cache. It runs synchronously behind the existing Job
 abstraction and never during startup.
 
+Available products cover the most recent winter image, annual winter mosaics
+for 2012-2025, three older winter periods, RGB 1979-1990 and panchromatic 1971.
+Historical products persist validity metadata and are deliberately excluded
+from configured-YOLO/current-GRB QA. `GET .../datasets/{dataset_id}/raster/image`
+is the constrained binary PNG endpoint used by the MapLibre image overlay.
+
 Settings: `ORTHOPHOTO_ENABLED`, `ORTHOPHOTO_WMS_URL`,
 `ORTHOPHOTO_WMS_LAYER`, `ORTHOPHOTO_RESOLUTION_M`,
 `ORTHOPHOTO_MIN_SIDE_M`, `ORTHOPHOTO_MAX_SIDE_M`,
 `ORTHOPHOTO_TIMEOUT_SECONDS`, `ORTHOPHOTO_MAX_RESPONSE_MB` and
 `ORTHOPHOTO_CACHE_TTL_HOURS`. Keep the official HTTPS URL and 1 m profile
 unless a separately verified deployment/model profile requires a change.
+
+## Waterinfo station histories
+
+Run the explicit operator after the regional workspace and Mol Area exist:
+
+```bash
+docker exec geointel python /app/scripts/provision_waterinfo_station_history.py \
+  --project-name "Kempen Regional Workbench" \
+  --area-name "Gemeente Mol" \
+  --from-year 2013 --to-year 2025
+```
+
+The command retains raw KiWIS JSON/checksums and imports only real annual
+observations through the canonical dataset upload API. Every station has its
+own temporal-series key. Water levels and discharges remain Point measurements;
+they are never averaged across stations or presented as municipal water volume.
+Use `--fetch-only` to prepare and audit artifacts without persistence.

backend/app/api/routes/datasets.py

@@ -6,10 +6,10 @@ from typing import Any
 from uuid import UUID
 from uuid import UUID as _UUID
 
-from fastapi import APIRouter, Depends, File, Form, HTTPException, Query
+from fastapi import APIRouter, Depends, File, Form, HTTPException, Query, Response
 from fastapi import UploadFile
 from sqlalchemy.orm import Session
-from app.models import Area
+from app.models import Area, Project
 
 from app.core.errors import AppError
 from app.db.session import get_db
@@ -143,6 +143,14 @@ def acquire_bounded_orthophoto(
     return envelope(job)
 
 
+@router.get("/datasets/orthophoto/products", response_model=dict)
+def list_orthophoto_products(project_id: UUID, db: Session = Depends(get_db)):
+    if not db.get(Project, project_id):
+        raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
+    items = OrthophotoAcquisitionService.list_products()
+    return envelope({"items": items, "total": len(items)})
+
+
 @router.get("/datasets", response_model=dict)
 def list_datasets(
     project_id: UUID,
@@ -426,6 +434,20 @@ def raster_preview_readiness(
     return envelope(RasterOperationsService.preview(db, dataset_id))
 
 
+@router.get("/datasets/{dataset_id}/raster/image")
+def raster_orthophoto_image(
+    project_id: UUID,
+    dataset_id: UUID,
+    db: Session = Depends(get_db),
+):
+    content = OrthophotoAcquisitionService.render_png(db, project_id, dataset_id)
+    return Response(
+        content=content,
+        media_type="image/png",
+        headers={"Cache-Control": "private, max-age=86400"},
+    )
+
+
 @router.get("/datasets/{dataset_id}/raster/stats", response_model=dict)
 def raster_stats(
     project_id: UUID,

backend/app/schemas/__init__.py

@@ -32,7 +32,7 @@ from .segmentation import (
 )
 from .health import HealthResponse, SystemCapabilities
 from .job import JobCreate, JobList, JobRead, JobStatus
-from .orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult
+from .orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult, OrthophotoProductRead
 from .external import (
     ExternalFetchRequest,
     ExternalFetchResponse,
@@ -134,6 +134,7 @@ __all__ = [
     "JobStatus",
     "OrthophotoAcquireRequest",
     "OrthophotoAcquisitionResult",
+    "OrthophotoProductRead",
     "VectorBBoxResponse",
     "VectorClipRequest",
     "VectorBufferRequest",

backend/app/schemas/orthophoto.py

@@ -10,13 +10,31 @@ from .operations import VectorSelectionBBox
 class OrthophotoAcquireRequest(BaseModel):
     bbox: VectorSelectionBBox
     area_id: UUID | None = None
+    product_key: str = "most_recent"
     force_refresh: bool = False
 
 
+class OrthophotoProductRead(BaseModel):
+    key: str
+    display_name: str
+    observation_label: str
+    temporal_granularity: str
+    native_resolution_m: float
+    supports_detection: bool
+    color_mode: str
+    catalog_url: str
+    limitation_message: str
+
+
 class OrthophotoAcquisitionResult(BaseModel):
     output_dataset_id: UUID
     reused: bool
     provider: str
+    product_key: str
+    display_name: str
+    observation_label: str
+    temporal_granularity: str
+    supports_detection: bool
     layer: str
     width: int
     height: int

backend/app/services/dataset_service.py

@@ -422,6 +422,11 @@ class DatasetService:
         source_metadata: dict[str, Any],
         provenance_metadata: dict[str, Any],
         area_id: UUID | None = None,
+        temporal_series_key: str | None = None,
+        observed_at: datetime | None = None,
+        valid_from: datetime | None = None,
+        valid_to: datetime | None = None,
+        temporal_granularity: str | None = None,
         source_version: str | None = None,
         content_type: str = "image/tiff",
     ) -> DatasetCreateResponse:
@@ -438,6 +443,14 @@ class DatasetService:
         safe_filename = DatasetService._validate_upload_filename(filename)
         if DatasetService._extension_for_path(safe_filename) not in DatasetService.RASTER_EXTENSIONS:
             raise AppError(code="INVALID_UPLOAD", message="Raster artifacts require a GeoTIFF filename", status_code=415)
+        temporal = DatasetService._validate_temporal_metadata(
+            temporal_series_key=temporal_series_key,
+            observed_at=observed_at,
+            valid_from=valid_from,
+            valid_to=valid_to,
+            temporal_granularity=temporal_granularity,
+            source_version=source_version,
+        )
 
         dataset_id = uuid.uuid4()
         storage_info = StorageService.persist_dataset_file(
@@ -462,7 +475,7 @@ class DatasetService:
                 source_metadata=source_metadata,
                 provenance_metadata=provenance_metadata,
                 imported_at=datetime.now(timezone.utc),
-                source_version=source_version,
+                **temporal,
                 storage_path=storage_info["storage_path"],
                 original_filename=storage_info["original_filename"],
                 stored_filename=storage_info["stored_filename"],
@@ -483,6 +496,9 @@ class DatasetService:
                     version=1,
                     storage_path=dataset.storage_path,
                     source_version=dataset.source_version,
+                    observed_at=dataset.observed_at,
+                    valid_from=dataset.valid_from,
+                    valid_to=dataset.valid_to,
                     checksum_sha256=dataset.checksum_sha256,
                     source_metadata=dataset.source_metadata,
                     provenance_metadata=dataset.provenance_metadata,

backend/app/services/geo_assistant_service.py

@@ -349,6 +349,7 @@ class GeoAssistantService:
                             "measurement_quality": (
                                 "schatting" if summary["is_estimate"] else "exact_binnen_bronrepresentatie"
                             ),
+                            "warning": summary.get("warning"),
                         }
                     )
                     if dataset.id not in source_dataset_ids:

backend/app/services/orthophoto_acquisition_service.py

@@ -1,9 +1,11 @@
 from __future__ import annotations
 
 import hashlib
+import io
 import json
 import math
 import warnings
+from dataclasses import dataclass
 from datetime import UTC, datetime, timedelta
 from pathlib import Path
 from typing import Any, Callable
@@ -20,21 +22,169 @@ 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.orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult
+from app.schemas.orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult, OrthophotoProductRead
 from app.services.dataset_service import DatasetService
 
 
+@dataclass(frozen=True)
+class OrthophotoProduct:
+    key: str
+    display_name: str
+    observation_label: str
+    temporal_granularity: str
+    native_resolution_m: float
+    wms_url: str
+    layer: str
+    catalog_url: str
+    limitation_message: str
+    supports_detection: bool = False
+    color_mode: str = "rgb"
+    observed_at: datetime | None = None
+    valid_from: datetime | None = None
+    valid_to: datetime | None = None
+
+
 class OrthophotoAcquisitionService:
     PROVIDER = "digitaal_vlaanderen_orthophoto"
     ATTRIBUTION = "Bron: Orthofotomozaiek Vlaanderen, Digitaal Vlaanderen"
     CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/orthofotomozaiek-middenschalig-winteropnamen-kleur-meest-recent-vlaanderen"
     LIMITATION = "Meest recente samengestelde winterorthofoto op het moment van de aanvraag; geen historische opnamedatum per pixel."
+    HISTORICAL_WINTER_WMS_URL = "https://geo.api.vlaanderen.be/OMW/wms"
+    HISTORICAL_WINTER_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/wmts-orthofotomozaiek-middenschalig-winteropnamen"
+    HISTORICAL_SUMMER_WMS_URL = "https://geo.api.vlaanderen.be/OKZ/wms"
+    HISTORICAL_SUMMER_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/orthofotomozaiek-kleinschalig-zomeropnamen"
+
+    @staticmethod
+    def _products(settings: Settings) -> dict[str, OrthophotoProduct]:
+        products: list[OrthophotoProduct] = [
+            OrthophotoProduct(
+                key="most_recent",
+                display_name="Meest recente winterluchtbeeld",
+                observation_label="Meest recent beschikbaar",
+                temporal_granularity="snapshot",
+                native_resolution_m=0.15,
+                wms_url=settings.orthophoto_wms_url,
+                layer=settings.orthophoto_wms_layer,
+                catalog_url=OrthophotoAcquisitionService.CATALOG_URL,
+                limitation_message=OrthophotoAcquisitionService.LIMITATION,
+                supports_detection=True,
+            )
+        ]
+        for year in range(2025, 2011, -1):
+            products.append(
+                OrthophotoProduct(
+                    key=str(year),
+                    display_name=f"Winterluchtbeeld {year}",
+                    observation_label=str(year),
+                    temporal_granularity="year",
+                    native_resolution_m=0.15 if year >= 2022 else 0.25,
+                    wms_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_WMS_URL,
+                    layer=f"OMWRGB{year % 100:02d}VL",
+                    catalog_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_CATALOG_URL,
+                    limitation_message=(
+                        "Officiële samengestelde winterorthofoto voor deze jaargang; de exacte opnamedatum kan per tegel verschillen. "
+                        "Historische beelden worden niet met de actuele GRB-toestand gevalideerd."
+                    ),
+                    observed_at=datetime(year, 1, 1, tzinfo=UTC),
+                    valid_from=datetime(year, 1, 1, tzinfo=UTC),
+                    valid_to=datetime(year, 12, 31, 23, 59, 59, tzinfo=UTC),
+                )
+            )
+        for key, start_year, end_year, layer in (
+            ("2008_2011", 2008, 2011, "OMWRGB08_11VL"),
+            ("2005_2007", 2005, 2007, "OMWRGB05_07VL"),
+            ("2000_2003", 2000, 2003, "OMWRGB00_03VL"),
+        ):
+            products.append(
+                OrthophotoProduct(
+                    key=key,
+                    display_name=f"Winterluchtbeeld {start_year}-{end_year}",
+                    observation_label=f"{start_year}-{end_year}",
+                    temporal_granularity="period",
+                    native_resolution_m=0.25,
+                    wms_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_WMS_URL,
+                    layer=layer,
+                    catalog_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_CATALOG_URL,
+                    limitation_message=(
+                        "Officiële samengestelde winterorthofoto uit een meerjarige opnameperiode; dit is geen exacte jaaropname. "
+                        "Historische beelden worden niet met de actuele GRB-toestand gevalideerd."
+                    ),
+                    observed_at=datetime(start_year, 1, 1, tzinfo=UTC),
+                    valid_from=datetime(start_year, 1, 1, tzinfo=UTC),
+                    valid_to=datetime(end_year, 12, 31, 23, 59, 59, tzinfo=UTC),
+                )
+            )
+        products.extend(
+            [
+                OrthophotoProduct(
+                    key="1979_1990",
+                    display_name="Zomerluchtbeeld 1979-1990",
+                    observation_label="1979-1990",
+                    temporal_granularity="period",
+                    native_resolution_m=1.0,
+                    wms_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_WMS_URL,
+                    layer="OKZRGB79_90VL",
+                    catalog_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_CATALOG_URL,
+                    limitation_message="Kleinschalig RGB-mozaïek uit meerdere zomervluchten tussen 1979 en 1990; geen exacte jaartoestand.",
+                    observed_at=datetime(1979, 1, 1, tzinfo=UTC),
+                    valid_from=datetime(1979, 1, 1, tzinfo=UTC),
+                    valid_to=datetime(1990, 12, 31, 23, 59, 59, tzinfo=UTC),
+                ),
+                OrthophotoProduct(
+                    key="1971",
+                    display_name="Zomerluchtbeeld 1971",
+                    observation_label="1971",
+                    temporal_granularity="year",
+                    native_resolution_m=1.0,
+                    wms_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_WMS_URL,
+                    layer="OKZPAN71VL",
+                    catalog_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_CATALOG_URL,
+                    limitation_message="Kleinschalig panchromatisch mozaïek uit 1971; zwart-wit en niet geschikt voor het huidige RGB-detectiemodel.",
+                    color_mode="panchromatic",
+                    observed_at=datetime(1971, 1, 1, tzinfo=UTC),
+                    valid_from=datetime(1971, 1, 1, tzinfo=UTC),
+                    valid_to=datetime(1971, 12, 31, 23, 59, 59, tzinfo=UTC),
+                ),
+            ]
+        )
+        return {product.key: product for product in products}
+
+    @staticmethod
+    def list_products(settings: Settings | None = None) -> list[dict[str, Any]]:
+        resolved_settings = settings or get_settings()
+        return [
+            OrthophotoProductRead(
+                key=product.key,
+                display_name=product.display_name,
+                observation_label=product.observation_label,
+                temporal_granularity=product.temporal_granularity,
+                native_resolution_m=product.native_resolution_m,
+                supports_detection=product.supports_detection,
+                color_mode=product.color_mode,
+                catalog_url=product.catalog_url,
+                limitation_message=product.limitation_message,
+            ).model_dump()
+            for product in OrthophotoAcquisitionService._products(resolved_settings).values()
+        ]
+
+    @staticmethod
+    def _product(product_key: str, settings: Settings) -> OrthophotoProduct:
+        product = OrthophotoAcquisitionService._products(settings).get(product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="ORTHOPHOTO_PRODUCT_NOT_SUPPORTED",
+                message="Select an orthophoto product from the official product registry",
+                details={"product_key": product_key},
+                status_code=422,
+            )
+        return product
 
     @staticmethod
     def _prepared_request(
         payload: OrthophotoAcquireRequest,
         settings: Settings,
     ) -> dict[str, Any]:
+        product = OrthophotoAcquisitionService._product(payload.product_key, settings)
         if payload.bbox.crs.upper() != "EPSG:4326":
             raise AppError(code="INVALID_CRS", message="Orthophoto selection bbox must use EPSG:4326", status_code=400)
         min_x = float(payload.bbox.min_x)
@@ -68,8 +218,9 @@ class OrthophotoAcquisitionService:
         bbox_31370 = [float(value) for value in lambert_bounds]
         request_identity = {
             "provider": OrthophotoAcquisitionService.PROVIDER,
-            "wms_url": settings.orthophoto_wms_url,
-            "layer": settings.orthophoto_wms_layer,
+            "product_key": product.key,
+            "wms_url": product.wms_url,
+            "layer": product.layer,
             "bbox_epsg4326": [round(value, 8) for value in bbox_4326],
             "bbox_epsg31370": [round(value, 3) for value in bbox_31370],
             "width": width,
@@ -77,11 +228,18 @@ class OrthophotoAcquisitionService:
             "resolution_m": settings.orthophoto_resolution_m,
         }
         request_hash = hashlib.sha256(json.dumps(request_identity, sort_keys=True).encode("utf-8")).hexdigest()
+        spatial_identity = {
+            "bbox_epsg4326": request_identity["bbox_epsg4326"],
+            "width": width,
+            "height": height,
+            "resolution_m": settings.orthophoto_resolution_m,
+        }
+        spatial_hash = hashlib.sha256(json.dumps(spatial_identity, sort_keys=True).encode("utf-8")).hexdigest()
         params = {
             "SERVICE": "WMS",
             "VERSION": "1.3.0",
             "REQUEST": "GetMap",
-            "LAYERS": settings.orthophoto_wms_layer,
+            "LAYERS": product.layer,
             "STYLES": "",
             "FORMAT": "image/tiff",
             "CRS": "EPSG:31370",
@@ -91,8 +249,10 @@ class OrthophotoAcquisitionService:
         }
         return {
             **request_identity,
+            "product": product,
+            "spatial_hash": spatial_hash,
             "request_hash": request_hash,
-            "request_url": f"{settings.orthophoto_wms_url}?{urlencode(params)}",
+            "request_url": f"{product.wms_url}?{urlencode(params)}",
             "params": params,
             "bbox_epsg4326": bbox_4326,
             "bbox_epsg31370": bbox_31370,
@@ -124,8 +284,14 @@ class OrthophotoAcquisitionService:
             )
 
     @staticmethod
-    def _cached_dataset(db, project_id: UUID, filename: str, settings: Settings) -> Dataset | None:
-        if settings.orthophoto_cache_ttl_hours <= 0:
+    def _cached_dataset(
+        db,
+        project_id: UUID,
+        filename: str,
+        settings: Settings,
+        product: OrthophotoProduct,
+    ) -> Dataset | None:
+        if product.key == "most_recent" and settings.orthophoto_cache_ttl_hours <= 0:
             return None
         candidate = (
             db.query(Dataset)
@@ -145,7 +311,7 @@ class OrthophotoAcquisitionService:
             return None
         if imported_at.tzinfo is None:
             imported_at = imported_at.replace(tzinfo=UTC)
-        if datetime.now(UTC) - imported_at > timedelta(hours=settings.orthophoto_cache_ttl_hours):
+        if product.key == "most_recent" and datetime.now(UTC) - imported_at > timedelta(hours=settings.orthophoto_cache_ttl_hours):
             return None
         return candidate
 
@@ -196,7 +362,9 @@ class OrthophotoAcquisitionService:
                 with warnings.catch_warnings():
                     warnings.simplefilter("ignore", NotGeoreferencedWarning)
                     with source_memory.open() as source:
-                        if source.width != prepared["width"] or source.height != prepared["height"] or source.count < 3:
+                        product: OrthophotoProduct = prepared["product"]
+                        minimum_band_count = 1 if product.color_mode == "panchromatic" else 3
+                        if source.width != prepared["width"] or source.height != prepared["height"] or source.count < minimum_band_count:
                             raise AppError(
                                 code="ORTHOPHOTO_PROVIDER_INVALID_RESPONSE",
                                 message="Official orthophoto dimensions or RGB bands do not match the bounded request",
@@ -216,7 +384,7 @@ class OrthophotoAcquisitionService:
                         with output_memory.open(**profile) as output:
                             output.write(image)
                             output.update_tags(
-                                source="Digitaal Vlaanderen OMWRGBMRVL WMS Ortho layer",
+                                source=f"Digitaal Vlaanderen WMS {product.layer}",
                                 source_url=prepared["request_url"],
                                 attribution=OrthophotoAcquisitionService.ATTRIBUTION,
                                 acquisition="explicit_bounded_map_selection",
@@ -245,44 +413,74 @@ class OrthophotoAcquisitionService:
         if not resolved_settings.orthophoto_enabled:
             raise AppError(code="ORTHOPHOTO_NOT_CONFIGURED", message="Official orthophoto acquisition is disabled", status_code=503)
         prepared = OrthophotoAcquisitionService._prepared_request(payload, resolved_settings)
+        product: OrthophotoProduct = prepared["product"]
         OrthophotoAcquisitionService._validate_area_scope(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
-        filename = f"orthofoto_selectie_{prepared['request_hash'][:12]}.tif"
+        filename = f"orthofoto_{product.key}_{prepared['request_hash'][:12]}.tif"
 
-        cached = None if payload.force_refresh else OrthophotoAcquisitionService._cached_dataset(db, project_id, filename, resolved_settings)
+        cached = None if payload.force_refresh else OrthophotoAcquisitionService._cached_dataset(
+            db,
+            project_id,
+            filename,
+            resolved_settings,
+            product,
+        )
         if cached is not None:
             return OrthophotoAcquisitionResult(
                 output_dataset_id=cached.id,
                 reused=True,
                 provider=OrthophotoAcquisitionService.PROVIDER,
-                layer=resolved_settings.orthophoto_wms_layer,
+                product_key=product.key,
+                display_name=product.display_name,
+                observation_label=product.observation_label,
+                temporal_granularity=product.temporal_granularity,
+                supports_detection=product.supports_detection,
+                layer=product.layer,
                 width=prepared["width"],
                 height=prepared["height"],
                 resolution_m=resolved_settings.orthophoto_resolution_m,
                 bbox_epsg4326=prepared["bbox_epsg4326"],
                 bbox_epsg31370=prepared["bbox_epsg31370"],
                 attribution=OrthophotoAcquisitionService.ATTRIBUTION,
-                limitation_message=OrthophotoAcquisitionService.LIMITATION,
+                limitation_message=product.limitation_message,
             ).model_dump(mode="json")
 
         raw_content, response_content_type = OrthophotoAcquisitionService._fetch(prepared["request_url"], resolved_settings, opener)
         geotiff_content = OrthophotoAcquisitionService._georeference_tiff(raw_content, prepared)
         acquired_at = datetime.now(UTC)
+        observed_at = product.observed_at or acquired_at
         dataset = DatasetService.import_raster_bytes(
             db,
             project_id=project_id,
             area_id=payload.area_id,
             filename=filename,
             content=geotiff_content,
-            source="Digitaal Vlaanderen OMWRGBMRVL WMS",
+            source=f"Digitaal Vlaanderen WMS {product.layer}",
             source_name=OrthophotoAcquisitionService.PROVIDER,
-            source_version=f"most_recent_at_{acquired_at.date().isoformat()}",
+            temporal_series_key=f"digitaal-vlaanderen:orthophoto:{prepared['spatial_hash'][:24]}",
+            observed_at=observed_at,
+            valid_from=product.valid_from or observed_at,
+            valid_to=product.valid_to,
+            temporal_granularity=product.temporal_granularity,
+            source_version=(
+                f"most_recent_at_{acquired_at.date().isoformat()}"
+                if product.key == "most_recent"
+                else product.key
+            ),
             content_type="image/tiff",
             source_metadata={
                 "provider": OrthophotoAcquisitionService.PROVIDER,
                 "service": "WMS",
                 "service_version": "1.3.0",
-                "layer": resolved_settings.orthophoto_wms_layer,
-                "catalog_url": OrthophotoAcquisitionService.CATALOG_URL,
+                "product_key": product.key,
+                "product_display_name": product.display_name,
+                "observation_label": product.observation_label,
+                "observation_date_precision": product.temporal_granularity,
+                "native_resolution_m": product.native_resolution_m,
+                "requested_resolution_m": resolved_settings.orthophoto_resolution_m,
+                "color_mode": product.color_mode,
+                "supports_detection": product.supports_detection,
+                "layer": product.layer,
+                "catalog_url": product.catalog_url,
                 "attribution": OrthophotoAcquisitionService.ATTRIBUTION,
                 "license_note": "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen.",
             },
@@ -290,6 +488,7 @@ class OrthophotoAcquisitionService:
                 "acquisition": "explicit_bounded_map_selection",
                 "acquired_at": acquired_at.isoformat(),
                 "request_hash": prepared["request_hash"],
+                "spatial_hash": prepared["spatial_hash"],
                 "request_url": prepared["request_url"],
                 "response_content_type": response_content_type,
                 "bbox_epsg4326": prepared["bbox_epsg4326"],
@@ -297,19 +496,70 @@ class OrthophotoAcquisitionService:
                 "width": prepared["width"],
                 "height": prepared["height"],
                 "resolution_m": resolved_settings.orthophoto_resolution_m,
-                "limitation_message": OrthophotoAcquisitionService.LIMITATION,
+                "limitation_message": product.limitation_message,
             },
         )
         return OrthophotoAcquisitionResult(
             output_dataset_id=dataset.id,
             reused=False,
             provider=OrthophotoAcquisitionService.PROVIDER,
-            layer=resolved_settings.orthophoto_wms_layer,
+            product_key=product.key,
+            display_name=product.display_name,
+            observation_label=product.observation_label,
+            temporal_granularity=product.temporal_granularity,
+            supports_detection=product.supports_detection,
+            layer=product.layer,
             width=prepared["width"],
             height=prepared["height"],
             resolution_m=resolved_settings.orthophoto_resolution_m,
             bbox_epsg4326=prepared["bbox_epsg4326"],
             bbox_epsg31370=prepared["bbox_epsg31370"],
             attribution=OrthophotoAcquisitionService.ATTRIBUTION,
-            limitation_message=OrthophotoAcquisitionService.LIMITATION,
+            limitation_message=product.limitation_message,
         ).model_dump(mode="json")
+
+    @staticmethod
+    def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1600) -> bytes:
+        dataset = db.get(Dataset, dataset_id)
+        if (
+            dataset is None
+            or dataset.project_id != project_id
+            or dataset.source_name != OrthophotoAcquisitionService.PROVIDER
+            or dataset.status != "ready"
+            or not dataset.storage_path
+            or not Path(dataset.storage_path).is_file()
+        ):
+            raise AppError(code="ORTHOPHOTO_NOT_FOUND", message="Orthophoto dataset not found", status_code=404)
+        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="Raster preview dependencies are unavailable", 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))
+                indexes = [1] if source.count == 1 else list(range(1, min(source.count, 3) + 1))
+                pixels = source.read(indexes, out_shape=(len(indexes), height, width), resampling=Resampling.bilinear)
+                if pixels.dtype != np.uint8:
+                    pixels = np.clip(pixels, 0, 255).astype(np.uint8)
+                if len(indexes) == 1:
+                    image = Image.fromarray(pixels[0])
+                else:
+                    image = Image.fromarray(np.moveaxis(pixels[:3], 0, 2))
+                output = io.BytesIO()
+                image.save(output, format="PNG", optimize=True)
+                return output.getvalue()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="ORTHOPHOTO_PREVIEW_FAILED",
+                message="The persisted orthophoto could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc

backend/app/services/vector_feature_service.py

@@ -23,6 +23,7 @@ FULL_AREA_CLIPPED_OPERATOR_TOOLS = {
     "provision_official_landuse_timeseries.py",
     "provision_regional_grb_buildings.py",
     "provision_regional_grb_context.py",
+    "provision_waterinfo_station_history.py",
 }
 
 
@@ -439,7 +440,7 @@ class VectorFeatureService:
             length_m = db.query(func.coalesce(func.sum(length_expression), 0.0)).filter(*metric_filter).scalar()
             divisor = 1_000.0 if unit == "km" else 1.0
             metric_value = float(length_m or 0.0) / divisor
-        elif method in {"sum", "area_weighted_sum"}:
+        elif method in {"sum", "mean", "area_weighted_sum"}:
             property_name = str(config.get("property") or "").strip()
             if not property_name:
                 raise AppError(
@@ -457,8 +458,9 @@ class VectorFeatureService:
                 )
                 coverage_ratio = intersection_area / func.nullif(source_area, 0.0)
                 value_expression = numeric_value * coverage_ratio
+            aggregate_function = func.avg if method == "mean" else func.sum
             aggregate_value = (
-                db.query(func.coalesce(func.sum(value_expression), 0.0))
+                db.query(func.coalesce(aggregate_function(value_expression), 0.0))
                 .filter(*selection_filter)
                 .filter(VectorFeature.properties_json.op("->>")(property_name).isnot(None))
                 .scalar()

backend/tests/test_sprint196_map_orthophoto_analysis.py

@@ -86,7 +86,13 @@ class FakeImageResponse:
         return self.content[:limit]
 
 
-def _selection_payload(*, side_m: float = 512.0, force_refresh: bool = True, area_id=None) -> OrthophotoAcquireRequest:
+def _selection_payload(
+    *,
+    side_m: float = 512.0,
+    force_refresh: bool = True,
+    area_id=None,
+    product_key: str = "most_recent",
+) -> OrthophotoAcquireRequest:
     west, south = 199_000.0, 210_000.0
     transformer = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True)
     min_lon, min_lat = transformer.transform(west, south)
@@ -100,6 +106,7 @@ def _selection_payload(*, side_m: float = 512.0, force_refresh: bool = True, are
             "crs": "EPSG:4326",
         },
         area_id=area_id,
+        product_key=product_key,
         force_refresh=force_refresh,
     )
 
@@ -136,6 +143,25 @@ def test_orthophoto_request_is_bounded_and_uses_official_wms_contract() -> None:
     assert len(prepared["request_hash"]) == 64
 
 
+def test_orthophoto_product_registry_exposes_only_governed_official_layers() -> None:
+    settings = Settings(_env_file=None)
+    products = OrthophotoAcquisitionService.list_products(settings)
+    keys = [item["key"] for item in products]
+
+    assert keys[0] == "most_recent"
+    assert {"2025", "2012", "2008_2011", "2000_2003", "1979_1990", "1971"}.issubset(keys)
+    assert next(item for item in products if item["key"] == "most_recent")["supports_detection"] is True
+    assert all(item["supports_detection"] is False for item in products if item["key"] != "most_recent")
+
+    prepared = OrthophotoAcquisitionService._prepared_request(_selection_payload(product_key="1971"), settings)
+    assert prepared["params"]["LAYERS"] == "OKZPAN71VL"
+    assert prepared["wms_url"] == "https://geo.api.vlaanderen.be/OKZ/wms"
+
+    with pytest.raises(AppError) as exc_info:
+        OrthophotoAcquisitionService._prepared_request(_selection_payload(product_key="arbitrary-layer"), settings)
+    assert exc_info.value.code == "ORTHOPHOTO_PRODUCT_NOT_SUPPORTED"
+
+
 @pytest.mark.parametrize(
     ("side_m", "expected_code"),
     [(64.0, "ORTHOPHOTO_SELECTION_TOO_SMALL"), (1_200.0, "ORTHOPHOTO_SELECTION_TOO_LARGE")],
@@ -240,7 +266,7 @@ def test_orthophoto_acquisition_reuses_fresh_exact_request_without_provider_call
     cached = Dataset(
         id=uuid4(),
         project_id=project_id,
-        name=f"orthofoto_selectie_{prepared['request_hash'][:12]}.tif",
+        name=f"orthofoto_most_recent_{prepared['request_hash'][:12]}.tif",
         dataset_type="raster",
         source="Digitaal Vlaanderen",
         source_name="digitaal_vlaanderen_orthophoto",
@@ -277,6 +303,54 @@ def test_orthophoto_provider_rejects_non_image_response() -> None:
     assert exc_info.value.code == "ORTHOPHOTO_PROVIDER_INVALID_RESPONSE"
 
 
+def test_historical_orthophoto_persists_temporal_product_provenance(tmp_path) -> None:
+    project_id = uuid4()
+    payload = _selection_payload(product_key="2020")
+    db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")})
+    settings = Settings(_env_file=None, storage_root=str(tmp_path), orthophoto_resolution_m=1.0)
+    prepared = OrthophotoAcquisitionService._prepared_request(payload, settings)
+    response = FakeImageResponse(_source_tiff(prepared["width"], prepared["height"]))
+
+    result = OrthophotoAcquisitionService.acquire(
+        db,
+        project_id,
+        payload,
+        settings=settings,
+        opener=lambda *_args, **_kwargs: response,
+    )
+
+    dataset = next(row for row in db.added if isinstance(row, Dataset))
+    assert result["product_key"] == "2020"
+    assert result["supports_detection"] is False
+    assert dataset.observed_at.year == 2020
+    assert dataset.temporal_granularity == "year"
+    assert dataset.source_metadata["layer"] == "OMWRGB20VL"
+    assert dataset.source_metadata["product_key"] == "2020"
+    assert dataset.provenance_metadata["spatial_hash"] == prepared["spatial_hash"]
+
+
+def test_persisted_orthophoto_renders_browser_png(tmp_path) -> None:
+    project_id = uuid4()
+    dataset_id = uuid4()
+    path = tmp_path / "ortho.tif"
+    path.write_bytes(_source_tiff(32, 24))
+    dataset = Dataset(
+        id=dataset_id,
+        project_id=project_id,
+        name="ortho.tif",
+        dataset_type="raster",
+        source="Digitaal Vlaanderen",
+        source_name="digitaal_vlaanderen_orthophoto",
+        status="ready",
+        storage_path=str(path),
+    )
+    db = FakeSession({(Dataset, dataset_id): dataset})
+
+    png = OrthophotoAcquisitionService.render_png(db, project_id, dataset_id)
+
+    assert png.startswith(b"\x89PNG\r\n\x1a\n")
+
+
 def test_orthophoto_endpoint_returns_canonical_job_envelope(monkeypatch) -> None:
     project_id = uuid4()
     output_dataset_id = uuid4()
@@ -307,6 +381,22 @@ def test_orthophoto_endpoint_returns_canonical_job_envelope(monkeypatch) -> None
     assert any(isinstance(row, Job) for row in db.added)
 
 
+def test_orthophoto_product_endpoint_returns_canonical_envelope() -> None:
+    project_id = uuid4()
+    db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")})
+    app.dependency_overrides[get_db] = lambda: db
+    try:
+        response = TestClient(app).get(f"/api/v1/projects/{project_id}/datasets/orthophoto/products")
+    finally:
+        app.dependency_overrides.clear()
+
+    assert response.status_code == 200
+    body = response.json()
+    assert set(body) == {"data"}
+    assert body["data"]["total"] == len(body["data"]["items"])
+    assert body["data"]["items"][0]["key"] == "most_recent"
+
+
 def test_frontend_connects_map_selection_to_existing_detection_and_qa_flows() -> None:
     app_source = (ROOT / "frontend" / "src" / "App.tsx").read_text(encoding="utf-8")
     hook_source = (ROOT / "frontend" / "src" / "hooks" / "useMapOrthophotoAnalysis.ts").read_text(encoding="utf-8")

backend/tests/test_sprint201_semantic_selection_metrics.py

@@ -114,6 +114,36 @@ def test_population_keeps_configured_metric_and_adds_sector_count() -> None:
     }
 
 
+def test_station_measurement_uses_numeric_mean_without_area_extrapolation() -> None:
+    dataset = themed_dataset("water", method="mean")
+    dataset.source_name = "waterinfo"
+    dataset.source_metadata.update(
+        {
+            "semantic_metrics": False,
+            "selection_aggregation": {
+                "metric_key": "water_level",
+                "method": "mean",
+                "property": "annual_mean_water_level_m",
+                "label": "Jaargemiddelde waterstand",
+                "unit": "m",
+                "warning": "Puntmeting; geen gebiedsdekkend watervolume.",
+            },
+        }
+    )
+
+    result = VectorFeatureService.summarize_features_by_bbox(
+        SequenceScalarSession([30.455]),
+        dataset=dataset,
+        bbox=BBOX,
+        total_feature_count=1,
+    )
+
+    assert result["metric_value"] == 30.455
+    assert result["aggregation_method"] == "mean"
+    assert result["metric_unit"] == "m"
+    assert result["warning"] == "Puntmeting; geen gebiedsdekkend watervolume."
+
+
 def test_future_regional_imports_persist_semantic_aggregation_configuration() -> None:
     buildings = (ROOT / "scripts/provision_regional_grb_buildings.py").read_text(encoding="utf-8")
     context = (ROOT / "scripts/provision_regional_grb_context.py").read_text(encoding="utf-8")

backend/tests/test_sprint203_waterinfo_history.py

@@ -0,0 +1,135 @@
+from __future__ import annotations
+
+import importlib.util
+import sys
+from pathlib import Path
+
+from shapely.geometry import box
+
+
+ROOT = Path(__file__).resolve().parents[2]
+
+
+def load_operator():
+    script_path = ROOT / "scripts" / "provision_waterinfo_station_history.py"
+    spec = importlib.util.spec_from_file_location("waterinfo_history_operator", script_path)
+    assert spec is not None
+    assert spec.loader is not None
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[spec.name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+class JsonResponse:
+    ok = True
+    status_code = 200
+    text = ""
+
+    def __init__(self, payload):
+        self.payload = payload
+
+    def raise_for_status(self):
+        return None
+
+    def json(self):
+        return self.payload
+
+
+class JsonSession:
+    def __init__(self, payloads):
+        self.payloads = iter(payloads)
+        self.calls = []
+
+    def get(self, url, *, params, timeout):
+        self.calls.append((url, params, timeout))
+        return JsonResponse(next(self.payloads))
+
+
+def test_waterinfo_station_discovery_filters_exact_area_and_uses_annual_group() -> None:
+    module = load_operator()
+    payload = {
+        "type": "FeatureCollection",
+        "features": [
+            {
+                "type": "Feature",
+                "geometry": {"type": "Point", "coordinates": [5.1, 51.2]},
+                "properties": {"ts_id": 5319042, "station_no": "L10_089", "station_name": "Mol/ScheppelijkeNete"},
+            },
+            {
+                "type": "Feature",
+                "geometry": {"type": "Point", "coordinates": [6.0, 52.0]},
+                "properties": {"ts_id": 999, "station_no": "outside", "station_name": "Outside"},
+            },
+        ],
+    }
+    session = JsonSession([payload])
+
+    raw, stations = module.discover_station_series(
+        session,
+        module.PARAMETERS["water_level"],
+        box(5.0, 51.0, 5.3, 51.4),
+        timeout=30,
+    )
+
+    assert raw == payload
+    assert [item["ts_id"] for item in stations] == ["5319042"]
+    assert session.calls[0][1]["timeseriesgroup_id"] == "192784"
+    assert session.calls[0][1]["request"] == "getTimeseriesValueLayer"
+
+
+def test_waterinfo_annual_values_reject_invalid_sentinel_and_keep_real_zero() -> None:
+    module = load_operator()
+    payload = [
+        {
+            "ts_id": 5319042,
+            "data": [
+                ["2013-01-01T00:00:00.000+01:00", 30.46],
+                ["2014-01-01T00:00:00.000+01:00", -9999],
+                ["2015-01-01T00:00:00.000+01:00", 0.0],
+                ["2026-01-01T00:00:00.000+01:00", 99.0],
+            ],
+        }
+    ]
+    session = JsonSession([payload])
+
+    raw, values = module.fetch_annual_values(session, "5319042", from_year=2013, to_year=2025, timeout=30)
+
+    assert raw == payload
+    assert values == {2013: 30.46, 2015: 0.0}
+    assert session.calls[0][1]["request"] == "getTimeseriesValues"
+
+
+def test_waterinfo_snapshot_and_series_keep_station_identity_and_honest_metric() -> None:
+    module = load_operator()
+    parameter = module.PARAMETERS["water_level"]
+    station = {
+        "ts_id": "5319042",
+        "geometry": {"type": "Point", "coordinates": [5.1, 51.2]},
+        "properties": {
+            "station_id": "123",
+            "station_no": "L10_089",
+            "station_name": "Mol/ScheppelijkeNete",
+            "ts_unitsymbol": "m",
+        },
+    }
+
+    snapshot = module.build_snapshot(parameter, station, 2025, 30.455)
+
+    feature = snapshot["features"][0]
+    assert module.series_key(parameter, station) == "waterinfo:water_level:annual:l10-089"
+    assert feature["geometry"]["type"] == "Point"
+    assert feature["properties"]["annual_mean_water_level_m"] == 30.455
+    assert feature["properties"]["timeseries_id"] == "5319042"
+    assert "volume" in parameter.limitation
+
+
+def test_waterinfo_operator_is_packaged_and_readiness_checked() -> None:
+    readiness = (ROOT / "scripts" / "run_readiness_check.sh").read_text(encoding="utf-8")
+    dockerfile = (ROOT / "deploy" / "unraid" / "Dockerfile.all-in-one").read_text(encoding="utf-8")
+    vector_service = (ROOT / "backend" / "app" / "services" / "vector_feature_service.py").read_text(encoding="utf-8")
+
+    assert "py_compile scripts/provision_waterinfo_station_history.py" in readiness
+    assert "COPY scripts/provision_waterinfo_station_history.py" in dockerfile
+    assert '"provision_waterinfo_station_history.py"' in vector_service
+    assert '"sum", "mean", "area_weighted_sum"' in vector_service

deploy/unraid/Dockerfile.all-in-one

@@ -77,6 +77,7 @@ COPY scripts/provision_mol_context_layers.py /app/scripts/provision_mol_context_
 COPY scripts/provision_mol_population_history.py /app/scripts/provision_mol_population_history.py
 COPY scripts/provision_mol_historical_landuse.py /app/scripts/provision_mol_historical_landuse.py
 COPY scripts/provision_official_landuse_timeseries.py /app/scripts/provision_official_landuse_timeseries.py
+COPY scripts/provision_waterinfo_station_history.py /app/scripts/provision_waterinfo_station_history.py
 COPY scripts/provision_regional_timeseries.py /app/scripts/provision_regional_timeseries.py
 COPY scripts/geographic_scopes.py /app/scripts/geographic_scopes.py
 COPY scripts/provision_geographic_scope.py /app/scripts/provision_geographic_scope.py

docs/API_CONTRACTS.md

@@ -191,15 +191,24 @@ Response: `DatasetRead` with extracted metadata if supported.
 
 Vector uploads remain stored as original files and are also persisted into `vector_features` as queryable PostGIS state.
 
+### GET `/api/v1/projects/{project_id}/datasets/orthophoto/products`
+
+Return the governed Digitaal Vlaanderen orthophoto product allowlist in the
+canonical envelope. Every product reports its key, display/observation label,
+temporal granularity, native resolution, colour mode, catalogue URL,
+limitations and whether current configured-YOLO detection is allowed.
+
 ### POST `/api/v1/projects/{project_id}/datasets/orthophoto/acquire`
 
-Explicitly acquire a bounded most-recent winter orthophoto selection from the
-official Digitaal Vlaanderen `OMWRGBMRVL` WMS `Ortho` layer.
+Explicitly acquire a bounded orthophoto selection from a governed official
+Digitaal Vlaanderen WMS product. Arbitrary WMS URLs and layer names are not
+accepted.
 
 ```json
 {
   "bbox": {"min_x": 5.10, "min_y": 51.17, "max_x": 5.11, "max_y": 51.18, "crs": "EPSG:4326"},
   "area_id": "optional-project-area-uuid",
+  "product_key": "most_recent",
   "force_refresh": false
 }
 ```
@@ -207,7 +216,8 @@ official Digitaal Vlaanderen `OMWRGBMRVL` WMS `Ortho` layer.
 The canonical envelope contains a synchronous Job. Its `output_dataset_id`
 identifies the raster Dataset; `result_json` contains provider, layer, pixel
 dimensions, EPSG:4326/EPSG:31370 bounds, sampling resolution, attribution,
-cache reuse and limitation text.
+cache reuse and limitation text. Historical products also persist their
+observation/validity period and a spatially scoped temporal-series key.
 
 Safety contract:
 
@@ -218,8 +228,11 @@ Safety contract:
   reuse;
 - WMS bytes are georeferenced to EPSG:31370 and persisted only through
   `DatasetService`; no fetch runs on startup;
-- this is the latest mosaic available at request time, not a historical
-  observation date for every pixel.
+- only `most_recent` can enter the current configured-YOLO plus GRB-QA path;
+  historical products are visual evidence and are never validated against the
+  current GRB state;
+- product periods such as `1979_1990` remain explicitly multi-year and are not
+  presented as exact annual observations.
 
 ### GET `/api/v1/projects/{project_id}/datasets`
 
@@ -274,6 +287,13 @@ If preview dependencies are unavailable:
 - code: `RASTER_PROCESSING_UNAVAILABLE`
 - message: `Raster preview unavailable...`
 
+### GET `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/image`
+
+Return a persisted orthophoto Dataset as a bounded browser-safe PNG. This is an
+explicit binary non-envelope endpoint used by the MapLibre image source. It
+accepts only ready datasets from the governed orthophoto provider and never
+reads arbitrary filesystem paths.
+
 ### POST `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/clip`
 
 Clip raster by selected area. Returns a `202`-style accepted job payload through the job wrapper (`jobs` create/read flow).

docs/CODEX_EXECUTION_LOG.md

@@ -8438,3 +8438,33 @@ Next:
 - Integrate Waterinfo/VMM station observations as point time series and add
   historical orthophoto acquisition, while preserving their spatial and
   methodological limitations.
+
+## Sprint 203 - Governed Waterinfo history and historical orthophotos (2026-07-15)
+
+Implemented:
+- Added `provision_waterinfo_station_history.py` using the documented Waterinfo
+  KiWIS annual water-level/discharge groups. The operator filters station
+  points against the exact persisted Area, retains raw JSON plus SHA256
+  manifests and writes only through the canonical upload API.
+- Persisted one Point Dataset per station/year and kept station identities in
+  separate temporal series. Added backend `mean` aggregation without combining
+  stations or inferring area-wide water level/volume.
+- Added a governed orthophoto allowlist for the most-recent product, annual
+  winter mosaics 2012-2025, older winter periods, RGB 1979-1990 and
+  panchromatic 1971. Arbitrary WMS URLs/layers remain impossible.
+- Added raster temporal metadata and a constrained PNG rendering endpoint for
+  persisted orthophoto datasets. The map can select and display official
+  historical imagery over the same bounded rectangle.
+- Historical products explicitly bypass configured-YOLO and current-GRB QA;
+  only `most_recent` retains that path.
+- Updated source inventory behavior and wrote a governed source backlog with
+  acceptance criteria for BWK/Natura 2000, agricultural parcels, Buildings
+  Register, DHMV and bathymetry.
+
+Validation evidence:
+- Focused Waterinfo/orthophoto/semantic metric tests and frontend typecheck/build
+  passed before the full repository gate.
+
+Remaining operational step:
+- Deploy the all-in-one image, provision the real Mol Waterinfo station series
+  and verify one historical image acquisition/overlay through the LAN browser.

docs/DATA_SOURCES.md

@@ -214,9 +214,6 @@ These sources are available from their public authorities but are not silently
 treated as loaded GeoIntel data. The Source inventory labels them separately
 until a governed operator import, provenance record and validation pass exist.
 
-- Historical orthophotos (Digitaal Vlaanderen): 1971 and 1979-1990 through
-  the `OKZ` WMS, with additional dated mosaics as separate products. Suitable
-  for visual/image evolution after a bounded acquisition contract is added.
 - Biologische Waarderingskaart / Natura 2000 (INBO), state 2025: suitable for
   habitat, biotope and ecological-value analysis, not a continuous annual
   series.
@@ -229,9 +226,35 @@ until a governed operator import, provenance record and validation pass exist.
 - DHMV II DTM/DSM (Digitaal Vlaanderen): 1 m/5 m elevation based on 2013-2015
   LiDAR, suitable for elevation, slope and drainage. It does not provide water
   depth.
-- Waterinfo/VMM: station time series for water level, flow and precipitation.
-  These can describe hydrological state, but do not provide area-wide water
-  volume without a compatible bottom profile/bathymetry model.
+
+## Governed Waterinfo station history
+
+`scripts/provision_waterinfo_station_history.py` uses the public Waterinfo
+KiWIS query service only after an explicit operator command. It discovers the
+documented annual water-level (`192784`) and discharge (`192895`) groups,
+filters station points against the exact persisted Area and retains the raw
+station/value JSON plus SHA256 manifests.
+
+Each station and year becomes one immutable reference Dataset through the
+normal upload API. Series keys include the station identity; different stations
+are never averaged into one municipal value. Selection aggregation is a numeric
+mean over the selected station records and remains labelled as a point
+measurement. Water level or discharge does not establish area-wide water
+volume without compatible depth, profile and coverage data.
+
+## Governed historical orthophotos
+
+The bounded map acquisition registry includes the official annual winter
+mosaics for 2012-2025, period products for 2000-2003, 2005-2007 and 2008-2011,
+the RGB 1979-1990 summer mosaic and the panchromatic 1971 mosaic. Requests stay
+within the configured 128-1,024 m safety envelope and are persisted as
+EPSG:31370 raster Datasets with explicit product, layer, observation period,
+attribution and request hashes.
+
+Historical mosaics are for visual comparison only in this phase. Current YOLO
+building detection and current GRB QA remain restricted to `most_recent`, since
+validating an old image against today's building state would produce dishonest
+quality metrics.
 
 ## OSM
 

docs/DATA_SPECIFICATION.md

@@ -198,6 +198,21 @@ Metadata:
 - nodata
 - resolution
 
+Temporal orthophoto rasters additionally require a governed product key, WMS
+layer, observation label, `observed_at`, optional `valid_from`/`valid_to`,
+temporal granularity, request/spatial hash, attribution and a limitation that
+states whether the product is annual, multi-year or merely most recent.
+
+### Hydrological station observations
+
+Waterinfo observations are persisted as EPSG:4326 Point features, one station
+and one observation year per immutable Dataset. Required properties are the
+station/timeseries identity, measurement type, numeric annual value, reported
+unit, observation year, provider/owner and attribution. Different station
+series may not be merged into one area-wide value. Point water level and
+discharge may not be converted to water volume without governed compatible
+depth/profile data.
+
 ### Vector
 
 Ondersteund:

docs/STORAGE_ARCHITECTURE.md

@@ -99,6 +99,16 @@ Mask files are provenance/debug artifacts. QA, map display and GeoJSON output mu
 
 Do not delete originals automatically. Derived outputs may be cleaned through explicit cache management.
 
+## Official temporal source artifacts
+
+Waterinfo raw station layers, timeseries responses and checksum manifests live
+under `storage/operator-data/waterinfo/<scope>/`. These are immutable source
+evidence; queryable annual Point snapshots are normal Dataset/vector_feature
+records. Bounded orthophotos are normal raster Dataset files. Their WMS URL,
+product/layer, request/spatial hash, temporal validity and limitations are held
+in source/provenance metadata. Browser PNG rendering is derived on request and
+does not replace the stored GeoTIFF.
+
 Offline demo export artifacts can be inspected and cleaned with:
 
 ```bash

docs/TODO.md

@@ -22,12 +22,55 @@
 - [x] Add a source inventory that separates loaded data from audited official follow-up sources.
 - [x] Add a local Ollama question window grounded in persisted GeoIntel metrics and installed server models.
 - [ ] Add a governed depth/bathymetry source before exposing water volume; never infer volume from 2D GRB water geometry.
-- [ ] Integrate one governed hydrology source (Waterinfo/VMM station series) without presenting point measurements as area-wide water volume.
-- [ ] Add bounded historical orthophoto acquisition and visual change analysis after validating layer/year coverage.
+- [x] Integrate governed Waterinfo/VMM annual station series without presenting point measurements as area-wide water volume.
+- [x] Add bounded historical orthophoto acquisition for official 1971-2025 products with a map overlay and no current-GRB QA on old imagery.
 - [ ] Add BWK/Natura 2000 and annual agricultural-use parcels through explicit provider/operator contracts.
 - [ ] Extend the official 1778/1873/1969 historical buildings, water and roads series from Mol to the approved regional scope with partitioned source audits.
 - [x] Connect a drawn rectangle to bounded official orthophoto acquisition, local configured-YOLO detection and persisted GRB QA.
 
+## Governed source expansion backlog
+
+Implement these in order. Every source must use an explicit operator/provider
+contract, retain raw checksummed evidence, persist through DatasetService and
+pass live Mol validation before regional expansion.
+
+### P1 - Biologische Waarderingskaart / Natura 2000
+
+- [ ] Confirm the official downloadable service/layer/version and licence for the 2025 state.
+- [ ] Define stable habitat/biotope/value fields and keep BWK valuation separate from Natura 2000 habitat classification.
+- [ ] Clip in EPSG:31370 to Mol, transform to EPSG:4326 and persist valid polygonal `vector_features` with source feature ids.
+- [ ] Add hectare metrics per governed class plus unknown/unmapped-class diagnostics; do not invent an annual trend from one state.
+- [ ] Add source manifest, unit tests, live row/count/geometry audit and map/source-inventory presentation.
+
+### P2 - Annual agricultural-use parcels
+
+- [ ] Inventory official annual editions and document schema/code-list changes before selecting a comparable year range.
+- [ ] Retain annual source files and crop code lists; normalize only fields whose meaning is stable across editions.
+- [ ] Persist separate annual Datasets with hectares by crop/use class and an explicit unstable-parcel-identity limitation.
+- [ ] Compare area/category totals over time; do not claim parcel lineage where identifiers or boundaries changed.
+- [ ] Validate Mol first, then partition all 28 Kempen municipalities with completeness/checksum manifests.
+
+### P3 - Buildings and Addresses Register
+
+- [ ] Define a governed snapshot operator and relation mapping between building unit, building object, address and GRB geometry.
+- [ ] Keep register lifecycle/status semantics separate from GRB footprint geometry and document match confidence/unmatched rows.
+- [ ] Add current building-status/address metrics without exposing personal data or treating addresses as households/population.
+- [ ] Add Mol reconciliation tests and a live mismatch report before regional import.
+
+### P4 - Digitaal Hoogtemodel Vlaanderen
+
+- [ ] Select official DTM/DSM product, service and native resolution; retain acquisition date and vertical reference.
+- [ ] Add bounded raster acquisition/clip storage with nodata, CRS, resolution and checksum validation.
+- [ ] Implement governed elevation, relief and slope statistics in metres/degrees; keep drainage interpretation explicitly derived.
+- [ ] Do not label terrain/surface height as water depth and do not enable volume from DHMV alone.
+
+### P5 - Water depth / bathymetry
+
+- [ ] Identify an authoritative source with compatible spatial coverage, vertical datum, date and uncertainty; otherwise keep volume unavailable.
+- [ ] Define waterbody linkage, surface elevation, bottom elevation and uncertainty propagation before adding any volume metric.
+- [ ] Validate coverage gaps and prohibit extrapolation outside measured/profiled waterbodies.
+- [ ] Add independent GIS review and golden-volume fixtures before exposing the result to users or Ollama.
+
 This file now starts with the current implementation status. Older preparation/backlog sections are preserved below as historical planning context and should not be treated as the live sprint board without checking `docs/CODEX_EXECUTION_LOG.md`.
 
 ## Release hardening status

frontend/README.md

@@ -404,6 +404,13 @@ their real observation ranges. A collapsed follow-up catalogue distinguishes
 official sources that exist from datasets that are already persisted in the
 active project.
 
+The Map result panel also loads a governed orthophoto product list. The most
+recent product retains the configured-YOLO plus current-GRB QA action. Official
+historical years/periods use the same bounded rectangle, persist as raster
+Datasets and render as MapLibre image overlays, but expose no misleading
+current-state AI/QA action. The Sources inventory moves Waterinfo and historical
+orthophotos from follow-up to loaded evidence only after such datasets exist.
+
 Evolution mode compares the exact selected persisted Area when the full
 municipality/region action is used. It shows the selected before/after values,
 all compatible supporting metrics and a chart/table for every observation in

frontend/src/App.tsx

@@ -1018,6 +1018,10 @@ function App(): JSX.Element {
               orthophotoAnalysisRunning={mapOrthophotoAnalysis.running}
               orthophotoAnalysisQuality={mapOrthophotoAnalysis.lastQuality}
               orthophotoAnalysisDetectionCount={mapOrthophotoAnalysis.lastDetectionCount}
+              orthophotoProducts={mapOrthophotoAnalysis.products}
+              selectedOrthophotoProductKey={mapOrthophotoAnalysis.selectedProductKey}
+              orthophotoResult={mapOrthophotoAnalysis.lastResult}
+              orthophotoImageUrl={mapOrthophotoAnalysis.imageUrl}
               availableMapDatasets={availableMapDatasets}
               selectedMapDatasetId={selectedDataset && isVectorDatasetType(selectedDataset.dataset_type) ? selectedDataset.id : ''}
               selectedFeature={selectedMapFeature}
@@ -1039,6 +1043,7 @@ function App(): JSX.Element {
               onRunMapSelectionQa={runMapSelectionQa}
               onOpenMapSelectionQualityEvidence={openMapSelectionQualityEvidence}
               onRunOrthophotoAnalysis={mapOrthophotoAnalysis.run}
+              onSelectOrthophotoProduct={mapOrthophotoAnalysis.setSelectedProductKey}
               onClearQualityEvidence={clearQualityEvidenceGeoJson}
             />
           ) : null}

frontend/src/components/GeoMap.tsx

@@ -3,7 +3,7 @@ import maplibregl from 'maplibre-gl'
 import 'maplibre-gl/dist/maplibre-gl.css'
 import { PRIMARY_FOCUS_CENTER } from '../config/primaryFocus'
 import { featureCollectionBounds } from '../lib/geojsonBounds'
-import type { MapViewportState, VectorSelectionBBox } from '../types'
+import type { MapImageOverlay, MapViewportState, VectorSelectionBBox } from '../types'
 
 interface GeoMapProps {
   data: GeoJSON.FeatureCollection | null
@@ -13,6 +13,7 @@ interface GeoMapProps {
   selectedFeature?: GeoJSON.Feature | null
   selectionData?: GeoJSON.FeatureCollection | null
   qaEvidenceData?: GeoJSON.FeatureCollection | null
+  imageOverlay?: MapImageOverlay | null
   selectionBbox?: { min_x: number; min_y: number; max_x: number; max_y: number } | null
   bboxSelectionMode?: boolean
   visible?: boolean
@@ -153,6 +154,7 @@ function GeoMap({
   selectedFeature = null,
   selectionData = null,
   qaEvidenceData = null,
+  imageOverlay = null,
   selectionBbox = null,
   bboxSelectionMode = false,
   visible = true,
@@ -346,6 +348,43 @@ function GeoMap({
     }
   }, [])
 
+  useEffect(() => {
+    const map = mapRef.current
+    if (!map || !mapStyleReady || !map.isStyleLoaded()) {
+      return
+    }
+    if (map.getLayer('bounded-orthophoto')) {
+      map.removeLayer('bounded-orthophoto')
+    }
+    if (map.getSource('bounded-orthophoto')) {
+      map.removeSource('bounded-orthophoto')
+    }
+    if (!imageOverlay) {
+      return
+    }
+    const [minX, minY, maxX, maxY] = imageOverlay.bbox
+    map.addSource('bounded-orthophoto', {
+      type: 'image',
+      url: imageOverlay.url,
+      coordinates: [
+        [minX, maxY],
+        [maxX, maxY],
+        [maxX, minY],
+        [minX, minY],
+      ],
+    })
+    const beforeLayer = ['area-fill', 'dataset-fill', 'selection-bbox-fill'].find((layerId) => map.getLayer(layerId))
+    map.addLayer(
+      {
+        id: 'bounded-orthophoto',
+        type: 'raster',
+        source: 'bounded-orthophoto',
+        paint: { 'raster-opacity': imageOverlay.opacity ?? 0.88 },
+      },
+      beforeLayer,
+    )
+  }, [imageOverlay, mapStyleReady])
+
   useEffect(() => {
     const map = mapRef.current
     if (!map || !mapStyleReady || !map.isStyleLoaded()) {

frontend/src/components/datasets/SourceCatalogPanel.tsx

@@ -15,13 +15,15 @@ const THEME_LABELS: Record<string, string> = {
 
 const AVAILABLE_SOURCES = [
   {
+    key: 'historical_orthophoto',
     name: 'Historische orthofoto’s',
     owner: 'Digitaal Vlaanderen',
-    coverage: '1971 en 1979-1990; aanvullende jaargangen bestaan afzonderlijk',
-    value: 'Visuele evolutie en toekomstige beeldvergelijking',
+    coverage: '1971, 1979-1990, 2000-2025',
+    value: 'Visuele evolutie via begrensde officiële luchtbeelden',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/orthofotomozaiek-kleinschalig-zomeropnamen',
   },
   {
+    key: 'bwk',
     name: 'Biologische Waarderingskaart / Natura 2000',
     owner: 'INBO',
     coverage: 'Toestand 2025',
@@ -29,6 +31,7 @@ const AVAILABLE_SOURCES = [
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/biologische-waarderingskaart-en-natura-2000-habitatkaart-toestand-2025',
   },
   {
+    key: 'agriculture',
     name: 'Landbouwgebruikspercelen',
     owner: 'Agentschap Landbouw en Zeevisserij',
     coverage: 'Jaarlijkse bestanden',
@@ -36,6 +39,7 @@ const AVAILABLE_SOURCES = [
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/open-geodata-landbouwgebruikspercelen',
   },
   {
+    key: 'buildings_register',
     name: 'Gebouwen- en adressenregister',
     owner: 'Digitaal Vlaanderen',
     coverage: 'Continu geactualiseerd',
@@ -43,6 +47,7 @@ const AVAILABLE_SOURCES = [
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/gebouwen-en-adressenregister',
   },
   {
+    key: 'elevation',
     name: 'Digitaal Hoogtemodel Vlaanderen II',
     owner: 'Digitaal Vlaanderen',
     coverage: 'LiDAR-opname 2013-2015, DTM/DSM 1 m en 5 m',
@@ -50,6 +55,7 @@ const AVAILABLE_SOURCES = [
     url: 'https://www.vlaanderen.be/digitaal-vlaanderen/onze-diensten-en-platformen/earth-observation-data-science-eodas/het-digitaal-hoogtemodel/digitaal-hoogtemodel-vlaanderen-ii',
   },
   {
+    key: 'waterinfo',
     name: 'Waterinfo en VMM-metingen',
     owner: 'Vlaamse Milieumaatschappij',
     coverage: 'Meetpunten en tijdreeksen voor waterstand, debiet en neerslag',
@@ -91,6 +97,17 @@ function timelineSummary(
 
 export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.Element {
   const ready = datasets.filter((dataset) => dataset.status === 'ready')
+  const waterinfoDatasets = ready.filter((dataset) => dataset.source_name === 'waterinfo')
+  const historicalOrthophotos = ready.filter(
+    (dataset) =>
+      dataset.source_name === 'digitaal_vlaanderen_orthophoto' &&
+      String(dataset.source_metadata?.['product_key'] ?? 'most_recent') !== 'most_recent',
+  )
+  const pendingSources = AVAILABLE_SOURCES.filter((source) => {
+    if (source.key === 'waterinfo') return waterinfoDatasets.length === 0
+    if (source.key === 'historical_orthophoto') return historicalOrthophotos.length === 0
+    return true
+  })
   const themes = Object.keys(THEME_LABELS).map((theme) => {
     const matches = ready.filter((dataset) => datasetTheme(dataset) === theme)
     const temporal = matches.filter((dataset) => dataset.temporal_series_key && dataset.observed_at)
@@ -148,10 +165,29 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E
         ))}
       </div>
 
+      {waterinfoDatasets.length > 0 || historicalOrthophotos.length > 0 ? (
+        <div className="source-catalog-loaded" aria-label="Aanvullende ingeladen bronnen">
+          {waterinfoDatasets.length > 0 ? (
+            <article>
+              <strong>Waterinfo meetreeksen</strong>
+              <span>{new Set(waterinfoDatasets.map((dataset) => dataset.temporal_series_key).filter(Boolean)).size} stationsreeksen · {waterinfoDatasets.length} jaarmetingen</span>
+              <p>Puntmetingen blijven afzonderlijk per station en worden niet als gebiedsgemiddelde of watervolume voorgesteld.</p>
+            </article>
+          ) : null}
+          {historicalOrthophotos.length > 0 ? (
+            <article>
+              <strong>Historische luchtbeelden</strong>
+              <span>{historicalOrthophotos.length} begrensde kaartselecties bewaard</span>
+              <p>Officiële jaargangen en periodes zijn via de kaart beschikbaar zonder actuele GRB-validatie.</p>
+            </article>
+          ) : null}
+        </div>
+      ) : null}
+
       <details className="source-opportunity-list">
         <summary>Officiële bronnen die hierna kunnen worden ingeladen</summary>
         <div>
-          {AVAILABLE_SOURCES.map((source) => (
+          {pendingSources.map((source) => (
             <article key={source.name}>
               <div>
                 <strong>{source.name}</strong>

frontend/src/components/map/MapWorkspace.tsx

@@ -1,6 +1,6 @@
 import { useEffect, useMemo, useState } from 'react'
 import GeoMap from '../GeoMap'
-import type { AreaRead, DatasetCreateResponse, DetectionQaResult, MapViewportState, ProjectRead, QaComparisonResult, VectorSelectionBBox, VectorSelectionMetric, VectorSelectionResponse } from '../../types'
+import type { AreaRead, DatasetCreateResponse, DetectionQaResult, MapViewportState, OrthophotoAcquisitionResult, OrthophotoProductRead, ProjectRead, QaComparisonResult, VectorSelectionBBox, VectorSelectionMetric, VectorSelectionResponse } from '../../types'
 import { featureCollectionBounds } from '../../lib/geojsonBounds'
 import { useMapThemeSelectionInsights } from '../../hooks/useMapThemeSelectionInsights'
 import { useTemporalComparison } from '../../hooks/useTemporalComparison'
@@ -459,6 +459,10 @@ interface MapWorkspaceProps {
   orthophotoAnalysisRunning: boolean
   orthophotoAnalysisQuality: DetectionQaResult | null
   orthophotoAnalysisDetectionCount: number | null
+  orthophotoProducts: OrthophotoProductRead[]
+  selectedOrthophotoProductKey: string
+  orthophotoResult: OrthophotoAcquisitionResult | null
+  orthophotoImageUrl: string | null
   availableMapDatasets: DatasetCreateResponse[]
   selectedMapDatasetId: string
   onSelectMapArea: (areaId: string) => void
@@ -479,6 +483,7 @@ interface MapWorkspaceProps {
   onRunMapSelectionQa: (candidateDataset?: DatasetCreateResponse | null) => Promise<QaComparisonResult | null>
   onOpenMapSelectionQualityEvidence: () => void
   onRunOrthophotoAnalysis: (bbox: VectorSelectionBBox) => Promise<boolean>
+  onSelectOrthophotoProduct: (productKey: string) => void
   onClearQualityEvidence?: () => void
 }
 
@@ -534,6 +539,10 @@ export function MapWorkspace({
   orthophotoAnalysisRunning,
   orthophotoAnalysisQuality,
   orthophotoAnalysisDetectionCount,
+  orthophotoProducts,
+  selectedOrthophotoProductKey,
+  orthophotoResult,
+  orthophotoImageUrl,
   availableMapDatasets,
   selectedMapDatasetId,
   onSelectMapArea,
@@ -554,6 +563,7 @@ export function MapWorkspace({
   onRunMapSelectionQa,
   onOpenMapSelectionQualityEvidence,
   onRunOrthophotoAnalysis,
+  onSelectOrthophotoProduct,
   onClearQualityEvidence,
 }: MapWorkspaceProps): JSX.Element {
   const [advancedMode, setAdvancedMode] = useState(false)
@@ -615,6 +625,15 @@ export function MapWorkspace({
   const activeTheme = DATA_THEMES.find((theme) => theme.id === activeThemeId) ?? DATA_THEMES[0]
   const activeThemeMapStyle = DATA_THEME_MAP_STYLES[activeTheme.id]
   const analysisOverlayActive = mapContentMode === 'analysis' && analysisLayerAvailable && Boolean(mapFeatureCollection)
+  const selectedOrthophotoProduct = orthophotoProducts.find((item) => item.key === selectedOrthophotoProductKey) ?? null
+  const orthophotoImageOverlay = orthophotoResult && orthophotoImageUrl && orthophotoResult.bbox_epsg4326.length === 4
+    ? {
+        url: orthophotoImageUrl,
+        bbox: orthophotoResult.bbox_epsg4326 as [number, number, number, number],
+        label: orthophotoResult.display_name,
+        opacity: 0.9,
+      }
+    : null
   const activeThemeDataset = themeDatasetMap[activeTheme.id]
   const activeScopeProject = projects.find((project) => project.id === selectedProjectId) ?? null
   const activeScopeLabel = activeScopeProject ? operationalScopeProjectLabel(activeScopeProject) : 'Werkgebied'
@@ -1173,6 +1192,7 @@ export function MapWorkspace({
                 areaData={areaFeatureCollection}
                 selectedFeature={selectedFeature}
                 selectionData={analysisMode === 'current' ? mapSelectionResult?.geojson ?? null : null}
+                imageOverlay={orthophotoImageOverlay}
                 selectionBbox={mapSelectionBbox}
                 bboxSelectionMode={bboxSelectionMode}
                 visible={mapLayerVisible}
@@ -1188,6 +1208,7 @@ export function MapWorkspace({
               />
               <div className="geo-map-legend" aria-label="Kaartlegende">
                 <span><i className="geo-legend-area" /> Werkgebied</span>
+                {orthophotoImageOverlay ? <span><i className="geo-legend-imagery" /> {orthophotoImageOverlay.label}</span> : null}
                 {analysisOverlayActive ? (
                   <>
                     <span><i className="geo-legend-layer geo-legend-layer-buildings" /> AI-kandidaten</span>
@@ -1233,9 +1254,25 @@ export function MapWorkspace({
               <div className={`geo-image-analysis geo-image-analysis-${orthophotoAnalysisStage}`}>
                 <div>
                   <span>Beeldanalyse</span>
-                  <strong>Gebouwen herkennen op luchtbeeld</strong>
-                  <small>Officieel luchtbeeld, lokaal AI-model en automatische controle met GRB.</small>
+                  <strong>{selectedOrthophotoProduct?.supports_detection ? 'Gebouwen herkennen op luchtbeeld' : 'Historisch luchtbeeld bekijken'}</strong>
+                  <small>
+                    {selectedOrthophotoProduct?.supports_detection
+                      ? 'Officieel luchtbeeld, lokaal AI-model en automatische controle met GRB.'
+                      : 'Officieel historisch mozaïek. Geen vergelijking met de actuele GRB-toestand.'}
+                  </small>
                 </div>
+                <label className="geo-orthophoto-product">
+                  <span>Luchtbeeld</span>
+                  <select
+                    value={selectedOrthophotoProductKey}
+                    onChange={(event) => onSelectOrthophotoProduct(event.target.value)}
+                    disabled={orthophotoAnalysisRunning}
+                  >
+                    {orthophotoProducts.map((product) => (
+                      <option key={product.key} value={product.key}>{product.display_name}</option>
+                    ))}
+                  </select>
+                </label>
                 <button
                   className="primary-action"
                   disabled={orthophotoAnalysisRunning}
@@ -1249,8 +1286,8 @@ export function MapWorkspace({
                       : orthophotoAnalysisStage === 'validating'
                         ? 'Controleren...'
                         : orthophotoAnalysisStage === 'complete'
-                          ? 'Opnieuw analyseren'
-                          : 'Herken gebouwen'}
+                          ? selectedOrthophotoProduct?.supports_detection ? 'Opnieuw analyseren' : 'Opnieuw tonen'
+                          : selectedOrthophotoProduct?.supports_detection ? 'Herken gebouwen' : 'Toon luchtbeeld'}
                 </button>
                 {orthophotoAnalysisStatus ? <p role="status">{orthophotoAnalysisStatus}</p> : null}
                 {orthophotoAnalysisError ? <p className="error" role="alert">{orthophotoAnalysisError}</p> : null}

frontend/src/hooks/useMapOrthophotoAnalysis.ts

@@ -5,6 +5,7 @@ import type {
   DetectionQaResult,
   DetectionRunResponse,
   OrthophotoAcquisitionResult,
+  OrthophotoProductRead,
   VectorSelectionBBox,
 } from '../types'
 import { formatError } from '../lib/formatError'
@@ -87,6 +88,33 @@ export function useMapOrthophotoAnalysis({
   const [lastQuality, setLastQuality] = useState<DetectionQaResult | null>(null)
   const [lastDetectionCount, setLastDetectionCount] = useState<number | null>(null)
   const [lastAnalysisRunId, setLastAnalysisRunId] = useState<string | null>(null)
+  const [products, setProducts] = useState<OrthophotoProductRead[]>([])
+  const [selectedProductKey, setSelectedProductKey] = useState('most_recent')
+
+  useEffect(() => {
+    if (!selectedProjectId) {
+      setProducts([])
+      return
+    }
+    let cancelled = false
+    void datasetsApi.listOrthophotoProducts(selectedProjectId)
+      .then((response) => {
+        if (!cancelled) {
+          setProducts(response.items)
+          setSelectedProductKey((current) =>
+            response.items.some((item) => item.key === current) ? current : response.items[0]?.key ?? 'most_recent',
+          )
+        }
+      })
+      .catch((caught) => {
+        if (!cancelled) {
+          setError(formatError(caught, 'De luchtbeeldcatalogus kon niet worden geladen'))
+        }
+      })
+    return () => {
+      cancelled = true
+    }
+  }, [selectedProjectId])
 
   useEffect(() => {
     setStage('idle')
@@ -96,7 +124,7 @@ export function useMapOrthophotoAnalysis({
     setLastQuality(null)
     setLastDetectionCount(null)
     setLastAnalysisRunId(null)
-  }, [selectionBbox?.min_x, selectionBbox?.min_y, selectionBbox?.max_x, selectionBbox?.max_y])
+  }, [selectionBbox?.min_x, selectionBbox?.min_y, selectionBbox?.max_x, selectionBbox?.max_y, selectedProductKey])
 
   const run = async (bbox: VectorSelectionBBox): Promise<boolean> => {
     if (!selectedProjectId) {
@@ -115,6 +143,7 @@ export function useMapOrthophotoAnalysis({
       const job = await datasetsApi.acquireOrthophoto(selectedProjectId, {
         bbox,
         area_id: selectedAreaId || undefined,
+        product_key: selectedProductKey,
       })
       const acquisition = job.result_json as unknown as OrthophotoAcquisitionResult | null
       const datasetId = job.output_dataset_id || acquisition?.output_dataset_id
@@ -124,6 +153,14 @@ export function useMapOrthophotoAnalysis({
       setLastResult(acquisition)
       await loadProjectData(selectedProjectId)
 
+      if (!acquisition.supports_detection) {
+        setStatus(
+          `${acquisition.display_name} is ingeladen en op de kaart geplaatst. ${acquisition.limitation_message}`,
+        )
+        setStage('complete')
+        return true
+      }
+
       setStage('detecting')
       setStatus('2/3 Lokaal AI-model herkent gebouwen...')
       const detection = await prepareAndRunDetection(datasetId)
@@ -169,9 +206,15 @@ export function useMapOrthophotoAnalysis({
     status,
     error,
     lastResult,
+    imageUrl: lastResult && selectedProjectId
+      ? datasetsApi.orthophotoImageUrl(selectedProjectId, lastResult.output_dataset_id)
+      : null,
     lastQuality,
     lastDetectionCount,
     lastAnalysisRunId,
+    products,
+    selectedProductKey,
+    setSelectedProductKey,
     running: stage === 'acquiring' || stage === 'detecting' || stage === 'validating',
     run,
   }

frontend/src/services/api/datasets.ts

@@ -17,6 +17,7 @@ import type {
   RasterNdwiRequest,
   RasterNdbiRequest,
   OrthophotoAcquireRequest,
+  OrthophotoProductRead,
 } from '../../types'
 
 export const datasetsApi = {
@@ -84,6 +85,10 @@ export const datasetsApi = {
   },
   acquireOrthophoto: (projectId: string, payload: OrthophotoAcquireRequest): Promise<JobRead> =>
     apiPost<JobRead>(`/api/v1/projects/${projectId}/datasets/orthophoto/acquire`, payload),
+  listOrthophotoProducts: (projectId: string): Promise<{ items: OrthophotoProductRead[]; total: number }> =>
+    apiGet<{ items: OrthophotoProductRead[]; total: number }>(`/api/v1/projects/${projectId}/datasets/orthophoto/products`),
+  orthophotoImageUrl: (projectId: string, datasetId: string): string =>
+    `/api/v1/projects/${projectId}/datasets/${datasetId}/raster/image`,
   refreshMetadata: (projectId: string, datasetId: string): Promise<DatasetCreateResponse> =>
     apiPost<DatasetCreateResponse>(`/api/v1/projects/${projectId}/datasets/${datasetId}/metadata/refresh`, {}),
   inspectRaster: (projectId: string, datasetId: string): Promise<RasterInspectResponse> =>

frontend/src/styles/app.css

@@ -5953,6 +5953,11 @@ section {
   background: rgba(107, 74, 170, 0.18);
 }
 
+.geo-map-legend .geo-legend-imagery {
+  border-color: #334155;
+  background: linear-gradient(135deg, #7a9b68 0 33%, #d1b37a 33% 66%, #8eb5cb 66%);
+}
+
 .geo-map-legend .geo-legend-added {
   border-color: #15803d;
   background: rgba(22, 163, 74, 0.2);
@@ -6085,6 +6090,23 @@ section {
   min-height: 2.35rem;
 }
 
+.geo-orthophoto-product {
+  display: grid;
+  gap: 0.3rem;
+}
+
+.geo-orthophoto-product select {
+  width: 100%;
+  min-height: 2.25rem;
+  border: 1px solid #b9cec7;
+  border-radius: 4px;
+  padding: 0.35rem 0.5rem;
+  background: #fff;
+  color: #263a34;
+  font: inherit;
+  font-size: 0.72rem;
+}
+
 .geo-image-analysis-acquiring,
 .geo-image-analysis-detecting,
 .geo-image-analysis-validating {
@@ -6675,6 +6697,28 @@ section {
   grid-column: 1 / -1;
 }
 
+.source-catalog-loaded {
+  display: grid;
+  grid-template-columns: repeat(2, minmax(0, 1fr));
+  gap: 0.6rem;
+  margin-top: 0.75rem;
+}
+
+.source-catalog-loaded article {
+  display: grid;
+  gap: 0.2rem;
+  border-left: 3px solid #287051;
+  padding: 0.6rem 0.75rem;
+  background: #f2f8f5;
+}
+
+.source-catalog-loaded span,
+.source-catalog-loaded p {
+  margin: 0;
+  color: #5d6d67;
+  font-size: 0.72rem;
+}
+
 .source-opportunity-list {
   margin-top: 0.8rem;
   border-top: 1px solid #dfe7e4;
@@ -6716,6 +6760,7 @@ section {
 
 @media (max-width: 980px) {
   .source-catalog-grid,
+  .source-catalog-loaded,
   .source-opportunity-list > div {
     grid-template-columns: repeat(2, minmax(0, 1fr));
   }
@@ -6723,6 +6768,7 @@ section {
 
 @media (max-width: 640px) {
   .source-catalog-grid,
+  .source-catalog-loaded,
   .source-opportunity-list > div {
     grid-template-columns: 1fr;
   }

frontend/src/types.ts

@@ -301,13 +301,31 @@ export interface VectorSelectionBBox {
 export interface OrthophotoAcquireRequest {
   bbox: VectorSelectionBBox
   area_id?: string
+  product_key?: string
   force_refresh?: boolean
 }
 
+export interface OrthophotoProductRead {
+  key: string
+  display_name: string
+  observation_label: string
+  temporal_granularity: string
+  native_resolution_m: number
+  supports_detection: boolean
+  color_mode: string
+  catalog_url: string
+  limitation_message: string
+}
+
 export interface OrthophotoAcquisitionResult {
   output_dataset_id: string
   reused: boolean
   provider: string
+  product_key: string
+  display_name: string
+  observation_label: string
+  temporal_granularity: string
+  supports_detection: boolean
   layer: string
   width: number
   height: number
@@ -318,6 +336,13 @@ export interface OrthophotoAcquisitionResult {
   limitation_message: string
 }
 
+export interface MapImageOverlay {
+  url: string
+  bbox: [number, number, number, number]
+  label: string
+  opacity?: number
+}
+
 export interface MapViewportState {
   bbox: VectorSelectionBBox
   zoom: number

scripts/README.md

@@ -1419,6 +1419,25 @@ artifacts without persistence using `--fetch-only`; bound a run with
 boundaries as resumable request partitions, preserve the native 10 m
 resolution and merge locally before exact clipping to the regional union.
 
+## Waterinfo station histories
+
+Provision real annual station observations for the persisted Mol Area:
+
+```bash
+docker exec geointel python /app/scripts/provision_waterinfo_station_history.py \
+  --project-name "Kempen Regional Workbench" \
+  --area-name "Gemeente Mol" \
+  --parameters water_level,discharge \
+  --from-year 2013 --to-year 2025
+```
+
+Use `--fetch-only` before first persistence or `--force` to refresh retained
+source JSON. The operator is idempotent for existing station/year Datasets,
+retains source checksums and refuses station sets above `--max-stations`. A
+missing discharge series is reported without synthesizing values. Different
+stations remain separate temporal series and may not be treated as area-wide
+water level or volume.
+
 ## Tower deployment
 
 Push the local branch to Gitea, then rebuild the Unraid/Tower Docker runtime:

scripts/audit_api_contracts.py

@@ -12,6 +12,7 @@ DOCS = ROOT / "docs" / "API_CONTRACTS.md"  # docs/API_CONTRACTS.md
 ALLOWED_NON_ENVELOPE_ENDPOINTS = {
     ("GET", "/health"),
     ("GET", "/api/v1/exports/{export_id}/download"),
+    ("GET", "/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/image"),
 }
 
 IGNORED_OPENAPI_PATHS = {

scripts/provision_waterinfo_station_history.py

@@ -0,0 +1,579 @@
+"""Provision official Waterinfo station time series for a persisted Area.
+
+This operator discovers annual station series through the public Waterinfo
+KiWIS service, filters station points against the exact persisted Area and
+imports one immutable GeoJSON point dataset per station and observation year.
+It never averages different stations and never interprets a point measurement
+as area-wide water volume.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import math
+import os
+import re
+import sys
+from dataclasses import dataclass
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any
+
+import requests
+from requests.adapters import HTTPAdapter
+from shapely.geometry import Point, mapping, shape
+from urllib3.util.retry import Retry
+
+
+DEFAULT_API_URL = "http://127.0.0.1:8000"
+DEFAULT_PROJECT_NAME = "Kempen Regional Workbench"
+DEFAULT_AREA_NAME = "Gemeente Mol"
+DEFAULT_OUTPUT_DIR = Path("/app/storage/operator-data/waterinfo/mol")
+KIWIS_URL = "https://download.waterinfo.be/tsmdownload/KiWIS/KiWIS"
+CATALOG_URL = "https://waterinfo.vlaanderen.be/"
+ATTRIBUTION = "Bron: Waterinfo Vlaanderen / Vlaamse Milieumaatschappij"
+
+
+@dataclass(frozen=True)
+class ParameterDefinition:
+    key: str
+    group_id: str
+    property_name: str
+    label: str
+    unit: str
+    reference_layer_name: str
+    limitation: str
+
+
+PARAMETERS = {
+    "water_level": ParameterDefinition(
+        key="water_level",
+        group_id="192784",
+        property_name="annual_mean_water_level_m",
+        label="Jaargemiddelde waterstand",
+        unit="m",
+        reference_layer_name="water_level_station",
+        limitation=(
+            "Dit is een jaargemiddelde op één meetstation. De waarde geldt niet voor het volledige geselecteerde gebied "
+            "en levert zonder profiel- of bathymetriegegevens geen watervolume op."
+        ),
+    ),
+    "discharge": ParameterDefinition(
+        key="discharge",
+        group_id="192895",
+        property_name="annual_mean_discharge_m3_s",
+        label="Jaargemiddeld debiet",
+        unit="m³/s",
+        reference_layer_name="discharge_station",
+        limitation=(
+            "Dit is een jaargemiddeld debiet op één meetstation. De waarde geldt niet voor alle waterlopen in het "
+            "geselecteerde gebied en is geen gebiedsdekkend watervolume."
+        ),
+    ),
+}
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Provision official annual Waterinfo station histories.")
+    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL))
+    parser.add_argument("--project-name", default=DEFAULT_PROJECT_NAME)
+    parser.add_argument("--area-name", default=DEFAULT_AREA_NAME)
+    parser.add_argument("--parameters", default="water_level,discharge")
+    parser.add_argument("--from-year", type=int, default=2013)
+    parser.add_argument("--to-year", type=int, default=datetime.now(timezone.utc).year)
+    parser.add_argument("--min-observations", type=int, default=2)
+    parser.add_argument("--max-stations", type=int, default=50)
+    parser.add_argument("--output-dir", type=Path, default=Path(os.environ.get("WATERINFO_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)))
+    parser.add_argument("--request-timeout", type=int, default=120)
+    parser.add_argument("--import-timeout", type=int, default=300)
+    parser.add_argument("--force", action="store_true", help="Refetch source artifacts; persisted datasets remain immutable.")
+    parser.add_argument("--fetch-only", action="store_true")
+    return parser.parse_args()
+
+
+def build_session() -> requests.Session:
+    retry = Retry(
+        total=5,
+        connect=5,
+        read=5,
+        status=5,
+        backoff_factor=1.0,
+        status_forcelist=(429, 500, 502, 503, 504),
+        allowed_methods=frozenset({"GET"}),
+        raise_on_status=True,
+    )
+    session = requests.Session()
+    session.headers.update({"User-Agent": "GeoIntel-Waterinfo-Operator/1.0"})
+    adapter = HTTPAdapter(max_retries=retry)
+    session.mount("https://", adapter)
+    session.mount("http://", adapter)
+    return session
+
+
+def response_data(response: requests.Response) -> Any:
+    try:
+        payload = response.json()
+    except ValueError as exc:
+        raise RuntimeError(f"GeoIntel API returned non-JSON ({response.status_code}): {response.text[:300]}") from exc
+    if not response.ok:
+        raise RuntimeError(f"GeoIntel API failed ({response.status_code}): {json.dumps(payload, ensure_ascii=False)[:800]}")
+    if not isinstance(payload, dict) or "data" not in payload:
+        raise RuntimeError("GeoIntel API response does not use the canonical data envelope")
+    return payload["data"]
+
+
+def list_paginated_items(session: requests.Session, url: str, *, timeout: int) -> list[dict[str, Any]]:
+    items: list[dict[str, Any]] = []
+    offset = 0
+    total: int | None = None
+    while total is None or offset < total:
+        page = response_data(session.get(url, params={"limit": 200, "offset": offset}, timeout=timeout))
+        page_items = page.get("items") if isinstance(page, dict) else None
+        if not isinstance(page_items, list):
+            raise RuntimeError(f"GeoIntel list response for {url} has no items array")
+        if total is None:
+            total = int(page.get("total", len(page_items)))
+        items.extend(page_items)
+        if not page_items:
+            break
+        offset += len(page_items)
+    if total is not None and len(items) != total:
+        raise RuntimeError(f"GeoIntel list response for {url} returned {len(items)} of {total} items")
+    return items
+
+
+def locate_workspace(session: requests.Session, base_url: str, args: argparse.Namespace):
+    projects = list_paginated_items(session, f"{base_url}/api/v1/projects", timeout=args.import_timeout)
+    project = next((item for item in projects if item.get("name") == args.project_name), None)
+    if not project:
+        raise RuntimeError(f"Project {args.project_name!r} is missing")
+    project_id = str(project["id"])
+    areas = list_paginated_items(
+        session,
+        f"{base_url}/api/v1/projects/{project_id}/areas",
+        timeout=args.import_timeout,
+    )
+    fragment = args.area_name.strip().casefold()
+    matches = [item for item in areas if fragment in str(item.get("name") or "").casefold()]
+    if len(matches) != 1 or not isinstance(matches[0].get("geometry"), dict):
+        raise RuntimeError(f"Expected one persisted Area with geometry matching {args.area_name!r}, received {len(matches)}")
+    area = matches[0]
+    datasets = list_paginated_items(
+        session,
+        f"{base_url}/api/v1/projects/{project_id}/datasets",
+        timeout=args.import_timeout,
+    )
+    return project_id, str(area["id"]), shape(area["geometry"]), datasets
+
+
+def fetch_json(session: requests.Session, params: dict[str, Any], *, timeout: int) -> Any:
+    response = session.get(KIWIS_URL, params=params, timeout=timeout)
+    response.raise_for_status()
+    try:
+        return response.json()
+    except ValueError as exc:
+        raise RuntimeError(f"Waterinfo returned non-JSON: {response.text[:300]}") from exc
+
+
+def discover_station_series(
+    session: requests.Session,
+    parameter: ParameterDefinition,
+    area_geometry,
+    *,
+    timeout: int,
+) -> tuple[dict[str, Any], list[dict[str, Any]]]:
+    payload = fetch_json(
+        session,
+        {
+            "service": "kisters",
+            "type": "queryServices",
+            "request": "getTimeseriesValueLayer",
+            "datasource": 1,
+            "format": "geojson",
+            "timeseriesgroup_id": parameter.group_id,
+            "metadata": "true",
+            "md_returnfields": (
+                "custom_attributes,station_id,station_no,station_name,ts_id,ts_name,"
+                "stationparameter_name,ts_unitsymbol,parametertype_name"
+            ),
+            "custattr_returnfields": "dataprovider,dataowner",
+            "invalidValue": -9999,
+            "invalidPeriod": "P2Y",
+        },
+        timeout=timeout,
+    )
+    if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection":
+        raise RuntimeError(f"Waterinfo group {parameter.group_id} did not return a GeoJSON FeatureCollection")
+    selected: list[dict[str, Any]] = []
+    for feature in payload.get("features") or []:
+        geometry = feature.get("geometry") if isinstance(feature, dict) else None
+        properties = feature.get("properties") if isinstance(feature, dict) else None
+        if not isinstance(geometry, dict) or geometry.get("type") != "Point" or not isinstance(properties, dict):
+            continue
+        point = shape(geometry)
+        ts_id = properties.get("ts_id")
+        if point.is_empty or ts_id in (None, "") or not area_geometry.covers(point):
+            continue
+        selected.append({"geometry": geometry, "properties": properties, "ts_id": str(ts_id)})
+    return payload, selected
+
+
+def fetch_annual_values(
+    session: requests.Session,
+    ts_id: str,
+    *,
+    from_year: int,
+    to_year: int,
+    timeout: int,
+) -> tuple[Any, dict[int, float]]:
+    payload = fetch_json(
+        session,
+        {
+            "service": "kisters",
+            "type": "queryServices",
+            "request": "getTimeseriesValues",
+            "datasource": 1,
+            "format": "json",
+            "ts_id": ts_id,
+            "metadata": "true",
+            "md_returnfields": "station_name,station_no,stationparameter_name,ts_id,ts_unitsymbol",
+            "from": f"{from_year}-01-01",
+            "to": f"{to_year}-12-31",
+        },
+        timeout=timeout,
+    )
+    entries = payload if isinstance(payload, list) else [payload]
+    values: dict[int, float] = {}
+    for entry in entries:
+        if not isinstance(entry, dict):
+            continue
+        for row in entry.get("data") or []:
+            if not isinstance(row, list) or len(row) < 2:
+                continue
+            try:
+                year = int(str(row[0])[:4])
+                value = float(row[1])
+            except (TypeError, ValueError):
+                continue
+            if year < from_year or year > to_year or not math.isfinite(value) or value <= -9999:
+                continue
+            if year in values and not math.isclose(values[year], value, rel_tol=0.0, abs_tol=1e-12):
+                raise RuntimeError(f"Waterinfo series {ts_id} contains multiple annual values for {year}")
+            values[year] = value
+    return payload, values
+
+
+def safe_slug(value: str) -> str:
+    normalized = re.sub(r"[^a-z0-9]+", "-", value.casefold()).strip("-")
+    return normalized[:80] or "station"
+
+
+def series_key(parameter: ParameterDefinition, station: dict[str, Any]) -> str:
+    properties = station["properties"]
+    station_identity = str(properties.get("station_no") or properties.get("station_id") or station["ts_id"])
+    return f"waterinfo:{parameter.key}:annual:{safe_slug(station_identity)}"
+
+
+def build_snapshot(
+    parameter: ParameterDefinition,
+    station: dict[str, Any],
+    year: int,
+    value: float,
+) -> dict[str, Any]:
+    properties = station["properties"]
+    station_name = str(properties.get("station_name") or "Waterinfo meetstation")
+    station_no = str(properties.get("station_no") or properties.get("station_id") or station["ts_id"])
+    feature_id = f"waterinfo-{parameter.key}-{safe_slug(station_no)}-{year}"
+    return {
+        "type": "FeatureCollection",
+        "name": f"{parameter.label} - {station_name} - {year}",
+        "crs": {"type": "name", "properties": {"name": "EPSG:4326"}},
+        "features": [
+            {
+                "type": "Feature",
+                "id": feature_id,
+                "geometry": station["geometry"],
+                "properties": {
+                    "source_feature_id": feature_id,
+                    "source_name": "waterinfo",
+                    "theme": "water",
+                    "measurement_type": parameter.key,
+                    "observation_year": year,
+                    parameter.property_name: value,
+                    "station_id": properties.get("station_id"),
+                    "station_no": properties.get("station_no"),
+                    "station_name": station_name,
+                    "timeseries_id": station["ts_id"],
+                    "timeseries_name": properties.get("ts_name"),
+                    "reported_unit": properties.get("ts_unitsymbol"),
+                    "data_provider": properties.get("dataprovider"),
+                    "data_owner": properties.get("dataowner"),
+                    "authority_level": "authoritative",
+                    "attribution": ATTRIBUTION,
+                },
+            }
+        ],
+    }
+
+
+def sha256_bytes(content: bytes) -> str:
+    return hashlib.sha256(content).hexdigest()
+
+
+def write_json(path: Path, payload: Any, *, pretty: bool = False) -> str:
+    content = json.dumps(
+        payload,
+        ensure_ascii=False,
+        indent=2 if pretty else None,
+        separators=None if pretty else (",", ":"),
+        sort_keys=pretty,
+    ).encode("utf-8")
+    temporary = path.with_suffix(f"{path.suffix}.partial")
+    temporary.write_bytes(content)
+    temporary.replace(path)
+    return sha256_bytes(content)
+
+
+def upload_snapshot(
+    session: requests.Session,
+    *,
+    base_url: str,
+    project_id: str,
+    area_id: str,
+    parameter: ParameterDefinition,
+    station: dict[str, Any],
+    year: int,
+    path: Path,
+    raw_sha256: str,
+    coverage: tuple[int, int, int],
+    timeout: int,
+) -> dict[str, Any]:
+    properties = station["properties"]
+    station_name = str(properties.get("station_name") or "Waterinfo meetstation")
+    key = series_key(parameter, station)
+    observed_at = f"{year}-01-01T00:00:00Z"
+    source_metadata = {
+        "provider": "Waterinfo Vlaanderen / Vlaamse Milieumaatschappij",
+        "authority_level": "authoritative",
+        "theme": "water",
+        "semantic_metrics": False,
+        "geometry_clipped_to_area": True,
+        "measurement_type": parameter.key,
+        "station_name": station_name,
+        "station_no": properties.get("station_no"),
+        "timeseries_id": station["ts_id"],
+        "attribution": ATTRIBUTION,
+        "catalog_url": CATALOG_URL,
+        "temporal_series_label": f"{station_name} - {parameter.label.lower()}",
+        "observation_date_precision": "year",
+        "selection_aggregation": {
+            "metric_key": parameter.key,
+            "method": "mean",
+            "property": parameter.property_name,
+            "label": parameter.label,
+            "unit": parameter.unit,
+            "warning": parameter.limitation,
+        },
+    }
+    provenance_metadata = {
+        "operator_tool": "provision_waterinfo_station_history.py",
+        "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
+        "kiwis_url": KIWIS_URL,
+        "timeseries_group_id": parameter.group_id,
+        "timeseries_id": station["ts_id"],
+        "station_id": properties.get("station_id"),
+        "station_no": properties.get("station_no"),
+        "raw_timeseries_sha256": raw_sha256,
+        "coverage_first_year": coverage[0],
+        "coverage_last_year": coverage[1],
+        "coverage_observation_count": coverage[2],
+        "generated_at": datetime.now(timezone.utc).isoformat(),
+        "limitation_message": parameter.limitation,
+    }
+    with path.open("rb") as handle:
+        response = session.post(
+            f"{base_url}/api/v1/projects/{project_id}/datasets/upload",
+            data={
+                "dataset_type": "vector",
+                "source": "operator_official_import",
+                "dataset_role": "reference",
+                "source_name": "waterinfo",
+                "reference_layer_name": parameter.reference_layer_name,
+                "source_metadata_json": json.dumps(source_metadata, ensure_ascii=False),
+                "provenance_metadata_json": json.dumps(provenance_metadata, ensure_ascii=False),
+                "area_id": area_id,
+                "temporal_series_key": key,
+                "observed_at": observed_at,
+                "valid_from": observed_at,
+                "valid_to": f"{year}-12-31T23:59:59Z",
+                "temporal_granularity": "year",
+                "source_version": f"{station['ts_id']}:{year}",
+            },
+            files={"file": (path.name, handle, "application/geo+json")},
+            timeout=timeout,
+        )
+    return response_data(response)
+
+
+def main() -> int:
+    args = parse_args()
+    if args.from_year > args.to_year or args.min_observations < 2 or args.max_stations < 1:
+        print(json.dumps({"status": "error", "message": "Invalid year, observation or station limits"}), file=sys.stderr)
+        return 2
+    requested_keys = [value.strip() for value in args.parameters.split(",") if value.strip()]
+    unsupported = [key for key in requested_keys if key not in PARAMETERS]
+    if unsupported or not requested_keys:
+        print(json.dumps({"status": "error", "message": f"Unsupported parameters: {unsupported}"}), file=sys.stderr)
+        return 2
+
+    args.output_dir.mkdir(parents=True, exist_ok=True)
+    base_url = args.base_url.rstrip("/")
+    results: list[dict[str, Any]] = []
+    try:
+        with requests.Session() as api_session:
+            project_id, area_id, area_geometry, existing = locate_workspace(api_session, base_url, args)
+            with build_session() as source_session:
+                for parameter_key in requested_keys:
+                    parameter = PARAMETERS[parameter_key]
+                    station_layer, stations = discover_station_series(
+                        source_session,
+                        parameter,
+                        area_geometry,
+                        timeout=args.request_timeout,
+                    )
+                    layer_path = args.output_dir / f"waterinfo_{parameter.key}_station_layer.json"
+                    layer_sha256 = write_json(layer_path, station_layer)
+                    if len(stations) > args.max_stations:
+                        raise RuntimeError(
+                            f"Waterinfo returned {len(stations)} in-area {parameter.key} stations; safety limit is {args.max_stations}"
+                        )
+                    for station in stations:
+                        raw_path = args.output_dir / f"waterinfo_{parameter.key}_{station['ts_id']}_annual.json"
+                        if args.force or not raw_path.exists():
+                            raw_payload, values = fetch_annual_values(
+                                source_session,
+                                station["ts_id"],
+                                from_year=args.from_year,
+                                to_year=args.to_year,
+                                timeout=args.request_timeout,
+                            )
+                            raw_sha256 = write_json(raw_path, raw_payload)
+                        else:
+                            raw_payload = json.loads(raw_path.read_text(encoding="utf-8"))
+                            raw_sha256 = sha256_bytes(raw_path.read_bytes())
+                            entries = raw_payload if isinstance(raw_payload, list) else [raw_payload]
+                            values = {}
+                            for entry in entries:
+                                for row in entry.get("data", []) if isinstance(entry, dict) else []:
+                                    try:
+                                        year = int(str(row[0])[:4])
+                                        value = float(row[1])
+                                    except (IndexError, TypeError, ValueError):
+                                        continue
+                                    if args.from_year <= year <= args.to_year and math.isfinite(value) and value > -9999:
+                                        values[year] = value
+                        if len(values) < args.min_observations:
+                            results.append(
+                                {
+                                    "parameter": parameter.key,
+                                    "timeseries_id": station["ts_id"],
+                                    "status": "insufficient_observations",
+                                    "observation_count": len(values),
+                                }
+                            )
+                            continue
+                        coverage = (min(values), max(values), len(values))
+                        key = series_key(parameter, station)
+                        for year, value in sorted(values.items()):
+                            path = args.output_dir / f"{safe_slug(key)}_{year}.geojson"
+                            snapshot = build_snapshot(parameter, station, year, value)
+                            write_json(path, snapshot)
+                            existing_dataset = next(
+                                (
+                                    item
+                                    for item in existing
+                                    if item.get("temporal_series_key") == key
+                                    and str(item.get("observed_at") or "").startswith(str(year))
+                                ),
+                                None,
+                            )
+                            if existing_dataset:
+                                results.append(
+                                    {
+                                        "parameter": parameter.key,
+                                        "timeseries_id": station["ts_id"],
+                                        "year": year,
+                                        "dataset_id": existing_dataset["id"],
+                                        "status": "existing",
+                                    }
+                                )
+                            elif args.fetch_only:
+                                results.append(
+                                    {
+                                        "parameter": parameter.key,
+                                        "timeseries_id": station["ts_id"],
+                                        "year": year,
+                                        "path": str(path),
+                                        "status": "prepared",
+                                    }
+                                )
+                            else:
+                                dataset = upload_snapshot(
+                                    api_session,
+                                    base_url=base_url,
+                                    project_id=project_id,
+                                    area_id=area_id,
+                                    parameter=parameter,
+                                    station=station,
+                                    year=year,
+                                    path=path,
+                                    raw_sha256=raw_sha256,
+                                    coverage=coverage,
+                                    timeout=args.import_timeout,
+                                )
+                                results.append(
+                                    {
+                                        "parameter": parameter.key,
+                                        "timeseries_id": station["ts_id"],
+                                        "year": year,
+                                        "dataset_id": dataset["id"],
+                                        "status": "imported",
+                                    }
+                                )
+                    manifest = {
+                        "schema_version": 1,
+                        "parameter": parameter.key,
+                        "timeseries_group_id": parameter.group_id,
+                        "area_name": args.area_name,
+                        "station_count": len(stations),
+                        "station_layer_path": str(layer_path),
+                        "station_layer_sha256": layer_sha256,
+                        "from_year": args.from_year,
+                        "to_year": args.to_year,
+                        "generated_at": datetime.now(timezone.utc).isoformat(),
+                    }
+                    write_json(args.output_dir / f"waterinfo_{parameter.key}_manifest.json", manifest, pretty=True)
+    except (OSError, RuntimeError, requests.RequestException, ValueError, KeyError) as exc:
+        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
+        return 1
+
+    print(
+        json.dumps(
+            {
+                "status": "ok",
+                "project": args.project_name,
+                "area": args.area_name,
+                "results": results,
+            },
+            ensure_ascii=False,
+            indent=2,
+        )
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

scripts/run_readiness_check.sh

@@ -47,6 +47,7 @@ ${PYTHON_BIN} -m py_compile scripts/provision_mol_context_layers.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_population_history.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_historical_landuse.py
 ${PYTHON_BIN} -m py_compile scripts/provision_official_landuse_timeseries.py
+${PYTHON_BIN} -m py_compile scripts/provision_waterinfo_station_history.py
 ${PYTHON_BIN} -m py_compile scripts/provision_regional_timeseries.py
 ${PYTHON_BIN} -m py_compile scripts/geographic_scopes.py
 ${PYTHON_BIN} -m py_compile scripts/provision_geographic_scope.py