feat: add cross-domain Mol data profile

.env.example

@@ -28,6 +28,13 @@ FLOOD_HAZARD_MAX_SIDE_M=20000
 FLOOD_HAZARD_MAX_PIXELS=12000000
 FLOOD_HAZARD_TIMEOUT_SECONDS=300
 FLOOD_HAZARD_MAX_RESPONSE_MB=160
+THEMATIC_RASTER_ENABLED=true
+THEMATIC_RASTER_WCS_URL=https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs
+THEMATIC_RASTER_MIN_SIDE_M=100
+THEMATIC_RASTER_MAX_SIDE_M=20000
+THEMATIC_RASTER_MAX_PIXELS=12000000
+THEMATIC_RASTER_TIMEOUT_SECONDS=300
+THEMATIC_RASTER_MAX_RESPONSE_MB=160
 YOLO_ENABLED=false
 YOLO_MODELS_DIR=/app/models
 YOLO_MODEL_PATH=

CHANGELOG.md

@@ -7,6 +7,30 @@
 
 # Changelog
 
+## Sprint 213-214 Cross-domain area profile (2026-07-16)
+
+- Implemented one allowlisted MercatorNet WCS registry for official Flemish
+  space occupation 2025, open space 2022, population density 2019, public
+  transport node value 2022 and total service level 2022 rasters.
+- Added bounded tiled acquisition, exact Area clipping in EPSG:31370, raster
+  value/unit validation, checksummed Dataset/DatasetVersion provenance and
+  source-correct selection metrics without accepting arbitrary service URLs or
+  coverage identifiers.
+- Added MapLibre image overlays, legends and current-state selection for all
+  five products. Raster cell values are presented as hectares, an explicitly
+  estimated population total/density or source scores, never as object counts.
+- Grounded local Ollama answers in the persisted thematic measurements and
+  retained unsupported-current-count, live-timetable and causal limitations.
+- Added the official DOV digital soil map as an explicit Mol operator. It
+  paginates all bounded `bodemkaart:bodemtypes` features, stores checksummed raw
+  evidence, clips exactly in Lambert 72 and imports through DatasetService.
+- Added a Soil map theme with mapped hectares and inspectable soil type,
+  texture and drainage fields. The 1949-1971 survey period and 1:20,000 scale
+  remain visible; current drainage is never inferred.
+- Added focused acquisition, GIS, analysis, API, AI-context, operator,
+  packaging and frontend-contract tests. No migration or direct database write
+  was introduced.
+
 ## Sprint 212 Platform-wide official source portfolio (2026-07-16)
 
 - Rebalanced the source strategy across six user-facing domains: space and

backend/README.md

@@ -1289,6 +1289,43 @@ Settings: `FLOOD_HAZARD_ENABLED`, `FLOOD_HAZARD_WCS_URL`,
 `FLOOD_HAZARD_MAX_SIDE_M`, `FLOOD_HAZARD_MAX_PIXELS`,
 `FLOOD_HAZARD_TIMEOUT_SECONDS` and `FLOOD_HAZARD_MAX_RESPONSE_MB`.
 
+## Cross-domain thematic rasters and DOV soil
+
+The governed thematic registry exposes five fixed MercatorNet products through
+`GET .../datasets/thematic-raster/products`. Acquisition uses
+`POST .../datasets/thematic-raster/acquire`; selection and PNG rendering use
+`POST .../raster/thematic/select` and `GET .../raster/thematic/image`.
+
+Provision every product for the exact persisted Mol Area:
+
+```bash
+docker exec geointel python /app/scripts/provision_thematic_rasters.py
+```
+
+Inspect the complete 28-municipality matrix without writes, then run it after
+the Mol source/runtime gate passes:
+
+```bash
+docker exec geointel python /app/scripts/provision_thematic_rasters.py \
+  --project-name "Kempen Regional Workbench" --all-municipalities --dry-run
+```
+
+Settings: `THEMATIC_RASTER_ENABLED`, `THEMATIC_RASTER_WCS_URL`,
+`THEMATIC_RASTER_MIN_SIDE_M`, `THEMATIC_RASTER_MAX_SIDE_M`,
+`THEMATIC_RASTER_MAX_PIXELS`, `THEMATIC_RASTER_TIMEOUT_SECONDS` and
+`THEMATIC_RASTER_MAX_RESPONSE_MB`.
+
+Provision the official DOV soil polygons for Mol through the existing vector
+upload path:
+
+```bash
+docker exec geointel python /app/scripts/provision_mol_soil_map.py
+```
+
+Use `--fetch-only` to retain and validate source evidence without importing.
+The operator never writes directly to PostGIS. Soil drainage and related map
+classes represent the 1949-1971 survey and are not current observations.
+
 ## Waterinfo station histories
 
 Run the explicit operator after the regional workspace and Mol Area exist:

backend/app/api/routes/datasets.py

@@ -26,6 +26,8 @@ from app.schemas import (
     TerrainSelectionRequest,
     FloodHazardAcquireRequest,
     FloodHazardSelectionRequest,
+    ThematicRasterAcquireRequest,
+    ThematicRasterSelectionRequest,
     VectorBBoxResponse,
     VectorBufferRequest,
     VectorClipRequest,
@@ -48,6 +50,8 @@ from app.services.dhmv_acquisition_service import DhmvAcquisitionService
 from app.services.terrain_analysis_service import TerrainAnalysisService
 from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
 from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService
+from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService
+from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService
 from app.utils.response import envelope
 
 router = APIRouter(prefix="/projects/{project_id}", tags=["datasets"])
@@ -207,6 +211,30 @@ def list_flood_hazard_products(project_id: UUID, db: Session = Depends(get_db)):
     return envelope({"items": items, "total": len(items)})
 
 
+@router.post("/datasets/thematic-raster/acquire", response_model=dict)
+def acquire_bounded_thematic_raster(
+    project_id: UUID,
+    payload: ThematicRasterAcquireRequest,
+    db: Session = Depends(get_db),
+):
+    job = JobService.run_sync_job(
+        db=db,
+        project_id=project_id,
+        job_type="raster.thematic.acquire",
+        parameters=payload.model_dump(mode="json"),
+        operation=lambda: ThematicRasterAcquisitionService.acquire(db, project_id, payload),
+    )
+    return envelope(job)
+
+
+@router.get("/datasets/thematic-raster/products", response_model=dict)
+def list_thematic_raster_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 = ThematicRasterAcquisitionService.list_products()
+    return envelope({"items": items, "total": len(items)})
+
+
 @router.get("/datasets", response_model=dict)
 def list_datasets(
     project_id: UUID,
@@ -560,6 +588,30 @@ def raster_flood_hazard_image(
     )
 
 
+@router.post("/datasets/{dataset_id}/raster/thematic/select", response_model=dict)
+def raster_thematic_selection(
+    project_id: UUID,
+    dataset_id: UUID,
+    payload: ThematicRasterSelectionRequest,
+    db: Session = Depends(get_db),
+):
+    return envelope(ThematicRasterAnalysisService.analyze(db, project_id, dataset_id, payload))
+
+
+@router.get("/datasets/{dataset_id}/raster/thematic/image")
+def raster_thematic_image(
+    project_id: UUID,
+    dataset_id: UUID,
+    db: Session = Depends(get_db),
+):
+    content = ThematicRasterAnalysisService.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/core/config.py

@@ -53,6 +53,16 @@ class Settings(BaseSettings):
     flood_hazard_max_pixels: int = Field(default=12_000_000, ge=1, validation_alias="FLOOD_HAZARD_MAX_PIXELS")
     flood_hazard_timeout_seconds: int = Field(default=300, ge=1, validation_alias="FLOOD_HAZARD_TIMEOUT_SECONDS")
     flood_hazard_max_response_mb: int = Field(default=160, ge=1, validation_alias="FLOOD_HAZARD_MAX_RESPONSE_MB")
+    thematic_raster_enabled: bool = Field(default=True, validation_alias="THEMATIC_RASTER_ENABLED")
+    thematic_raster_wcs_url: str = Field(
+        default="https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs",
+        validation_alias="THEMATIC_RASTER_WCS_URL",
+    )
+    thematic_raster_min_side_m: float = Field(default=100.0, gt=0, validation_alias="THEMATIC_RASTER_MIN_SIDE_M")
+    thematic_raster_max_side_m: float = Field(default=20_000.0, gt=0, validation_alias="THEMATIC_RASTER_MAX_SIDE_M")
+    thematic_raster_max_pixels: int = Field(default=12_000_000, ge=1, validation_alias="THEMATIC_RASTER_MAX_PIXELS")
+    thematic_raster_timeout_seconds: int = Field(default=300, ge=1, validation_alias="THEMATIC_RASTER_TIMEOUT_SECONDS")
+    thematic_raster_max_response_mb: int = Field(default=160, ge=1, validation_alias="THEMATIC_RASTER_MAX_RESPONSE_MB")
     redis_url: str | None = Field(default=None, validation_alias="REDIS_URL")
     log_level: str = Field(default="INFO", validation_alias="GEOINTEL_LOG_LEVEL")
     database_statement_timeout_ms: int = Field(default=5_000, validation_alias="DATABASE_STATEMENT_TIMEOUT_MS")

backend/app/schemas/__init__.py

@@ -51,6 +51,15 @@ from .flood_hazard import (
     FloodHazardSelectionResponse,
     FloodHazardSelectionSummary,
 )
+from .thematic_raster import (
+    ThematicRasterAcquireRequest,
+    ThematicRasterAcquisitionResult,
+    ThematicRasterMetric,
+    ThematicRasterProductRead,
+    ThematicRasterSelectionRequest,
+    ThematicRasterSelectionResponse,
+    ThematicRasterSelectionSummary,
+)
 from .external import (
     ExternalFetchRequest,
     ExternalFetchResponse,
@@ -167,6 +176,13 @@ __all__ = [
     "FloodHazardSelectionRequest",
     "FloodHazardSelectionResponse",
     "FloodHazardSelectionSummary",
+    "ThematicRasterAcquireRequest",
+    "ThematicRasterAcquisitionResult",
+    "ThematicRasterMetric",
+    "ThematicRasterProductRead",
+    "ThematicRasterSelectionRequest",
+    "ThematicRasterSelectionResponse",
+    "ThematicRasterSelectionSummary",
     "VectorBBoxResponse",
     "VectorClipRequest",
     "VectorBufferRequest",

backend/app/schemas/thematic_raster.py

@@ -0,0 +1,96 @@
+from __future__ import annotations
+
+from uuid import UUID
+
+from pydantic import BaseModel
+
+from .operations import VectorSelectionBBox
+
+
+class ThematicRasterAcquireRequest(BaseModel):
+    bbox: VectorSelectionBBox
+    area_id: UUID | None = None
+    product_key: str
+    force_refresh: bool = False
+
+
+class ThematicRasterProductRead(BaseModel):
+    key: str
+    display_name: str
+    theme: str
+    metric_kind: str
+    coverage_id: str
+    native_resolution_m: float
+    source_crs: str
+    source_value_unit: str
+    observation_year: int
+    source_version: str
+    catalog_url: str
+    attribution: str
+    license_note: str
+    legend_min_label: str
+    legend_max_label: str
+    limitation_message: str
+
+
+class ThematicRasterAcquisitionResult(BaseModel):
+    output_dataset_id: UUID
+    reused: bool
+    provider: str
+    product_key: str
+    display_name: str
+    theme: str
+    metric_kind: str
+    coverage_id: str
+    resolution_m: float
+    width: int
+    height: int
+    valid_pixel_count: int
+    bbox_epsg4326: list[float]
+    bbox_epsg31370: list[float]
+    observation_year: int
+    source_value_unit: str
+    attribution: str
+    limitation_message: str
+
+
+class ThematicRasterSelectionRequest(BaseModel):
+    bbox: VectorSelectionBBox
+    area_id: UUID | None = None
+
+
+class ThematicRasterMetric(BaseModel):
+    metric_key: str
+    metric_label: str
+    metric_value: float
+    metric_unit: str
+    aggregation_method: str
+    derived: bool = True
+    is_estimate: bool = True
+
+
+class ThematicRasterSelectionSummary(BaseModel):
+    metric_label: str
+    metric_value: float
+    metric_unit: str
+    aggregation_method: str
+    primary_metric_key: str
+    metrics: list[ThematicRasterMetric]
+
+
+class ThematicRasterSelectionResponse(BaseModel):
+    dataset_id: UUID
+    product_key: str
+    theme: str
+    metric_kind: str
+    selection_bbox: VectorSelectionBBox
+    selection_area_id: UUID | None = None
+    selected_cell_count: int
+    valid_cell_count: int
+    coverage_ratio: float
+    resolution_m: float
+    observation_year: int
+    summary: ThematicRasterSelectionSummary
+    unsupported_metrics: list[str]
+    limitation_message: str
+    generated_at: str

backend/app/services/geo_assistant_service.py

@@ -22,8 +22,11 @@ from app.schemas.assistant import (
     AssistantTemporalSeries,
 )
 from app.schemas.flood_hazard import FloodHazardSelectionRequest
+from app.schemas.thematic_raster import ThematicRasterSelectionRequest
 from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
 from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService
+from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService
+from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService
 from app.services.vector_feature_service import VectorFeatureService
 
 
@@ -42,9 +45,17 @@ class GeoAssistantService:
     )
     ESTIMATE_TOPIC_TERMS = {
         "population": ("bevolk", "inwoner"),
+        "space_occupation": ("ruimtebeslag",),
+        "open_space": ("open ruimte",),
+        "accessibility": ("bereikbaar", "knooppunt"),
+        "services": ("voorziening",),
     }
     ESTIMATE_TOPIC_LABELS = {
         "population": "bevolkingswaarden",
+        "space_occupation": "ruimtebeslagoppervlakten",
+        "open_space": "openruimte-oppervlakten",
+        "accessibility": "bereikbaarheidsscores",
+        "services": "voorzieningenscores",
     }
 
     @classmethod
@@ -264,6 +275,19 @@ class GeoAssistantService:
             if dataset.dataset_type == "raster" and dataset.source_name == FloodHazardAcquisitionService.PROVIDER
             and (area is None or dataset.area_id is None or dataset.area_id == area.id)
         ]
+        thematic_candidates = [
+            dataset
+            for dataset in datasets
+            if dataset.dataset_type == "raster" and dataset.source_name == ThematicRasterAcquisitionService.PROVIDER
+            and (area is None or dataset.area_id is None or dataset.area_id == area.id)
+        ]
+        thematic_by_product: dict[str, Dataset] = {}
+        for dataset in thematic_candidates:
+            product_key = str((dataset.source_metadata or {}).get("product_key") or "")
+            current = thematic_by_product.get(product_key)
+            if product_key and (current is None or (dataset.imported_at or datetime.min.replace(tzinfo=timezone.utc)) > (current.imported_at or datetime.min.replace(tzinfo=timezone.utc))):
+                thematic_by_product[product_key] = dataset
+        thematic_datasets = list(thematic_by_product.values())
         warnings: list[str] = []
         context_metrics: list[AssistantContextMetric] = []
         source_dataset_ids: list[UUID] = []
@@ -323,6 +347,47 @@ class GeoAssistantService:
                     }
                 )
 
+            for dataset in sorted(thematic_datasets, key=lambda item: str((item.source_metadata or {}).get("product_key") or item.name)):
+                try:
+                    result = ThematicRasterAnalysisService.analyze(
+                        db,
+                        project_id,
+                        dataset.id,
+                        ThematicRasterSelectionRequest(bbox=bbox, area_id=area.id if area is not None else None),
+                        settings=self.settings,
+                    )
+                except AppError as exc:
+                    warnings.append(f"{dataset.name}: {exc.message}")
+                    continue
+                serialized_metrics: list[dict[str, Any]] = []
+                for metric in result["summary"]["metrics"]:
+                    item = AssistantContextMetric(
+                        theme=result["theme"],
+                        label=str(metric["metric_label"]),
+                        value=float(metric["metric_value"]),
+                        unit=str(metric["metric_unit"]),
+                        source=ThematicRasterAcquisitionService.ATTRIBUTION,
+                        dataset_id=dataset.id,
+                        observed_at=dataset.observed_at,
+                        is_estimate=bool(metric.get("is_estimate", True)),
+                    )
+                    context_metrics.append(item)
+                    serialized_metrics.append(item.model_dump(mode="json"))
+                    serialized_metrics[-1]["measurement_quality"] = "resolutiegebonden_bronmeting"
+                source_dataset_ids.append(dataset.id)
+                current_context.append(
+                    {
+                        "dataset_name": dataset.name,
+                        "dataset_id": str(dataset.id),
+                        "theme": result["theme"],
+                        "source": ThematicRasterAcquisitionService.ATTRIBUTION,
+                        "observed_at": dataset.observed_at.isoformat() if dataset.observed_at else None,
+                        "metrics": serialized_metrics,
+                        "unsupported_metrics": result["unsupported_metrics"],
+                        "warning": result["limitation_message"],
+                    }
+                )
+
             for dataset in sorted(
                 flood_hazard_datasets,
                 key=lambda item: str((item.source_metadata or {}).get("product_key") or item.name),
@@ -425,6 +490,7 @@ class GeoAssistantService:
                 "water_volume_available": False,
                 "water_volume_reason": "Geen bathymetrie gekoppeld voor de permanente inhoud van waterlichamen.",
                 "flood_hazard_scenarios_available": bool(flood_hazard_datasets),
+                "thematic_policy_rasters_available": bool(thematic_datasets),
                 "flood_depth_area_integral_is_concurrent_volume": False,
                 "object_counts_are_supporting_metrics": True,
                 "causal_explanations_available": False,

backend/app/services/thematic_raster_acquisition_service.py

@@ -0,0 +1,597 @@
+from __future__ import annotations
+
+import hashlib
+import json
+import math
+import time
+from dataclasses import dataclass
+from datetime import UTC, datetime
+from pathlib import Path
+from typing import Any, Callable
+from urllib.error import HTTPError, URLError
+from urllib.parse import urlencode
+from urllib.request import Request, urlopen
+from uuid import UUID
+
+from geoalchemy2.shape import to_shape
+from pyproj import Transformer
+from shapely.geometry import box, mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.config import Settings, get_settings
+from app.core.errors import AppError
+from app.models import Area, Dataset, Project
+from app.schemas.thematic_raster import (
+    ThematicRasterAcquireRequest,
+    ThematicRasterAcquisitionResult,
+    ThematicRasterProductRead,
+)
+from app.services.dataset_service import DatasetService
+
+
+@dataclass(frozen=True)
+class ThematicRasterProduct:
+    key: str
+    display_name: str
+    theme: str
+    metric_kind: str
+    coverage_id: str
+    native_resolution_m: float
+    source_value_unit: str
+    observation_year: int
+    source_version: str
+    catalog_url: str
+    legend_min_label: str
+    legend_max_label: str
+    limitation_message: str
+
+
+class ThematicRasterAcquisitionService:
+    """Acquire bounded, allowlisted policy rasters from MercatorNet WCS."""
+
+    PROVIDER = "department_omgeving_thematic_raster"
+    SOURCE_CRS = "EPSG:31370"
+    WCS_VERSION = "1.0.0"
+    NODATA = -9999.0
+    WCS_TILE_SIDE_M = 10_000.0
+    WCS_REQUEST_INTERVAL_SECONDS = 0.5
+    ATTRIBUTION = "Bron: Departement Omgeving, MercatorNet"
+    LICENSE_NOTE = "Publieke GDI-Vlaanderen bron; bronvermelding en productspecifieke gebruiksvoorwaarden blijven van toepassing."
+
+    @staticmethod
+    def _products() -> dict[str, ThematicRasterProduct]:
+        products = (
+            ThematicRasterProduct(
+                key="space_occupation_2025",
+                display_name="Ruimtebeslag Vlaanderen 2025",
+                theme="space_occupation",
+                metric_kind="binary_area",
+                coverage_id="lu:lu_ruibes_vlaa_2025_v3",
+                native_resolution_m=10.0,
+                source_value_unit="class_0_1",
+                observation_year=2025,
+                source_version="Toestand 2025 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/ruimtebeslag-vlaanderen-toestand-2025",
+                legend_min_label="Geen ruimtebeslag",
+                legend_max_label="Ruimtebeslag",
+                limitation_message=(
+                    "Binaire 10 m-kaart volgens de beleidsdefinitie van ruimtebeslag. Celgebaseerde oppervlakte is een "
+                    "resolutiegebonden schatting en is niet gelijk aan uitsluitend bebouwde oppervlakte of verharding."
+                ),
+            ),
+            ThematicRasterProduct(
+                key="open_space_2022",
+                display_name="Open ruimte Vlaanderen 2022",
+                theme="open_space",
+                metric_kind="binary_area",
+                coverage_id="lu:lu_openruimte_vlaa_2022_v3",
+                native_resolution_m=10.0,
+                source_value_unit="class_0_1",
+                observation_year=2022,
+                source_version="Toestand 2022 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/open-ruimte-vlaanderen-toestand-2022",
+                legend_min_label="Geen open ruimte",
+                legend_max_label="Open ruimte",
+                limitation_message=(
+                    "Binaire 10 m-beleidskaart afgeleid uit landgebruik, ruimtebeslag en kernen. Open ruimte is niet "
+                    "synoniem met natuur, bos, publieke toegankelijkheid of planologische bestemming."
+                ),
+            ),
+            ThematicRasterProduct(
+                key="population_density_2019",
+                display_name="Inwonersdichtheid per hectare 2019",
+                theme="population",
+                metric_kind="population_density",
+                coverage_id="ni:ni_inw_ha_vlaa_2019",
+                native_resolution_m=100.0,
+                source_value_unit="inhabitants_per_hectare",
+                observation_year=2019,
+                source_version="Toestand 2019",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/inwonersdichtheid-per-ha-vlaanderen-toestand-2019",
+                legend_min_label="0 inwoners/ha",
+                legend_max_label="Hogere dichtheid",
+                limitation_message=(
+                    "Statistische 1 ha-rasterinschatting voor 2019, gecorrigeerd op statistische-sectorbasis. De som "
+                    "binnen een getekende grens is een rasterraming en geen actuele registertelling."
+                ),
+            ),
+            ThematicRasterProduct(
+                key="node_value_2022",
+                display_name="Knooppuntwaarde collectief vervoer 2022",
+                theme="accessibility",
+                metric_kind="index_score",
+                coverage_id="lu:lu_knptw_ha_2022_v3",
+                native_resolution_m=100.0,
+                source_value_unit="source_index_score",
+                observation_year=2022,
+                source_version="Toestand 2022 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/knooppuntwaarde-per-ha-toestand-2022",
+                legend_min_label="Lagere knooppuntwaarde",
+                legend_max_label="Hogere knooppuntwaarde",
+                limitation_message=(
+                    "Bronindex per hectare op basis van collectief-vervoerknooppunten en afstandsverval. De score is "
+                    "geen percentage, reistijd, dienstregeling van vandaag of garantie op bereikbaarheid."
+                ),
+            ),
+            ThematicRasterProduct(
+                key="service_level_2022",
+                display_name="Totaal voorzieningenniveau 2022",
+                theme="services",
+                metric_kind="normalized_score",
+                coverage_id="lu:lu_totvznv_ha_2022_v3",
+                native_resolution_m=100.0,
+                source_value_unit="score_0_1",
+                observation_year=2022,
+                source_version="Toestand 2022 versie 3",
+                catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/totaal-voorzieningenniveau-toestand-2022",
+                legend_min_label="Lager voorzieningenniveau",
+                legend_max_label="Hoger voorzieningenniveau",
+                limitation_message=(
+                    "Genormaliseerde 0-1 nabijheidsscore voor basis-, regionale en metropolitane voorzieningen in "
+                    "referentiejaar 2022. Dit is geen objecttelling, openingsurencontrole of actuele reistijd."
+                ),
+            ),
+        )
+        return {product.key: product for product in products}
+
+    @staticmethod
+    def list_products() -> list[dict[str, Any]]:
+        return [
+            ThematicRasterProductRead(
+                key=product.key,
+                display_name=product.display_name,
+                theme=product.theme,
+                metric_kind=product.metric_kind,
+                coverage_id=product.coverage_id,
+                native_resolution_m=product.native_resolution_m,
+                source_crs=ThematicRasterAcquisitionService.SOURCE_CRS,
+                source_value_unit=product.source_value_unit,
+                observation_year=product.observation_year,
+                source_version=product.source_version,
+                catalog_url=product.catalog_url,
+                attribution=ThematicRasterAcquisitionService.ATTRIBUTION,
+                license_note=ThematicRasterAcquisitionService.LICENSE_NOTE,
+                legend_min_label=product.legend_min_label,
+                legend_max_label=product.legend_max_label,
+                limitation_message=product.limitation_message,
+            ).model_dump()
+            for product in ThematicRasterAcquisitionService._products().values()
+        ]
+
+    @staticmethod
+    def _product(product_key: str) -> ThematicRasterProduct:
+        product = ThematicRasterAcquisitionService._products().get(product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="THEMATIC_RASTER_PRODUCT_NOT_SUPPORTED",
+                message="Select a product from the governed Flemish thematic raster registry",
+                details={"product_key": product_key},
+                status_code=422,
+            )
+        return product
+
+    @staticmethod
+    def _prepared_request(payload: ThematicRasterAcquireRequest, settings: Settings) -> dict[str, Any]:
+        if not settings.thematic_raster_enabled:
+            raise AppError(code="THEMATIC_RASTER_NOT_CONFIGURED", message="Official thematic raster acquisition is disabled", status_code=503)
+        product = ThematicRasterAcquisitionService._product(payload.product_key)
+        values = (payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+        if payload.bbox.crs.upper() != "EPSG:4326":
+            raise AppError(code="INVALID_BBOX_CRS", message="Thematic raster selection requires EPSG:4326", status_code=400)
+        if not all(math.isfinite(value) for value in values) or values[0] >= values[2] or values[1] >= values[3]:
+            raise AppError(code="INVALID_BBOX", message="Thematic raster selection must be a finite non-empty rectangle", status_code=400)
+
+        transformer = Transformer.from_crs("EPSG:4326", ThematicRasterAcquisitionService.SOURCE_CRS, always_xy=True)
+        raw_bounds = transformer.transform_bounds(*values, densify_pts=21)
+        resolution = product.native_resolution_m
+        lambert_bounds = (
+            math.floor(raw_bounds[0] / resolution) * resolution,
+            math.floor(raw_bounds[1] / resolution) * resolution,
+            math.ceil(raw_bounds[2] / resolution) * resolution,
+            math.ceil(raw_bounds[3] / resolution) * resolution,
+        )
+        width_m = lambert_bounds[2] - lambert_bounds[0]
+        height_m = lambert_bounds[3] - lambert_bounds[1]
+        if width_m < settings.thematic_raster_min_side_m or height_m < settings.thematic_raster_min_side_m:
+            raise AppError(
+                code="THEMATIC_RASTER_SELECTION_TOO_SMALL",
+                message=f"Select an area of at least {settings.thematic_raster_min_side_m:g} by {settings.thematic_raster_min_side_m:g} metres",
+                status_code=422,
+            )
+        if width_m > settings.thematic_raster_max_side_m or height_m > settings.thematic_raster_max_side_m:
+            raise AppError(
+                code="THEMATIC_RASTER_SELECTION_TOO_LARGE",
+                message=f"Select an area no larger than {settings.thematic_raster_max_side_m:g} by {settings.thematic_raster_max_side_m:g} metres",
+                details={"width_m": width_m, "height_m": height_m},
+                status_code=422,
+            )
+        width = max(1, round(width_m / resolution))
+        height = max(1, round(height_m / resolution))
+        if width * height > settings.thematic_raster_max_pixels:
+            raise AppError(
+                code="THEMATIC_RASTER_SELECTION_TOO_LARGE",
+                message="Thematic raster selection exceeds the configured cell limit",
+                details={"pixel_count": width * height, "max_pixels": settings.thematic_raster_max_pixels},
+                status_code=422,
+            )
+        request_identity = {
+            "provider": ThematicRasterAcquisitionService.PROVIDER,
+            "coverage_id": product.coverage_id,
+            "bbox_epsg4326": [round(float(value), 8) for value in values],
+            "bbox_epsg31370": [round(float(value), 3) for value in lambert_bounds],
+            "resolution_m": resolution,
+            "area_id": str(payload.area_id) if payload.area_id else None,
+        }
+        request_hash = hashlib.sha256(json.dumps(request_identity, sort_keys=True).encode()).hexdigest()
+        return {
+            **request_identity,
+            "product": product,
+            "request_hash": request_hash,
+            "width": width,
+            "height": height,
+        }
+
+    @staticmethod
+    def _wcs_request_url(settings: Settings, product: ThematicRasterProduct, bounds: tuple[float, float, float, float]) -> str:
+        query = {
+            "SERVICE": "WCS",
+            "VERSION": ThematicRasterAcquisitionService.WCS_VERSION,
+            "REQUEST": "GetCoverage",
+            "COVERAGE": product.coverage_id,
+            "CRS": ThematicRasterAcquisitionService.SOURCE_CRS,
+            "BBOX": ",".join(f"{value:.3f}" for value in bounds),
+            "RESX": f"{product.native_resolution_m:g}",
+            "RESY": f"{product.native_resolution_m:g}",
+            "FORMAT": "image/tiff",
+            "RESPONSE_CRS": ThematicRasterAcquisitionService.SOURCE_CRS,
+        }
+        return f"{settings.thematic_raster_wcs_url}?{urlencode(query)}"
+
+    @staticmethod
+    def _tile_bounds(prepared: dict[str, Any]) -> list[tuple[float, float, float, float]]:
+        min_x, min_y, max_x, max_y = prepared["bbox_epsg31370"]
+        resolution = prepared["product"].native_resolution_m
+        side = max(resolution, math.floor(ThematicRasterAcquisitionService.WCS_TILE_SIDE_M / resolution) * resolution)
+        tiles: list[tuple[float, float, float, float]] = []
+        y = min_y
+        while y < max_y:
+            tile_max_y = min(y + side, max_y)
+            x = min_x
+            while x < max_x:
+                tile_max_x = min(x + side, max_x)
+                tiles.append((x, y, tile_max_x, tile_max_y))
+                x = tile_max_x
+            y = tile_max_y
+        return tiles
+
+    @staticmethod
+    def _scope_geometry(db, project_id: UUID, area_id: UUID | None, bbox_epsg4326: list[float]):
+        if not db.get(Project, project_id):
+            raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
+        selection = box(*bbox_epsg4326)
+        if area_id is None:
+            return selection
+        area = db.get(Area, area_id)
+        if not area:
+            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+        if area.project_id != project_id:
+            raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400)
+        intersection = to_shape(area.geometry).intersection(selection)
+        if intersection.is_empty or intersection.area <= 0:
+            raise AppError(code="THEMATIC_RASTER_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+        return intersection
+
+    @staticmethod
+    def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]:
+        request = Request(request_url, headers={"Accept": "image/tiff,*/*", "User-Agent": "GeoIntel/0.1 bounded-thematic-raster"})
+        max_bytes = settings.thematic_raster_max_response_mb * 1024 * 1024
+        try:
+            with (opener or urlopen)(request, timeout=settings.thematic_raster_timeout_seconds) as response:
+                content_type = str(response.headers.get("Content-Type", ""))
+                content_length = response.headers.get("Content-Length")
+                if content_length and int(content_length) > max_bytes:
+                    raise AppError(code="THEMATIC_RASTER_RESPONSE_TOO_LARGE", message="Official raster response exceeds the configured size limit", status_code=502)
+                content = response.read(max_bytes + 1)
+        except AppError:
+            raise
+        except HTTPError as exc:
+            preview = exc.read(300).decode("utf-8", errors="replace")
+            raise AppError(
+                code="THEMATIC_RASTER_PROVIDER_UNAVAILABLE",
+                message="The official MercatorNet WCS could not complete the bounded request",
+                details={"reason": str(exc), "provider_status_code": int(exc.code), "response_preview": preview},
+                status_code=502,
+            ) from exc
+        except (URLError, TimeoutError, OSError) as exc:
+            raise AppError(
+                code="THEMATIC_RASTER_PROVIDER_UNAVAILABLE",
+                message="The official MercatorNet WCS could not complete the bounded request",
+                details={"reason": str(exc)},
+                status_code=502,
+            ) from exc
+        if len(content) > max_bytes:
+            raise AppError(code="THEMATIC_RASTER_RESPONSE_TOO_LARGE", message="Official raster response exceeds the configured size limit", status_code=502)
+        if not content.startswith((b"II*\x00", b"MM\x00*")):
+            preview = content[:300].decode("utf-8", errors="replace")
+            raise AppError(
+                code="THEMATIC_RASTER_PROVIDER_INVALID_RESPONSE",
+                message="The official MercatorNet service did not return a GeoTIFF coverage",
+                details={"content_type": content_type, "response_preview": preview},
+                status_code=502,
+            )
+        return content, content_type
+
+    @staticmethod
+    def _mosaic(coverages: list[bytes], product: ThematicRasterProduct) -> bytes:
+        if len(coverages) == 1:
+            return coverages[0]
+        try:
+            import rasterio
+            from rasterio.io import MemoryFile
+            from rasterio.merge import merge
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio is required to assemble thematic raster tiles", status_code=503) from exc
+        memories = [MemoryFile(content) for content in coverages]
+        sources = []
+        try:
+            sources = [memory.open() for memory in memories]
+            for source in sources:
+                if source.crs is None or source.crs.to_epsg() != 31370 or source.count != 1:
+                    raise AppError(code="THEMATIC_RASTER_TILE_MISMATCH", message="Thematic raster tiles have incompatible CRS or bands", status_code=502)
+                if not all(math.isclose(abs(float(value)), product.native_resolution_m, abs_tol=0.05) for value in source.res):
+                    raise AppError(code="THEMATIC_RASTER_TILE_MISMATCH", message="Thematic raster tile resolution differs from the registry", status_code=502)
+            mosaic, transform = merge(sources, res=(product.native_resolution_m, product.native_resolution_m), nodata=ThematicRasterAcquisitionService.NODATA, dtype="float32")
+            profile = sources[0].profile.copy()
+            profile.pop("blockxsize", None)
+            profile.pop("blockysize", None)
+            profile.update(driver="GTiff", width=mosaic.shape[2], height=mosaic.shape[1], count=1, dtype="float32", crs=ThematicRasterAcquisitionService.SOURCE_CRS, transform=transform, nodata=ThematicRasterAcquisitionService.NODATA, compress="deflate", predictor=3)
+            with MemoryFile() as output_memory:
+                with output_memory.open(**profile) as output:
+                    output.write(mosaic)
+                return output_memory.read()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(code="THEMATIC_RASTER_TILE_MOSAIC_FAILED", message="Thematic raster tiles could not be assembled", details={"reason": str(exc)}, status_code=502) from exc
+        finally:
+            for source in sources:
+                source.close()
+            for memory in memories:
+                memory.close()
+
+    @staticmethod
+    def _fetch_coverage(prepared: dict[str, Any], settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, dict[str, Any]]:
+        product: ThematicRasterProduct = prepared["product"]
+        request_urls = [ThematicRasterAcquisitionService._wcs_request_url(settings, product, bounds) for bounds in ThematicRasterAcquisitionService._tile_bounds(prepared)]
+        coverages: list[bytes] = []
+        digest = hashlib.sha256()
+        content_types: list[str] = []
+        for index, request_url in enumerate(request_urls):
+            if index and opener is None:
+                time.sleep(ThematicRasterAcquisitionService.WCS_REQUEST_INTERVAL_SECONDS)
+            content, content_type = ThematicRasterAcquisitionService._fetch(request_url, settings, opener)
+            digest.update(len(content).to_bytes(8, "big"))
+            digest.update(content)
+            coverages.append(content)
+            content_types.append(content_type)
+        return ThematicRasterAcquisitionService._mosaic(coverages, product), {
+            "tile_count": len(request_urls),
+            "request_urls": request_urls,
+            "response_content_types": content_types,
+            "coverage_sha256": digest.hexdigest(),
+        }
+
+    @staticmethod
+    def _validate_values(values, product: ThematicRasterProduct) -> None:
+        import numpy as np
+
+        if values.size == 0:
+            raise AppError(code="THEMATIC_RASTER_NO_VALID_DATA", message="The official product contains no valid cells in this selection", status_code=422)
+        minimum = float(values.min())
+        maximum = float(values.max())
+        if minimum < 0:
+            raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Official thematic raster contains unexpected negative values", details={"minimum": minimum}, status_code=502)
+        if product.metric_kind == "binary_area" and not set(np.unique(values).tolist()).issubset({0.0, 1.0}):
+            raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Binary thematic raster contains classes outside 0 and 1", status_code=502)
+        if product.metric_kind == "normalized_score" and maximum > 1.0001:
+            raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Normalized thematic score falls outside the documented 0-1 range", details={"maximum": maximum}, status_code=502)
+
+    @staticmethod
+    def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]:
+        try:
+            import numpy as np
+            from rasterio.io import MemoryFile
+            from rasterio.mask import mask
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for thematic raster validation", status_code=503) from exc
+        product: ThematicRasterProduct = prepared["product"]
+        try:
+            with MemoryFile(content) as source_memory, source_memory.open() as source:
+                if source.crs is None or source.crs.to_epsg() != 31370:
+                    raise AppError(code="THEMATIC_RASTER_INVALID_CRS", message="Official thematic raster must use EPSG:31370", status_code=502)
+                if source.count != 1:
+                    raise AppError(code="THEMATIC_RASTER_INVALID_BANDS", message="Official thematic raster must contain one band", status_code=502)
+                if not all(math.isclose(abs(float(value)), product.native_resolution_m, abs_tol=0.05) for value in source.res):
+                    raise AppError(code="THEMATIC_RASTER_INVALID_RESOLUTION", message="Official thematic raster resolution differs from the registry", status_code=502)
+                transformer = Transformer.from_crs("EPSG:4326", ThematicRasterAcquisitionService.SOURCE_CRS, always_xy=True)
+                scope_metric = shapely_transform(transformer.transform, scope_geometry_4326)
+                clipped, transform = mask(source, [mapping(scope_metric)], crop=True, filled=False, indexes=[1])
+                band = np.ma.asarray(clipped[0], dtype="float32")
+                raw = np.asarray(band.filled(np.nan), dtype="float32")
+                invalid = np.ma.getmaskarray(band) | ~np.isfinite(raw)
+                if source.nodata is not None:
+                    invalid |= np.isclose(raw, float(source.nodata))
+                normalized = np.ma.array(raw, mask=invalid)
+                values = normalized.compressed().astype("float64")
+                ThematicRasterAcquisitionService._validate_values(values, product)
+                profile = source.profile.copy()
+                profile.pop("blockxsize", None)
+                profile.pop("blockysize", None)
+                profile.update(driver="GTiff", width=normalized.shape[1], height=normalized.shape[0], count=1, dtype="float32", crs=ThematicRasterAcquisitionService.SOURCE_CRS, transform=transform, nodata=ThematicRasterAcquisitionService.NODATA, compress="deflate", predictor=3)
+                with MemoryFile() as output_memory:
+                    with output_memory.open(**profile) as output:
+                        output.write(normalized.filled(ThematicRasterAcquisitionService.NODATA), 1)
+                    normalized_content = output_memory.read()
+                return normalized_content, {
+                    "width": int(normalized.shape[1]),
+                    "height": int(normalized.shape[0]),
+                    "valid_pixel_count": int(values.size),
+                    "nodata_value": ThematicRasterAcquisitionService.NODATA,
+                    "resolution_m": product.native_resolution_m,
+                    "minimum_value": float(values.min()),
+                    "maximum_value": float(values.max()),
+                    "p02_value": float(np.percentile(values, 2)),
+                    "p98_value": float(np.percentile(values, 98)),
+                }
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(code="THEMATIC_RASTER_INVALID", message="The official thematic raster could not be validated", details={"reason": str(exc)}, status_code=502) from exc
+
+    @staticmethod
+    def _cached_dataset(db, project_id: UUID, filename: str) -> Dataset | None:
+        candidate = (
+            db.query(Dataset)
+            .filter(Dataset.project_id == project_id, Dataset.name == filename, Dataset.source_name == ThematicRasterAcquisitionService.PROVIDER, Dataset.status == "ready")
+            .order_by(Dataset.imported_at.desc())
+            .first()
+        )
+        return candidate if candidate and candidate.storage_path and Path(candidate.storage_path).is_file() else None
+
+    @staticmethod
+    def acquire(db, project_id: UUID, payload: ThematicRasterAcquireRequest, *, settings: Settings | None = None, opener: Callable[..., Any] | None = None) -> dict[str, Any]:
+        resolved_settings = settings or get_settings()
+        prepared = ThematicRasterAcquisitionService._prepared_request(payload, resolved_settings)
+        product: ThematicRasterProduct = prepared["product"]
+        scope_geometry = ThematicRasterAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
+        filename = f"thematic_{product.key}_{prepared['request_hash'][:12]}.tif"
+        if not payload.force_refresh:
+            cached = ThematicRasterAcquisitionService._cached_dataset(db, project_id, filename)
+            if cached is not None:
+                metadata = cached.source_metadata or {}
+                raster_metadata = cached.metadata_json or {}
+                return ThematicRasterAcquisitionResult(
+                    output_dataset_id=cached.id,
+                    reused=True,
+                    provider=ThematicRasterAcquisitionService.PROVIDER,
+                    product_key=product.key,
+                    display_name=product.display_name,
+                    theme=product.theme,
+                    metric_kind=product.metric_kind,
+                    coverage_id=product.coverage_id,
+                    resolution_m=product.native_resolution_m,
+                    width=int(raster_metadata.get("width", prepared["width"])),
+                    height=int(raster_metadata.get("height", prepared["height"])),
+                    valid_pixel_count=int(metadata.get("valid_pixel_count", 0)),
+                    bbox_epsg4326=prepared["bbox_epsg4326"],
+                    bbox_epsg31370=prepared["bbox_epsg31370"],
+                    observation_year=product.observation_year,
+                    source_value_unit=product.source_value_unit,
+                    attribution=ThematicRasterAcquisitionService.ATTRIBUTION,
+                    limitation_message=product.limitation_message,
+                ).model_dump(mode="json")
+
+        coverage, transfer = ThematicRasterAcquisitionService._fetch_coverage(prepared, resolved_settings, opener)
+        normalized, validation = ThematicRasterAcquisitionService._normalize_raster(coverage, scope_geometry, prepared)
+        acquired_at = datetime.now(UTC)
+        observed_at = datetime(product.observation_year, 12, 31, 23, 59, 59, tzinfo=UTC)
+        scope_key = str(payload.area_id) if payload.area_id else prepared["request_hash"][:24]
+        dataset = DatasetService.import_raster_bytes(
+            db,
+            project_id=project_id,
+            area_id=payload.area_id,
+            filename=filename,
+            content=normalized,
+            source=f"Departement Omgeving MercatorNet WCS {product.coverage_id}",
+            source_name=ThematicRasterAcquisitionService.PROVIDER,
+            temporal_series_key=f"department-omgeving:thematic-raster:{product.key}:{scope_key}",
+            observed_at=observed_at,
+            valid_from=datetime(product.observation_year, 1, 1, tzinfo=UTC),
+            valid_to=observed_at,
+            temporal_granularity="year",
+            source_version=product.source_version,
+            source_metadata={
+                "provider": ThematicRasterAcquisitionService.PROVIDER,
+                "service": "WCS",
+                "service_version": ThematicRasterAcquisitionService.WCS_VERSION,
+                "product_key": product.key,
+                "product_display_name": product.display_name,
+                "theme": product.theme,
+                "metric_kind": product.metric_kind,
+                "coverage_id": product.coverage_id,
+                "native_resolution_m": product.native_resolution_m,
+                "analysis_resolution_m": product.native_resolution_m,
+                "source_crs": ThematicRasterAcquisitionService.SOURCE_CRS,
+                "source_value_unit": product.source_value_unit,
+                "observation_year": product.observation_year,
+                "observation_date_precision": "year",
+                "valid_pixel_count": validation["valid_pixel_count"],
+                "minimum_value": validation["minimum_value"],
+                "maximum_value": validation["maximum_value"],
+                "render_min_value": validation["p02_value"],
+                "render_max_value": validation["p98_value"],
+                "bbox_epsg4326": prepared["bbox_epsg4326"],
+                "bbox_epsg31370": prepared["bbox_epsg31370"],
+                "catalog_url": product.catalog_url,
+                "attribution": ThematicRasterAcquisitionService.ATTRIBUTION,
+                "license_note": ThematicRasterAcquisitionService.LICENSE_NOTE,
+                "legend_min_label": product.legend_min_label,
+                "legend_max_label": product.legend_max_label,
+                "coverage_scope": "municipality" if payload.area_id else "bounded_selection",
+            },
+            provenance_metadata={
+                "acquisition": "explicit_bounded_tiled_wcs_coverage",
+                "acquired_at": acquired_at.isoformat(),
+                "request_hash": prepared["request_hash"],
+                "tile_count": transfer["tile_count"],
+                "tile_request_urls": transfer["request_urls"],
+                "response_content_types": transfer["response_content_types"],
+                "coverage_sha256": transfer["coverage_sha256"],
+                "normalized_sha256": hashlib.sha256(normalized).hexdigest(),
+                "bbox_epsg4326": prepared["bbox_epsg4326"],
+                "bbox_epsg31370": prepared["bbox_epsg31370"],
+                "clipped_to_area_id": str(payload.area_id) if payload.area_id else None,
+                "validation": validation,
+                "limitation_message": product.limitation_message,
+            },
+        )
+        return ThematicRasterAcquisitionResult(
+            output_dataset_id=dataset.id,
+            reused=False,
+            provider=ThematicRasterAcquisitionService.PROVIDER,
+            product_key=product.key,
+            display_name=product.display_name,
+            theme=product.theme,
+            metric_kind=product.metric_kind,
+            coverage_id=product.coverage_id,
+            resolution_m=product.native_resolution_m,
+            width=validation["width"],
+            height=validation["height"],
+            valid_pixel_count=validation["valid_pixel_count"],
+            bbox_epsg4326=prepared["bbox_epsg4326"],
+            bbox_epsg31370=prepared["bbox_epsg31370"],
+            observation_year=product.observation_year,
+            source_value_unit=product.source_value_unit,
+            attribution=ThematicRasterAcquisitionService.ATTRIBUTION,
+            limitation_message=product.limitation_message,
+        ).model_dump(mode="json")

backend/app/services/thematic_raster_analysis_service.py

@@ -0,0 +1,260 @@
+from __future__ import annotations
+
+import io
+import math
+from datetime import UTC, datetime
+from pathlib import Path
+from uuid import UUID
+
+from geoalchemy2.shape import to_shape
+from pyproj import Transformer
+from shapely.geometry import box, mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.config import Settings, get_settings
+from app.core.errors import AppError
+from app.models import Area, Dataset
+from app.schemas.thematic_raster import (
+    ThematicRasterMetric,
+    ThematicRasterSelectionRequest,
+    ThematicRasterSelectionResponse,
+    ThematicRasterSelectionSummary,
+)
+from app.services.thematic_raster_acquisition_service import (
+    ThematicRasterAcquisitionService,
+    ThematicRasterProduct,
+)
+
+
+class ThematicRasterAnalysisService:
+    @staticmethod
+    def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> Dataset:
+        dataset = db.get(Dataset, dataset_id)
+        if not dataset or dataset.project_id != project_id:
+            raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
+        if dataset.dataset_type != "raster" or dataset.source_name != ThematicRasterAcquisitionService.PROVIDER:
+            raise AppError(
+                code="INVALID_THEMATIC_RASTER_DATASET",
+                message="Thematic analysis requires a governed Departement Omgeving raster dataset",
+                status_code=400,
+            )
+        if dataset.status != "ready" or not dataset.storage_path or not Path(dataset.storage_path).is_file():
+            raise AppError(code="DATASET_FILE_MISSING", message="Persisted thematic raster file is unavailable", status_code=404)
+        return dataset
+
+    @staticmethod
+    def _product(dataset: Dataset) -> ThematicRasterProduct:
+        source_metadata = dataset.source_metadata or {}
+        product = ThematicRasterAcquisitionService._products().get(str(source_metadata.get("product_key") or ""))
+        if product is None or source_metadata.get("coverage_id") != product.coverage_id:
+            raise AppError(code="INVALID_THEMATIC_RASTER_METADATA", message="Thematic raster provenance is incomplete", status_code=409)
+        return product
+
+    @staticmethod
+    def _selection_geometry(db, project_id: UUID, payload: ThematicRasterSelectionRequest):
+        selection = box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+        if payload.area_id is None:
+            return selection
+        area = db.get(Area, payload.area_id)
+        if not area:
+            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+        if area.project_id != project_id:
+            raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400)
+        selection = selection.intersection(to_shape(area.geometry))
+        if selection.is_empty or selection.area <= 0:
+            raise AppError(code="THEMATIC_RASTER_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
+        return selection
+
+    @staticmethod
+    def _unsupported_metrics(product: ThematicRasterProduct) -> list[str]:
+        if product.metric_kind == "binary_area":
+            return ["object_count", "parcel_area", "current_land_use"]
+        if product.metric_kind == "population_density":
+            return ["current_population", "household_count", "address_level_population"]
+        if product.metric_kind == "index_score":
+            return ["travel_time_minutes", "current_timetable", "stop_count"]
+        return ["facility_count", "opening_hours", "current_service_availability"]
+
+    @staticmethod
+    def analyze(
+        db,
+        project_id: UUID,
+        dataset_id: UUID,
+        payload: ThematicRasterSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        dataset = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        product = ThematicRasterAnalysisService._product(dataset)
+        selection_4326 = ThematicRasterAnalysisService._selection_geometry(db, project_id, payload)
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.features import geometry_mask
+            from rasterio.mask import mask
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for thematic raster analysis", status_code=503) from exc
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                if source.crs is None or source.crs.to_epsg() != 31370:
+                    raise AppError(code="INVALID_DATASET_CRS", message="Thematic raster CRS must be EPSG:31370", status_code=409)
+                transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True)
+                selection_metric = shapely_transform(transformer.transform, selection_4326)
+                analysis_geometry = selection_metric.intersection(box(*source.bounds))
+                if analysis_geometry.is_empty or analysis_geometry.area <= 0:
+                    raise AppError(code="THEMATIC_RASTER_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted thematic raster", status_code=422)
+                min_x, min_y, max_x, max_y = analysis_geometry.bounds
+                expected_cells = math.ceil((max_x - min_x) / abs(source.res[0])) * math.ceil((max_y - min_y) / abs(source.res[1]))
+                if expected_cells > resolved_settings.thematic_raster_max_pixels:
+                    raise AppError(
+                        code="THEMATIC_RASTER_SELECTION_TOO_LARGE",
+                        message="Thematic raster analysis exceeds the configured cell limit",
+                        details={"pixel_count": expected_cells, "max_pixels": resolved_settings.thematic_raster_max_pixels},
+                        status_code=422,
+                    )
+                clipped, clipped_transform = mask(source, [mapping(analysis_geometry)], crop=True, filled=False, indexes=[1])
+                band = np.ma.asarray(clipped[0], dtype="float64")
+                raw = band.filled(np.nan)
+                selected = geometry_mask([mapping(analysis_geometry)], out_shape=band.shape, transform=clipped_transform, invert=True)
+                valid = selected & ~np.ma.getmaskarray(band) & np.isfinite(raw)
+                if source.nodata is not None:
+                    valid &= ~np.isclose(raw, float(source.nodata))
+                values = raw[valid]
+                ThematicRasterAcquisitionService._validate_values(values, product)
+                selected_cell_count = int(selected.sum())
+                valid_cell_count = int(values.size)
+                resolution_x = abs(float(source.res[0]))
+                resolution_y = abs(float(source.res[1]))
+                cell_area_m2 = resolution_x * resolution_y
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="THEMATIC_RASTER_ANALYSIS_FAILED",
+                message="The persisted thematic raster could not be analysed",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc
+
+        def metric(key: str, label: str, value: float, unit: str, method: str, *, estimate: bool = True) -> ThematicRasterMetric:
+            return ThematicRasterMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+                is_estimate=estimate,
+            )
+
+        if product.metric_kind == "binary_area":
+            positive_count = int(np.count_nonzero(values >= 0.5))
+            positive_area_ha = positive_count * cell_area_m2 / 10_000.0
+            positive_share = positive_count / max(1, valid_cell_count) * 100.0
+            label = "Ruimtebeslag" if product.theme == "space_occupation" else "Open ruimte"
+            metrics = [
+                metric(f"{product.theme}_area_ha", f"{label} in selectie", positive_area_ha, "ha", "positive_source_cells_times_cell_area"),
+                metric(f"{product.theme}_share_pct", f"Aandeel {label.lower()}", positive_share, "%", "positive_source_cells_divided_by_valid_selected_cells"),
+                metric("valid_raster_area_ha", "Rasteroppervlakte met bronwaarde", valid_cell_count * cell_area_m2 / 10_000.0, "ha", "valid_selected_cells_times_cell_area"),
+            ]
+        elif product.metric_kind == "population_density":
+            estimated_population = float(values.sum() * (cell_area_m2 / 10_000.0))
+            metrics = [
+                metric("estimated_inhabitants", "Geraamd aantal inwoners (2019)", estimated_population, "inwoners", "sum_density_times_selected_cell_area_hectares"),
+                metric("population_density_mean_per_ha", "Gemiddelde inwonersdichtheid", values.mean(), "inwoners/ha", "mean_valid_one_hectare_source_cells"),
+                metric("population_density_p90_per_ha", "90e percentiel inwonersdichtheid", np.percentile(values, 90), "inwoners/ha", "percentile_90_valid_source_cells"),
+            ]
+        else:
+            unit = "score" if product.metric_kind == "index_score" else "score (0-1)"
+            label = "Knooppuntwaarde" if product.metric_kind == "index_score" else "Voorzieningenniveau"
+            metrics = [
+                metric(f"{product.theme}_mean", f"Gemiddelde {label.lower()}", values.mean(), unit, "mean_valid_source_cells"),
+                metric(f"{product.theme}_p10", f"10e percentiel {label.lower()}", np.percentile(values, 10), unit, "percentile_10_valid_source_cells"),
+                metric(f"{product.theme}_median", f"Mediaan {label.lower()}", np.percentile(values, 50), unit, "median_valid_source_cells"),
+                metric(f"{product.theme}_p90", f"90e percentiel {label.lower()}", np.percentile(values, 90), unit, "percentile_90_valid_source_cells"),
+            ]
+
+        primary = metrics[0]
+        response = ThematicRasterSelectionResponse(
+            dataset_id=dataset.id,
+            product_key=product.key,
+            theme=product.theme,
+            metric_kind=product.metric_kind,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            selected_cell_count=selected_cell_count,
+            valid_cell_count=valid_cell_count,
+            coverage_ratio=round(valid_cell_count / max(1, selected_cell_count), 6),
+            resolution_m=round(max(resolution_x, resolution_y), 4),
+            observation_year=product.observation_year,
+            summary=ThematicRasterSelectionSummary(
+                metric_label=primary.metric_label,
+                metric_value=primary.metric_value,
+                metric_unit=primary.metric_unit,
+                aggregation_method=primary.aggregation_method,
+                primary_metric_key=primary.metric_key,
+                metrics=metrics,
+            ),
+            unsupported_metrics=ThematicRasterAnalysisService._unsupported_metrics(product),
+            limitation_message=product.limitation_message,
+            generated_at=datetime.now(UTC).isoformat(),
+        )
+        return response.model_dump(mode="json")
+
+    @staticmethod
+    def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
+        dataset = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        product = ThematicRasterAnalysisService._product(dataset)
+        try:
+            import numpy as np
+            import rasterio
+            from PIL import Image
+            from rasterio.enums import Resampling
+        except ImportError as exc:
+            raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio, numpy and Pillow are required for thematic raster rendering", status_code=503) from exc
+        palettes = {
+            "space_occupation": np.asarray([[251, 231, 211], [190, 62, 51]], dtype="float64"),
+            "open_space": np.asarray([[221, 238, 219], [38, 122, 70]], dtype="float64"),
+            "population": np.asarray([[238, 231, 246], [103, 58, 151]], dtype="float64"),
+            "accessibility": np.asarray([[233, 241, 244], [15, 118, 110]], dtype="float64"),
+            "services": np.asarray([[255, 244, 191], [182, 109, 22]], dtype="float64"),
+        }
+        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))
+                resampling = Resampling.nearest if product.metric_kind == "binary_area" else Resampling.bilinear
+                data = source.read(1, out_shape=(height, width), masked=True, resampling=resampling)
+                values = np.asarray(data.filled(np.nan), dtype="float64")
+                valid = np.isfinite(values) & ~np.ma.getmaskarray(data)
+                if source.nodata is not None:
+                    valid &= ~np.isclose(values, float(source.nodata))
+                if product.metric_kind == "binary_area":
+                    valid &= values >= 0.5
+                    normalized = np.where(valid, 1.0, 0.0)
+                else:
+                    source_metadata = dataset.source_metadata or {}
+                    lower = float(source_metadata.get("render_min_value", np.nanpercentile(values[valid], 2) if valid.any() else 0.0))
+                    upper = float(source_metadata.get("render_max_value", np.nanpercentile(values[valid], 98) if valid.any() else 1.0))
+                    if upper <= lower:
+                        upper = lower + 1.0
+                    normalized = np.clip((values - lower) / (upper - lower), 0.0, 1.0)
+                colors = palettes[product.theme]
+                rgba = np.zeros((height, width, 4), dtype="uint8")
+                for channel in range(3):
+                    rgba[:, :, channel] = (colors[0, channel] + normalized * (colors[1, channel] - colors[0, channel])).astype("uint8")
+                rgba[:, :, 3] = np.where(valid, 205, 0).astype("uint8")
+                output = io.BytesIO()
+                Image.fromarray(rgba).save(output, format="PNG", optimize=True)
+                return output.getvalue()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="THEMATIC_RASTER_PREVIEW_FAILED",
+                message="The persisted thematic raster could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc

backend/app/services/vector_feature_service.py

@@ -29,6 +29,7 @@ FULL_AREA_CLIPPED_OPERATOR_TOOLS = {
     "provision_regional_bwk_natura2000.py",
     "provision_agricultural_parcel_history.py",
     "provision_buildings_addresses_register.py",
+    "provision_mol_soil_map.py",
 }
 
 SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS = {
@@ -101,6 +102,16 @@ SEMANTIC_SELECTION_METRICS: dict[str, tuple[dict[str, Any], ...]] = {
     ),
     "nature_value": (),
     "agriculture": (),
+    "soil": (
+        {
+            "metric_key": "soil_mapped_area",
+            "method": "intersection_area",
+            "label": "Bodemkaartoppervlakte",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "warning": "Historische bodemkartering op schaal 1:20.000; actuele lokale bodem- en drainagetoestand kan afwijken.",
+        },
+    ),
 }
 
 SEMANTIC_COUNT_LABELS = {
@@ -112,6 +123,7 @@ SEMANTIC_COUNT_LABELS = {
     "parcels": "Percelen",
     "nature_value": "BWK-kaartvlakken",
     "agriculture": "Landbouwgebruikspercelen",
+    "soil": "Bodemkaartvlakken",
 }
 
 # Sprint 205 initially normalized two official comma-separated ALZ group labels
@@ -158,6 +170,8 @@ class VectorFeatureService:
             "landbouw": "agriculture",
             "landbouwgebruik": "agriculture",
             "building_registry": "buildings",
+            "soil_map": "soil",
+            "bodem": "soil",
         }
         for candidate in candidates:
             if not isinstance(candidate, str) or not candidate.strip():

backend/tests/test_sprint213_thematic_rasters.py

@@ -0,0 +1,339 @@
+from __future__ import annotations
+
+from pathlib import Path
+from types import SimpleNamespace
+from uuid import uuid4
+
+import numpy as np
+import pytest
+import rasterio
+from fastapi.testclient import TestClient
+from pyproj import Transformer
+from rasterio.io import MemoryFile
+from rasterio.transform import from_origin
+
+from app.core.config import Settings
+from app.core.errors import AppError
+from app.db.session import get_db
+from app.main import app
+from app.models import Dataset, Job, Project
+from app.schemas.thematic_raster import ThematicRasterAcquireRequest, ThematicRasterSelectionRequest
+from app.schemas.assistant import AssistantQueryRequest
+from app.services.geo_assistant_service import GeoAssistantService
+from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService
+from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService
+from app.services.dataset_service import DatasetService
+
+
+ROOT = Path(__file__).resolve().parents[2]
+
+
+class FakeQuery:
+    def __init__(self, result=None):
+        self.result = result
+
+    def filter(self, *_args):
+        return self
+
+    def order_by(self, *_args):
+        return self
+
+    def first(self):
+        return self.result
+
+    def all(self):
+        return self.result if isinstance(self.result, list) else []
+
+
+class FakeSession:
+    def __init__(self, rows=None, query_result=None):
+        self.rows = rows or {}
+        self.query_result = query_result
+        self.added = []
+
+    def get(self, model, row_id):
+        row = self.rows.get((model, row_id))
+        if row is not None:
+            return row
+        return next((item for item in self.added if isinstance(item, model) and item.id == row_id), None)
+
+    def add(self, row):
+        self.added.append(row)
+
+    def commit(self):
+        return None
+
+    def rollback(self):
+        return None
+
+    def refresh(self, row):
+        return row
+
+    def query(self, _model):
+        return FakeQuery(self.query_result)
+
+
+class FakeResponse:
+    def __init__(self, content: bytes):
+        self.content = content
+        self.headers = {"Content-Type": "image/tiff", "Content-Length": str(len(content))}
+
+    def __enter__(self):
+        return self
+
+    def __exit__(self, *_args):
+        return None
+
+    def read(self, limit: int):
+        return self.content[:limit]
+
+
+def payload(product_key: str = "space_occupation_2025", *, side_m: float = 1000.0) -> ThematicRasterAcquireRequest:
+    transformer = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True)
+    min_x, min_y = transformer.transform(200_000, 210_000)
+    max_x, max_y = transformer.transform(200_000 + side_m, 210_000 + side_m)
+    return ThematicRasterAcquireRequest(
+        bbox={"min_x": min_x, "min_y": min_y, "max_x": max_x, "max_y": max_y, "crs": "EPSG:4326"},
+        product_key=product_key,
+        force_refresh=True,
+    )
+
+
+def raster_bytes(values: np.ndarray, resolution: float, *, nodata: float = -9999.0) -> bytes:
+    with MemoryFile() as memory:
+        with memory.open(
+            driver="GTiff",
+            width=values.shape[1],
+            height=values.shape[0],
+            count=1,
+            dtype=str(values.dtype),
+            crs="EPSG:31370",
+            transform=from_origin(200_000, 210_000 + values.shape[0] * resolution, resolution, resolution),
+            nodata=nodata,
+        ) as output:
+            output.write(values, 1)
+        return memory.read()
+
+
+def test_registry_contains_five_governed_non_water_policy_products() -> None:
+    products = ThematicRasterAcquisitionService.list_products()
+
+    assert [item["key"] for item in products] == [
+        "space_occupation_2025",
+        "open_space_2022",
+        "population_density_2019",
+        "node_value_2022",
+        "service_level_2022",
+    ]
+    assert {item["theme"] for item in products} == {"space_occupation", "open_space", "population", "accessibility", "services"}
+    assert {item["native_resolution_m"] for item in products} == {10.0, 100.0}
+    assert all(item["coverage_id"].startswith(("lu:", "ni:")) for item in products)
+    assert all(item["source_crs"] == "EPSG:31370" for item in products)
+    assert all(item["attribution"] and item["license_note"] and item["limitation_message"] for item in products)
+
+
+def test_request_is_bounded_allowlisted_and_uses_native_wcs_resolution() -> None:
+    settings = Settings(_env_file=None)
+    prepared = ThematicRasterAcquisitionService._prepared_request(payload("population_density_2019"), settings)
+    url = ThematicRasterAcquisitionService._wcs_request_url(settings, prepared["product"], tuple(prepared["bbox_epsg31370"]))
+
+    assert "VERSION=1.0.0" in url
+    assert "COVERAGE=ni%3Ani_inw_ha_vlaa_2019" in url
+    assert "RESX=100" in url and "RESY=100" in url
+    assert prepared["width"] * prepared["height"] <= settings.thematic_raster_max_pixels
+
+    with pytest.raises(AppError) as exc_info:
+        ThematicRasterAcquisitionService._prepared_request(payload("arbitrary_remote_layer"), settings)
+    assert exc_info.value.code == "THEMATIC_RASTER_PRODUCT_NOT_SUPPORTED"
+
+
+def test_binary_and_normalized_products_fail_closed_on_invalid_values() -> None:
+    binary = ThematicRasterAcquisitionService._product("space_occupation_2025")
+    score = ThematicRasterAcquisitionService._product("service_level_2022")
+
+    with pytest.raises(AppError, match="Binary"):
+        ThematicRasterAcquisitionService._validate_values(np.asarray([0.0, 2.0]), binary)
+    with pytest.raises(AppError, match="0-1"):
+        ThematicRasterAcquisitionService._validate_values(np.asarray([0.2, 1.2]), score)
+
+
+def test_acquisition_clips_validates_and_delegates_persistence(monkeypatch) -> None:
+    project_id, output_dataset_id = uuid4(), uuid4()
+    db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")})
+    content = raster_bytes(np.ones((100, 100), dtype="float32"), 10.0)
+    captured: dict = {}
+
+    def fake_import(_db, **kwargs):
+        captured.update(kwargs)
+        return SimpleNamespace(id=output_dataset_id)
+
+    monkeypatch.setattr(DatasetService, "import_raster_bytes", fake_import)
+    result = ThematicRasterAcquisitionService.acquire(
+        db,
+        project_id,
+        payload(side_m=1000.0),
+        settings=Settings(_env_file=None),
+        opener=lambda *_args, **_kwargs: FakeResponse(content),
+    )
+
+    assert result["output_dataset_id"] == str(output_dataset_id)
+    assert captured["source_name"] == ThematicRasterAcquisitionService.PROVIDER
+    assert captured["source_metadata"]["product_key"] == "space_occupation_2025"
+    assert captured["source_metadata"]["metric_kind"] == "binary_area"
+    assert captured["source_metadata"]["valid_pixel_count"] > 9_800
+    assert captured["provenance_metadata"]["acquisition"] == "explicit_bounded_tiled_wcs_coverage"
+    assert len(captured["provenance_metadata"]["normalized_sha256"]) == 64
+
+
+def test_binary_area_analysis_returns_hectares_and_share(tmp_path) -> None:
+    project_id, dataset_id = uuid4(), uuid4()
+    values = np.zeros((10, 10), dtype="float32")
+    values[:, :5] = 1.0
+    path = tmp_path / "space.tif"
+    path.write_bytes(raster_bytes(values, 10.0))
+    dataset = Dataset(
+        id=dataset_id,
+        project_id=project_id,
+        name="space.tif",
+        dataset_type="raster",
+        source="official",
+        source_name=ThematicRasterAcquisitionService.PROVIDER,
+        source_metadata={"product_key": "space_occupation_2025", "coverage_id": "lu:lu_ruibes_vlaa_2025_v3"},
+        status="ready",
+        storage_path=str(path),
+    )
+    db = FakeSession({(Dataset, dataset_id): dataset})
+    result = ThematicRasterAnalysisService.analyze(
+        db,
+        project_id,
+        dataset_id,
+        ThematicRasterSelectionRequest(bbox=payload(side_m=100.0).bbox),
+    )
+    metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]}
+
+    assert result["valid_cell_count"] == 100
+    assert metrics["space_occupation_area_ha"]["metric_value"] == pytest.approx(0.5)
+    assert metrics["space_occupation_share_pct"]["metric_value"] == pytest.approx(50.0)
+    assert "object_count" in result["unsupported_metrics"]
+
+
+def test_population_analysis_sums_one_hectare_density_cells_without_claiming_current_counts(tmp_path) -> None:
+    project_id, dataset_id = uuid4(), uuid4()
+    values = np.asarray([[10.0, 20.0], [30.0, 40.0]], dtype="float32")
+    path = tmp_path / "population.tif"
+    path.write_bytes(raster_bytes(values, 100.0))
+    dataset = Dataset(
+        id=dataset_id,
+        project_id=project_id,
+        name="population.tif",
+        dataset_type="raster",
+        source="official",
+        source_name=ThematicRasterAcquisitionService.PROVIDER,
+        source_metadata={"product_key": "population_density_2019", "coverage_id": "ni:ni_inw_ha_vlaa_2019"},
+        status="ready",
+        storage_path=str(path),
+    )
+    db = FakeSession({(Dataset, dataset_id): dataset})
+    result = ThematicRasterAnalysisService.analyze(
+        db,
+        project_id,
+        dataset_id,
+        ThematicRasterSelectionRequest(bbox=payload("population_density_2019", side_m=200.0).bbox),
+    )
+    metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]}
+
+    assert metrics["estimated_inhabitants"]["metric_value"] == pytest.approx(100.0)
+    assert metrics["population_density_mean_per_ha"]["metric_value"] == pytest.approx(25.0)
+    assert metrics["estimated_inhabitants"]["is_estimate"] is True
+    assert "current_population" in result["unsupported_metrics"]
+
+
+def test_assistant_context_receives_persisted_thematic_metrics(tmp_path) -> None:
+    project_id, dataset_id = uuid4(), uuid4()
+    values = np.asarray([[10.0, 20.0], [30.0, 40.0]], dtype="float32")
+    path = tmp_path / "assistant-population.tif"
+    path.write_bytes(raster_bytes(values, 100.0))
+    project = Project(id=project_id, name="Mol", region="Mol")
+    dataset = Dataset(
+        id=dataset_id,
+        project_id=project_id,
+        name="population.tif",
+        dataset_type="raster",
+        source="official",
+        source_name=ThematicRasterAcquisitionService.PROVIDER,
+        source_metadata={"product_key": "population_density_2019", "coverage_id": "ni:ni_inw_ha_vlaa_2019"},
+        status="ready",
+        storage_path=str(path),
+    )
+    db = FakeSession({(Project, project_id): project, (Dataset, dataset_id): dataset}, query_result=[dataset])
+    context, metrics, _series, dataset_ids, _warnings, _scope = GeoAssistantService(Settings(_env_file=None))._build_context(
+        db,
+        project_id=project_id,
+        payload=AssistantQueryRequest(question="Hoeveel inwoners?", bbox=payload("population_density_2019", side_m=200.0).bbox),
+    )
+
+    assert any(metric.theme == "population" and metric.label.startswith("Geraamd aantal") for metric in metrics)
+    assert dataset_id in dataset_ids
+    assert context["rules"]["thematic_policy_rasters_available"] is True
+
+
+def test_index_renderer_returns_browser_png(tmp_path) -> None:
+    project_id, dataset_id = uuid4(), uuid4()
+    values = np.linspace(0.1, 4.0, 100, dtype="float32").reshape((10, 10))
+    path = tmp_path / "node.tif"
+    path.write_bytes(raster_bytes(values, 100.0))
+    dataset = Dataset(
+        id=dataset_id,
+        project_id=project_id,
+        name="node.tif",
+        dataset_type="raster",
+        source="official",
+        source_name=ThematicRasterAcquisitionService.PROVIDER,
+        source_metadata={
+            "product_key": "node_value_2022",
+            "coverage_id": "lu:lu_knptw_ha_2022_v3",
+            "render_min_value": 0.1,
+            "render_max_value": 4.0,
+        },
+        status="ready",
+        storage_path=str(path),
+    )
+    db = FakeSession({(Dataset, dataset_id): dataset})
+
+    assert ThematicRasterAnalysisService.render_png(db, project_id, dataset_id).startswith(b"\x89PNG\r\n\x1a\n")
+
+
+def test_api_uses_canonical_envelopes(monkeypatch) -> None:
+    project_id, dataset_id = uuid4(), uuid4()
+    db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")})
+    monkeypatch.setattr(
+        ThematicRasterAcquisitionService,
+        "acquire",
+        lambda *_args, **_kwargs: {"output_dataset_id": str(dataset_id), "provider": ThematicRasterAcquisitionService.PROVIDER},
+    )
+    monkeypatch.setattr(
+        ThematicRasterAnalysisService,
+        "analyze",
+        lambda *_args, **_kwargs: {"dataset_id": str(dataset_id), "theme": "population", "summary": {"metric_value": 10.0}},
+    )
+    app.dependency_overrides[get_db] = lambda: db
+    try:
+        client = TestClient(app)
+        products = client.get(f"/api/v1/projects/{project_id}/datasets/thematic-raster/products")
+        acquisition = client.post(
+            f"/api/v1/projects/{project_id}/datasets/thematic-raster/acquire",
+            json=payload().model_dump(mode="json"),
+        )
+        selection = client.post(
+            f"/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select",
+            json={"bbox": payload().bbox.model_dump()},
+        )
+    finally:
+        app.dependency_overrides.clear()
+
+    assert products.status_code == 200 and set(products.json()) == {"data"}
+    assert products.json()["data"]["total"] == 5
+    assert acquisition.status_code == 200 and set(acquisition.json()) == {"data"}
+    assert acquisition.json()["data"]["job_type"] == "raster.thematic.acquire"
+    assert selection.status_code == 200 and selection.json()["data"]["theme"] == "population"
+    assert any(isinstance(item, Job) for item in db.added)

backend/tests/test_sprint214_dov_soil_map.py

@@ -0,0 +1,183 @@
+from __future__ import annotations
+
+import importlib.util
+from pathlib import Path
+from uuid import uuid4
+
+import pytest
+from shapely.geometry import box, mapping, shape
+from shapely.ops import transform as transform_geometry
+
+from app.models import Dataset
+from app.services.vector_feature_service import VectorFeatureService
+
+
+ROOT = Path(__file__).resolve().parents[2]
+
+
+def load_operator():
+    path = ROOT / "scripts" / "provision_mol_soil_map.py"
+    spec = importlib.util.spec_from_file_location("dov_soil_map_operator", path)
+    assert spec and spec.loader
+    module = importlib.util.module_from_spec(spec)
+    spec.loader.exec_module(module)
+    return module
+
+
+class FakeResponse:
+    def __init__(self, payload: dict, url: str):
+        self._payload = payload
+        self.url = url
+        self.content = b'{"type":"FeatureCollection"}'
+
+    def raise_for_status(self) -> None:
+        return None
+
+    def json(self) -> dict:
+        return self._payload
+
+
+class FakeSession:
+    def __init__(self, pages: list[dict]):
+        self.pages = pages
+        self.calls: list[dict] = []
+
+    def get(self, _url: str, *, params: dict, timeout: int):
+        self.calls.append({"params": dict(params), "timeout": timeout})
+        return FakeResponse(self.pages[len(self.calls) - 1], f"https://example.test/page/{len(self.calls)}")
+
+
+def soil_feature(module, feature_id: str = "bodemtypes.1") -> dict:
+    return {
+        "type": "Feature",
+        "id": feature_id,
+        "geometry": mapping(box(5.0, 51.0, 5.02, 51.02)),
+        "properties": {
+            "gid": 1,
+            "id_kaartvlak": 10,
+            "Bodemtype": "Zeg",
+            "Unibodemtype": "Zeg",
+            "Bodemserie": "Zeg",
+            "Beknopte_omschrijving_bodemserie": "Natte zandbodem",
+            "Gegeneraliseerde_legende": "Nat zand",
+            "Textuurklasse_code": "Z",
+            "Textuurklasse": "zand",
+            "Drainageklasse_code": "e",
+            "Drainageklasse": "nat",
+            "Profielontwikkelingsgroep_code": "g",
+            "Profielontwikkelingsgroep": "humus B horizont",
+            "Eenduidige_legende_titel": "bodemserie Zeg",
+        },
+    }
+
+
+def test_wfs_pagination_is_bounded_complete_and_deterministic() -> None:
+    module = load_operator()
+    feature = soil_feature(module)
+    pages = [
+        {
+            "type": "FeatureCollection",
+            "numberMatched": 3,
+            "numberReturned": 2,
+            "features": [feature, {**feature, "id": "bodemtypes.2"}],
+        },
+        {
+            "type": "FeatureCollection",
+            "numberMatched": 3,
+            "numberReturned": 1,
+            "features": [{**feature, "id": "bodemtypes.3"}],
+        },
+    ]
+    session = FakeSession(pages)
+
+    result = list(
+        module.iter_wfs_pages(
+            session,
+            (196000.0, 205000.0, 211000.0, 224000.0),
+            page_limit=2,
+            timeout=30,
+        )
+    )
+
+    assert len(result) == 2
+    assert [call["params"]["startIndex"] for call in session.calls] == ["0", "2"]
+    assert all(call["params"]["typeNames"] == "bodemkaart:bodemtypes" for call in session.calls)
+    assert all(call["params"]["bbox"].endswith("EPSG:31370") for call in session.calls)
+    assert all(call["params"]["sortBy"] == "gid" for call in session.calls)
+
+
+def test_soil_feature_is_exactly_clipped_and_keeps_governed_properties() -> None:
+    module = load_operator()
+    boundary_wgs84 = box(5.005, 51.005, 5.015, 51.015)
+    boundary_lambert72 = transform_geometry(module.TO_LAMBERT72.transform, boundary_wgs84)
+
+    normalized, was_clipped = module.normalize_feature(soil_feature(module), boundary_lambert72)
+
+    assert normalized is not None and was_clipped is True
+    persisted_geometry = shape(normalized["geometry"])
+    assert persisted_geometry.within(boundary_wgs84.buffer(1e-7))
+    properties = normalized["properties"]
+    assert properties["source_name"] == "dov_soil_map"
+    assert properties["soil_texture_class"] == "zand"
+    assert properties["soil_drainage_class"] == "nat"
+    assert properties["survey_period"] == "1949-1971"
+    assert properties["clipped_area_ha"] > 0
+    assert "may differ today" in properties["historical_drainage_limitation"]
+
+
+def test_soil_map_uses_existing_semantic_selection_architecture() -> None:
+    dataset = Dataset(
+        id=uuid4(),
+        project_id=uuid4(),
+        name="dov_soil_map_mol.geojson",
+        dataset_type="vector",
+        source="operator_official_import",
+        source_name="dov_soil_map",
+        reference_layer_name="soil",
+        source_metadata={
+            "theme": "soil",
+            "selection_aggregation": {
+                "method": "intersection_area",
+                "label": "Bodemkaartoppervlakte",
+                "unit": "ha",
+            },
+        },
+        status="ready",
+    )
+
+    assert VectorFeatureService._dataset_theme(dataset) == "soil"
+    assert VectorFeatureService.supports_selection_summary(dataset) is True
+    assert VectorFeatureService.can_use_full_area_fast_path(dataset, None) is False
+
+
+def test_soil_operator_contract_has_no_direct_persistence_and_is_packaged() -> None:
+    operator = (ROOT / "scripts" / "provision_mol_soil_map.py").read_text(encoding="utf-8")
+    dockerfile = (ROOT / "deploy" / "unraid" / "Dockerfile.all-in-one").read_text(encoding="utf-8")
+    readiness = (ROOT / "scripts" / "run_readiness_check.sh").read_text(encoding="utf-8")
+    map_workspace = (ROOT / "frontend" / "src" / "components" / "map" / "MapWorkspace.tsx").read_text(encoding="utf-8")
+
+    assert "/datasets/upload" in operator
+    assert "vector_features" in operator
+    assert "does not write directly" in " ".join(operator.split())
+    assert "SessionLocal" not in operator and "INSERT INTO" not in operator
+    assert "COPY scripts/provision_mol_soil_map.py" in dockerfile
+    assert "py_compile scripts/provision_mol_soil_map.py" in readiness
+    assert "id: 'soil'" in map_workspace
+    assert "dataset.source_name === 'dov_soil_map'" in map_workspace
+
+
+def test_incomplete_wfs_pagination_fails_closed() -> None:
+    module = load_operator()
+    session = FakeSession(
+        [
+            {
+                "type": "FeatureCollection",
+                "numberMatched": 2,
+                "numberReturned": 0,
+                "features": [],
+            }
+        ]
+    )
+
+    with pytest.raises(RuntimeError, match="returned 0 of 2"):
+        list(module.iter_wfs_pages(session, (0.0, 0.0, 1.0, 1.0), page_limit=100, timeout=30))

deploy/unraid/Dockerfile.all-in-one

@@ -77,6 +77,8 @@ COPY scripts/provision_mol_context_layers.py /app/scripts/provision_mol_context_
 COPY scripts/provision_mol_dhmv.py /app/scripts/provision_mol_dhmv.py
 COPY scripts/provision_mol_flood_hazards.py /app/scripts/provision_mol_flood_hazards.py
 COPY scripts/provision_regional_flood_hazards.py /app/scripts/provision_regional_flood_hazards.py
+COPY scripts/provision_thematic_rasters.py /app/scripts/provision_thematic_rasters.py
+COPY scripts/provision_mol_soil_map.py /app/scripts/provision_mol_soil_map.py
 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_regional_historical_landuse.py /app/scripts/provision_regional_historical_landuse.py

deploy/unraid/geointel-unraid-template.xml

@@ -43,6 +43,10 @@
   <Config Name="VMM Flood Hazard WCS URL" Target="FLOOD_HAZARD_WCS_URL" Default="https://geoservice.waterinfo.be/OGRK/wcs" Mode="" Description="Official VMM OGRK WCS endpoint for governed flood-depth scenarios." Type="Variable" Display="advanced" Required="true" Mask="false">https://geoservice.waterinfo.be/OGRK/wcs</Config>
   <Config Name="Flood Hazard Analysis Resolution (m)" Target="FLOOD_HAZARD_RESOLUTION_M" Default="5.0" Mode="" Description="Stored analysis grid resolution; official source values are converted from centimetres to metres." Type="Variable" Display="advanced" Required="true" Mask="false">5.0</Config>
   <Config Name="Flood Hazard Maximum Cells" Target="FLOOD_HAZARD_MAX_PIXELS" Default="12000000" Mode="" Description="Maximum raster cells per flood-hazard acquisition or selection analysis." Type="Variable" Display="advanced" Required="true" Mask="false">12000000</Config>
+  <Config Name="Official Thematic Raster Acquisition" Target="THEMATIC_RASTER_ENABLED" Default="true" Mode="" Description="Allow bounded official Departement Omgeving rasters for space, population, accessibility and services." Type="Variable" Display="advanced" Required="true" Mask="false">true</Config>
+  <Config Name="Thematic Raster WCS URL" Target="THEMATIC_RASTER_WCS_URL" Default="https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs" Mode="" Description="Official public MercatorNet WCS endpoint. Product identifiers remain server allowlisted." Type="Variable" Display="advanced" Required="true" Mask="false">https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs</Config>
+  <Config Name="Thematic Raster Maximum Side (m)" Target="THEMATIC_RASTER_MAX_SIDE_M" Default="20000" Mode="" Description="Maximum side length for one persisted thematic raster scope." Type="Variable" Display="advanced" Required="true" Mask="false">20000</Config>
+  <Config Name="Thematic Raster Maximum Cells" Target="THEMATIC_RASTER_MAX_PIXELS" Default="12000000" Mode="" Description="Maximum raster cells per acquisition or selection analysis." Type="Variable" Display="advanced" Required="true" Mask="false">12000000</Config>
   <Config Name="Local Ollama Assistant" Target="OLLAMA_ENABLED" Default="true" Mode="" Description="Enable the source-grounded GeoIntel assistant backed by Ollama on the Unraid host." Type="Variable" Display="always" Required="true" Mask="false">true</Config>
   <Config Name="Ollama Base URL" Target="OLLAMA_BASE_URL" Default="http://host.docker.internal:11434" Mode="" Description="Ollama API reachable from the container. The deployment maps host.docker.internal to the Unraid host gateway." Type="Variable" Display="always" Required="true" Mask="false">http://host.docker.internal:11434</Config>
   <Config Name="Default Ollama Model" Target="OLLAMA_DEFAULT_MODEL" Default="qwen3.5:9b" Mode="" Description="Preferred locally installed Ollama model. Users can select another installed model in GeoIntel." Type="Variable" Display="always" Required="true" Mask="false">qwen3.5:9b</Config>

deploy/unraid/run-dockerman-container.sh

@@ -43,6 +43,13 @@ FLOOD_HAZARD_MAX_SIDE_M="${FLOOD_HAZARD_MAX_SIDE_M:-20000}"
 FLOOD_HAZARD_MAX_PIXELS="${FLOOD_HAZARD_MAX_PIXELS:-12000000}"
 FLOOD_HAZARD_TIMEOUT_SECONDS="${FLOOD_HAZARD_TIMEOUT_SECONDS:-300}"
 FLOOD_HAZARD_MAX_RESPONSE_MB="${FLOOD_HAZARD_MAX_RESPONSE_MB:-160}"
+THEMATIC_RASTER_ENABLED="${THEMATIC_RASTER_ENABLED:-true}"
+THEMATIC_RASTER_WCS_URL="${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs}"
+THEMATIC_RASTER_MIN_SIDE_M="${THEMATIC_RASTER_MIN_SIDE_M:-100}"
+THEMATIC_RASTER_MAX_SIDE_M="${THEMATIC_RASTER_MAX_SIDE_M:-20000}"
+THEMATIC_RASTER_MAX_PIXELS="${THEMATIC_RASTER_MAX_PIXELS:-12000000}"
+THEMATIC_RASTER_TIMEOUT_SECONDS="${THEMATIC_RASTER_TIMEOUT_SECONDS:-300}"
+THEMATIC_RASTER_MAX_RESPONSE_MB="${THEMATIC_RASTER_MAX_RESPONSE_MB:-160}"
 YOLO_ENABLED="${YOLO_ENABLED:-false}"
 YOLO_MODELS_DIR="${YOLO_MODELS_DIR:-/app/models}"
 YOLO_MODEL_PATH="${YOLO_MODEL_PATH:-}"
@@ -135,6 +142,13 @@ docker run -d \
   -e FLOOD_HAZARD_MAX_PIXELS="$FLOOD_HAZARD_MAX_PIXELS" \
   -e FLOOD_HAZARD_TIMEOUT_SECONDS="$FLOOD_HAZARD_TIMEOUT_SECONDS" \
   -e FLOOD_HAZARD_MAX_RESPONSE_MB="$FLOOD_HAZARD_MAX_RESPONSE_MB" \
+  -e THEMATIC_RASTER_ENABLED="$THEMATIC_RASTER_ENABLED" \
+  -e THEMATIC_RASTER_WCS_URL="$THEMATIC_RASTER_WCS_URL" \
+  -e THEMATIC_RASTER_MIN_SIDE_M="$THEMATIC_RASTER_MIN_SIDE_M" \
+  -e THEMATIC_RASTER_MAX_SIDE_M="$THEMATIC_RASTER_MAX_SIDE_M" \
+  -e THEMATIC_RASTER_MAX_PIXELS="$THEMATIC_RASTER_MAX_PIXELS" \
+  -e THEMATIC_RASTER_TIMEOUT_SECONDS="$THEMATIC_RASTER_TIMEOUT_SECONDS" \
+  -e THEMATIC_RASTER_MAX_RESPONSE_MB="$THEMATIC_RASTER_MAX_RESPONSE_MB" \
   -e YOLO_ENABLED="$YOLO_ENABLED" \
   -e YOLO_MODELS_DIR="$YOLO_MODELS_DIR" \
   -e YOLO_MODEL_PATH="$YOLO_MODEL_PATH" \

docker-compose.unraid.yml

@@ -41,6 +41,13 @@ services:
       FLOOD_HAZARD_MAX_PIXELS: ${FLOOD_HAZARD_MAX_PIXELS:-12000000}
       FLOOD_HAZARD_TIMEOUT_SECONDS: ${FLOOD_HAZARD_TIMEOUT_SECONDS:-300}
       FLOOD_HAZARD_MAX_RESPONSE_MB: ${FLOOD_HAZARD_MAX_RESPONSE_MB:-160}
+      THEMATIC_RASTER_ENABLED: ${THEMATIC_RASTER_ENABLED:-true}
+      THEMATIC_RASTER_WCS_URL: ${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs}
+      THEMATIC_RASTER_MIN_SIDE_M: ${THEMATIC_RASTER_MIN_SIDE_M:-100}
+      THEMATIC_RASTER_MAX_SIDE_M: ${THEMATIC_RASTER_MAX_SIDE_M:-20000}
+      THEMATIC_RASTER_MAX_PIXELS: ${THEMATIC_RASTER_MAX_PIXELS:-12000000}
+      THEMATIC_RASTER_TIMEOUT_SECONDS: ${THEMATIC_RASTER_TIMEOUT_SECONDS:-300}
+      THEMATIC_RASTER_MAX_RESPONSE_MB: ${THEMATIC_RASTER_MAX_RESPONSE_MB:-160}
       YOLO_ENABLED: ${YOLO_ENABLED:-false}
       YOLO_MODELS_DIR: ${YOLO_MODELS_DIR:-/app/models}
       YOLO_MODEL_PATH: ${YOLO_MODEL_PATH:-}

docker-compose.yml

@@ -46,6 +46,13 @@ services:
       FLOOD_HAZARD_MAX_PIXELS: ${FLOOD_HAZARD_MAX_PIXELS:-12000000}
       FLOOD_HAZARD_TIMEOUT_SECONDS: ${FLOOD_HAZARD_TIMEOUT_SECONDS:-300}
       FLOOD_HAZARD_MAX_RESPONSE_MB: ${FLOOD_HAZARD_MAX_RESPONSE_MB:-160}
+      THEMATIC_RASTER_ENABLED: ${THEMATIC_RASTER_ENABLED:-true}
+      THEMATIC_RASTER_WCS_URL: ${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs}
+      THEMATIC_RASTER_MIN_SIDE_M: ${THEMATIC_RASTER_MIN_SIDE_M:-100}
+      THEMATIC_RASTER_MAX_SIDE_M: ${THEMATIC_RASTER_MAX_SIDE_M:-20000}
+      THEMATIC_RASTER_MAX_PIXELS: ${THEMATIC_RASTER_MAX_PIXELS:-12000000}
+      THEMATIC_RASTER_TIMEOUT_SECONDS: ${THEMATIC_RASTER_TIMEOUT_SECONDS:-300}
+      THEMATIC_RASTER_MAX_RESPONSE_MB: ${THEMATIC_RASTER_MAX_RESPONSE_MB:-160}
       YOLO_ENABLED: ${YOLO_ENABLED:-false}
       YOLO_MODELS_DIR: ${YOLO_MODELS_DIR:-/app/models}
       YOLO_MODEL_PATH: ${YOLO_MODEL_PATH:-}

docs/API_CONTRACTS.md

@@ -322,6 +322,53 @@ Returns a constrained transparent PNG for a persisted governed VMM flood-depth
 Dataset. It never accepts an arbitrary path or coverage id and is used by the
 existing MapLibre image-overlay path.
 
+### GET `/api/v1/projects/{project_id}/datasets/thematic-raster/products`
+
+Returns the fixed MercatorNet registry for `space_occupation_2025`,
+`open_space_2022`, `population_density_2019`, `node_value_2022` and
+`service_level_2022`. Every item includes the governed WCS coverage id, native
+resolution, source unit, observation year, legend, attribution and limitation.
+
+### POST `/api/v1/projects/{project_id}/datasets/thematic-raster/acquire`
+
+Acquires one allowlisted official coverage behind the synchronous Job
+abstraction. The request accepts only an EPSG:4326 bbox, optional project Area,
+one registry product key and an explicit refresh flag:
+
+```json
+{
+  "bbox": {"min_x": 5.03, "min_y": 51.15, "max_x": 5.25, "max_y": 51.33, "crs": "EPSG:4326"},
+  "area_id": "optional-project-area-uuid",
+  "product_key": "population_density_2019",
+  "force_refresh": false
+}
+```
+
+The backend uses native 10 m or 100 m resolution, splits requests into bounded
+WCS 1.0 tiles, validates EPSG:31370 and documented source values, masks the
+exact Area and persists an ordinary raster Dataset and DatasetVersion. It does
+not accept arbitrary URLs, coverage ids, resolutions or expressions.
+
+### POST `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select`
+
+Returns source-correct metrics for a bbox and optional exact Area mask:
+
+- occupied/open hectares, share and valid raster area for binary products;
+- estimated inhabitants plus mean/P90 inhabitants per hectare for the 2019
+  population raster;
+- mean, P10, median and P90 source score for node value and service level.
+
+The response names estimate status, aggregation method, source unit,
+observation year, attribution, unsupported metrics and product limitation.
+Current register population, live public-transport availability and causal
+interpretations are not produced.
+
+### GET `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/image`
+
+Returns a constrained transparent PNG generated from the persisted governed
+raster. It accepts neither an arbitrary file path nor a provider URL and feeds
+the existing MapLibre image-overlay path.
+
 ### GET `/api/v1/projects/{project_id}/datasets`
 
 List datasets.

docs/CODEX_EXECUTION_LOG.md

@@ -1,3 +1,29 @@
+## Sprint 213-214 Cross-domain thematic rasters and DOV soil map (2026-07-16)
+
+Changed:
+- Added a governed five-product MercatorNet thematic-raster registry covering
+  space, people, accessibility and services with fixed coverage ids, native
+  resolutions, units, years, legends and limitations.
+- Added safe WCS tiling, exact Area masking, canonical raster persistence,
+  MapLibre rendering, selection metrics and persisted-metric Ollama context.
+- Added an explicit DOV soil-map operator for Mol with deterministic WFS 2.0
+  pagination, exact EPSG:31370 clipping, raw checksums and canonical vector
+  upload. No source writes directly to PostGIS.
+- Added end-user themes for space occupation, open space, population,
+  accessibility, services and soil. Current-only official states no longer
+  make Evolution appear available without at least two observations.
+
+Validated during implementation:
+- Focused thematic and soil suites passed with 14 tests.
+- Frontend TypeScript typecheck passed after map and API integration.
+- Live source contracts were checked against MercatorNet WCS and the DOV
+  production WFS; live runtime provisioning follows after deployment.
+
+Next:
+- Run the complete readiness gate, deploy Tower, provision all six Mol layers,
+  verify live PostGIS metrics and inspect the map at desktop, widescreen and
+  mobile widths before starting regional rollout.
+
 ## Sprint 195 Guided raster-to-detection workflow (2026-07-14)
 
 Changed:

docs/DATAVINDPLAATS_SOURCE_ROADMAP.md

@@ -28,8 +28,10 @@ CRS/unit validation, persisted provenance and selection-metric tests pass.
 - Historical context: official 1778, 1873 and 1969 mapped land-use classes and
   modern land-use editions where persisted.
 
-This is a strong base, but population/services, soil classes and modeled
+The cross-domain Wave 1 implementation now adds population/services, modeled
-accessibility are still genuine coverage gaps.
+accessibility, space occupation, open space and the DOV soil classes to Mol.
+Regional completeness beyond Mol remains a controlled operator rollout rather
+than an automatic startup fetch.
 
 ## Domain A - Space, Buildings And Economy
 
@@ -142,6 +144,13 @@ DatasetService/VectorFeatureService path. This immediately gives every drawn
 rectangle a balanced profile across space, soil, people, services and
 accessibility.
 
+Implementation status on 2026-07-16: complete for Mol. The five raster
+products use one fixed MercatorNet registry and the soil map uses the official
+`bodemkaart:bodemtypes` DOV WFS. All outputs are persisted through existing
+Dataset services, carry source year/period and limitations, and are selectable
+on the map. Regional raster provisioning is supported per persisted
+municipality; regional soil partitioning remains the next scale-out step.
+
 ### Wave 2 - Economic and mobility objects
 
 Add governed vector operators for business parks, Hoppinpoints and cycle

docs/DATA_SOURCES.md

@@ -26,6 +26,43 @@ aanvraag. GeoIntel verzint geen historische pixelopnamedatum. Bronnen:
 
 Dit document verzamelt concrete databronnen voor GeoIntel Kempen.
 
+## Cross-domain official area profile
+
+GeoIntel uses one governed MercatorNet WCS registry for five non-water policy
+rasters. The operator is explicit and bounded; no source fetch occurs during
+startup or directly from the browser.
+
+| Product | Coverage | Resolution | Selection output |
+| --- | --- | ---: | --- |
+| Ruimtebeslag 2025 | `lu:lu_ruibes_vlaa_2025_v3` | 10 m | occupied hectares and share |
+| Open ruimte 2022 | `lu:lu_openruimte_vlaa_2022_v3` | 10 m | open-space hectares and share |
+| Inwonersdichtheid 2019 | `ni:ni_inw_ha_vlaa_2019` | 100 m | estimated inhabitants and inhabitants/ha |
+| Knooppuntwaarde 2022 | `lu:lu_knptw_ha_2022_v3` | 100 m | source-score distribution |
+| Voorzieningenniveau 2022 | `lu:lu_totvznv_ha_2022_v3` | 100 m | normalized source-score distribution |
+
+All coverages are stored in EPSG:31370 with their native cells, exact Area
+mask, request/checksum provenance and reference year. Space occupation is not
+the same as buildings or paving. Open space is not automatically nature or
+publicly accessible land. Population is a 2019 raster estimate, not a current
+register count. Accessibility and service scores are not live travel times or
+object counts.
+
+## DOV digital soil map
+
+- Service: `https://www.dov.vlaanderen.be/geoserver/wfs`
+- Layer: `bodemkaart:bodemtypes`
+- Source CRS and metric clipping: EPSG:31370
+- Persisted geometry: EPSG:4326 through DatasetService/VectorFeatureService
+- Scale: 1:20,000
+- Survey evidence: field data collected between 1949 and 1971
+
+The Mol operator follows all WFS pages, keeps checksummed raw responses and
+persists soil type, series, generalized legend, texture, drainage, profile and
+substrate fields. Selection output reports intersected mapped hectares and
+fixed sand/anthropogenic context classes. The drainage field remains historical
+baseline evidence; GeoIntel does not claim it describes current parcel
+drainage or replace a site investigation.
+
 ## GRB — Basiskaart Vlaanderen
 
 - Naam: Basiskaart Vlaanderen / GRB

docs/STORAGE_ARCHITECTURE.md

@@ -205,6 +205,29 @@ The smoke never passes `--apply`. It fails if the cleanup summary is not a
 dry-run, if any export/file deletion is reported, or if the dry-run candidate
 fields are missing.
 
+## Governed cross-domain raster and soil evidence
+
+MercatorNet thematic products use the ordinary raster Dataset and immutable
+DatasetVersion paths under `STORAGE_ROOT`. `source_metadata` records product,
+coverage, native resolution/unit, render range, reference year and request
+bounds. `provenance_metadata` records every tiled WCS URL, transfer checksum,
+normalized checksum, exact Area id and validation result. No new storage table
+or provider-side path is introduced.
+
+The DOV Mol soil operator retains evidence under:
+
+```text
+storage/operator-evidence/dov-soil-map/mol/
+  raw/dov_soil_map_page_*.json
+  dov_soil_map_mol.geojson
+  dov_soil_map_mol.manifest.json
+```
+
+Only the normalized, exactly clipped GeoJSON enters DatasetService and
+`vector_features`. The manifest binds the persisted artifact, Mol boundary,
+source page URLs/checksums, feature completeness, class-area summaries and
+historical limitations. Raw source responses remain operator evidence.
+
 ## Model storage
 
 Model artifacts live under:

docs/TODO.md

@@ -30,6 +30,10 @@
 - [x] Add annual agricultural-use parcels through an explicit provider/operator contract.
 - [x] 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.
+- [x] Add governed official space-occupation, open-space, population-density, accessibility and service-level rasters with semantic map metrics.
+- [x] Add the DOV digital soil map for Mol with exact clipping, source attributes and historical survey limitations.
+- [ ] Execute the governed thematic-raster operator for every persisted Kempen municipality after the Mol live gate passes.
+- [ ] Generalize the DOV soil-map operator to all 28 approved Kempen municipality partitions with one regional snapshot manifest.
 
 ## Governed source expansion backlog
 

frontend/README.md

@@ -555,6 +555,14 @@ provider technology: space/buildings, nature/agriculture, soil/relief,
 mobility/accessibility, population/services and climate/living environment.
 The central definitions live in `src/lib/sourcePortfolio.ts`.
 
+The current Map explorer also recognizes the governed `Ruimtebeslag`, `Open
+ruimte`, `Bevolking`, `Bereikbaarheid`, `Voorzieningen` and `Bodem` Datasets.
+The first five render persisted raster PNGs with product-specific legends and
+return cell-based semantic metrics; `Bodem` renders the persisted DOV polygons
+and exposes soil attributes on feature selection. A single official snapshot
+never activates Evolution by itself. Evolution is enabled only when at least
+two comparable observations exist.
+
 The domain cards count only matching `ready` Datasets as operational. Audited
 official sources that have not passed acquisition, persistence and metric
 validation remain inside the collapsed follow-up list with an explicit

frontend/src/App.tsx

@@ -191,7 +191,10 @@ function App(): JSX.Element {
       const floodHazards = datasets.filter(
         (dataset) => dataset.dataset_type === 'raster' && dataset.source_name === 'vmm_flood_hazard' && dataset.status === 'ready',
       )
-      return [...vectors, ...terrain, ...floodHazards]
+      const thematicRasters = datasets.filter(
+        (dataset) => dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster' && dataset.status === 'ready',
+      )
+      return [...vectors, ...terrain, ...floodHazards, ...thematicRasters]
     },
     [datasets],
   )

frontend/src/components/map/MapWorkspace.tsx

@@ -8,6 +8,7 @@ import { getDatasetDisplayName, getDatasetSourceDisplayName } from '../../lib/da
 import { TemporalTrendChart } from './TemporalTrendChart'
 import { terrainImageUrl } from '../../lib/terrainImage'
 import { floodHazardImageUrl } from '../../lib/floodHazardImage'
+import { thematicRasterImageUrl } from '../../lib/thematicRaster'
 
 const DEFAULT_SELECTED_FEATURE_FILENAME = 'selected-feature.geojson'
 const DEFAULT_AREA_SELECTION_FILENAME = 'area-selection.geojson'
@@ -15,7 +16,7 @@ const EMPTY_TEMPORAL_SERIES: DatasetCreateResponse[] = []
 const MOL_PROJECT_NAME = 'Mol Municipality Workbench'
 const KEMPEN_PROJECT_NAME = 'Kempen Regional Workbench'
 
-type DataThemeId = 'buildings' | 'population' | 'forest' | 'nature_value' | 'agriculture' | 'water' | 'flood_hazard' | 'elevation' | 'roads' | 'parcels'
+type DataThemeId = 'buildings' | 'space_occupation' | 'open_space' | 'population' | 'forest' | 'nature_value' | 'agriculture' | 'soil' | 'water' | 'flood_hazard' | 'elevation' | 'accessibility' | 'services' | 'roads' | 'parcels'
 
 interface DataTheme {
   id: DataThemeId
@@ -39,6 +40,20 @@ const DATA_THEMES: DataTheme[] = [
     description: 'Gebouwen en gebouwcontouren uit GRB of een andere persistente bron.',
     tokens: ['buildings', 'building', 'gebouwen', 'gebouw', 'bebouwing', 'gbg'],
   },
+  {
+    id: 'space_occupation',
+    label: 'Ruimtebeslag',
+    shortLabel: 'Ruimtebeslag',
+    description: 'Officiële 10 m-beleidskaart van ruimte ingenomen door wonen, economie, infrastructuur en recreatie.',
+    tokens: ['space_occupation', 'ruimtebeslag', 'ruibes'],
+  },
+  {
+    id: 'open_space',
+    label: 'Open ruimte',
+    shortLabel: 'Open ruimte',
+    description: 'Officiële 10 m-beleidskaart van open ruimte buiten kernen en ruimtebeslag.',
+    tokens: ['open_space', 'open ruimte', 'openruimte'],
+  },
   {
     id: 'population',
     label: 'Bevolking',
@@ -67,6 +82,13 @@ const DATA_THEMES: DataTheme[] = [
     description: 'Jaarlijkse officiële landbouwgebruikspercelen en hoofdteeltgroepen.',
     tokens: ['agriculture', 'agricultural', 'landbouw', 'landbouwgebruik', 'agpa'],
   },
+  {
+    id: 'soil',
+    label: 'Bodem',
+    shortLabel: 'Bodemkaart',
+    description: 'Historische DOV-bodemkartering met bodemtype, textuur en drainageklasse voor Mol.',
+    tokens: ['soil', 'bodem', 'bodemkaart', 'bodemtype', 'dov_soil_map'],
+  },
   {
     id: 'water',
     label: 'Water',
@@ -88,6 +110,20 @@ const DATA_THEMES: DataTheme[] = [
     description: 'Maaiveld- of oppervlaktehoogte, reliëf en helling uit DHMV II.',
     tokens: ['dhmv', 'elevation', 'height', 'hoogte', 'terrain', 'surface', 'dtm', 'dsm', 'reliëf'],
   },
+  {
+    id: 'accessibility',
+    label: 'Bereikbaarheid',
+    shortLabel: 'Knooppuntwaarde',
+    description: 'Knooppuntwaarde van collectief vervoer per hectare voor referentiejaar 2022.',
+    tokens: ['accessibility', 'bereikbaarheid', 'knooppuntwaarde', 'knptw'],
+  },
+  {
+    id: 'services',
+    label: 'Voorzieningen',
+    shortLabel: 'Voorzieningenniveau',
+    description: 'Genormaliseerde nabijheid van basis-, regionale en metropolitane voorzieningen in 2022.',
+    tokens: ['services', 'voorzieningen', 'voorzieningenniveau', 'totvznv'],
+  },
   {
     id: 'roads',
     label: 'Wegen',
@@ -106,13 +142,18 @@ const DATA_THEMES: DataTheme[] = [
 
 const DATA_THEME_MAP_STYLES: Record<DataThemeId, { fill: string; line: string }> = {
   buildings: { fill: '#d45f3d', line: '#9f3e24' },
+  space_occupation: { fill: '#be3e33', line: '#8f2c24' },
+  open_space: { fill: '#267a46', line: '#175c32' },
   population: { fill: '#7559a6', line: '#5b3f88' },
   forest: { fill: '#347950', line: '#225f3b' },
   nature_value: { fill: '#9a4f64', line: '#74364a' },
   agriculture: { fill: '#7b8f32', line: '#53671d' },
+  soil: { fill: '#9a7040', line: '#6f4c27' },
   water: { fill: '#2676a8', line: '#155b85' },
   flood_hazard: { fill: '#1597c2', line: '#075985' },
   elevation: { fill: '#a57a4b', line: '#315f59' },
+  accessibility: { fill: '#0f766e', line: '#115e59' },
+  services: { fill: '#b66d16', line: '#854d0e' },
   roads: { fill: '#6b7280', line: '#4b5563' },
   parcels: { fill: '#a7792f', line: '#7d571f' },
 }
@@ -126,6 +167,12 @@ function datasetAvailabilityLabel(dataset: DatasetCreateResponse): string {
     const resolution = Number(dataset.source_metadata?.['analysis_resolution_m'])
     return `${Number.isFinite(resolution) ? `${resolution.toLocaleString('nl-BE')} m` : 'Raster'} overstromingsscenario`
   }
+  if (dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster') {
+    const resolution = Number(dataset.source_metadata?.['analysis_resolution_m'])
+    const year = Number(dataset.source_metadata?.['observation_year'])
+    const resolutionLabel = Number.isFinite(resolution) ? `${resolution.toLocaleString('nl-BE')} m` : 'Raster'
+    return `${resolutionLabel} officiële bron${Number.isFinite(year) ? ` · ${year}` : ''}`
+  }
   return `${(dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0).toLocaleString('nl-BE')} objecten beschikbaar`
 }
 
@@ -156,6 +203,12 @@ function datasetMatchesTheme(dataset: DatasetCreateResponse, theme: DataTheme):
   if (dataset.source_name === 'digitaal_vlaanderen_dhmv') {
     return theme.id === 'elevation'
   }
+  if (dataset.source_name === 'department_omgeving_thematic_raster') {
+    return dataset.source_metadata?.['theme'] === theme.id
+  }
+  if (dataset.source_name === 'dov_soil_map') {
+    return theme.id === 'soil'
+  }
   const searchText = datasetSearchText(dataset)
   return theme.tokens.some((token) => searchText.includes(token))
 }
@@ -184,6 +237,7 @@ function pickThemeDataset(
       (dataset.source_name === 'department_omgeving_land_use' ? 90_000 : 0) +
       (dataset.source_name === 'inbo_bwk_natura2000' ? 95_000 : 0) +
       (dataset.source_name === 'agentschap_landbouw_zeevisserij_agricultural_parcels' ? 98_000 : 0) +
+      (dataset.source_name === 'department_omgeving_thematic_raster' ? 5_000_000 : 0) +
       (dataset.source_name === 'digitaal_vlaanderen_buildings_addresses_register' ? 5_000_000 : 0) +
       (dataset.source_name === 'digitaal_vlaanderen_dhmv' ? 5_000_000 : 0) +
       (dataset.source_name === 'vmm_flood_hazard' ? 5_000_000 : 0) +
@@ -763,7 +817,20 @@ export function MapWorkspace({
         opacity: 0.82,
       }
     : null
-  const activeImageOverlay = floodHazardImageOverlay ?? terrainImageOverlay ?? orthophotoImageOverlay
+  const thematicRasterBounds = activeThemeDataset?.source_name === 'department_omgeving_thematic_raster'
+    ? activeThemeDataset.source_metadata?.['bbox_epsg4326']
+    : null
+  const thematicRasterImageOverlay = activeThemeDataset?.source_name === 'department_omgeving_thematic_raster' && selectedProjectId && Array.isArray(thematicRasterBounds) && thematicRasterBounds.length === 4
+    ? {
+        url: thematicRasterImageUrl(selectedProjectId, activeThemeDataset.id),
+        bbox: thematicRasterBounds.map(Number) as [number, number, number, number],
+        label: getDatasetDisplayName(activeThemeDataset),
+        opacity: 0.78,
+      }
+    : null
+  const thematicLegendMin = String(activeThemeDataset?.source_metadata?.['legend_min_label'] ?? 'Lagere waarde')
+  const thematicLegendMax = String(activeThemeDataset?.source_metadata?.['legend_max_label'] ?? 'Hogere waarde')
+  const activeImageOverlay = thematicRasterImageOverlay ?? floodHazardImageOverlay ?? terrainImageOverlay ?? orthophotoImageOverlay
   const activeScopeProject = projects.find((project) => project.id === selectedProjectId) ?? null
   const activeScopeLabel = activeScopeProject ? operationalScopeProjectLabel(activeScopeProject) : 'Werkgebied'
   const municipalityAreaCount = areas.filter((area) => /^Gemeente\s/i.test(area.name)).length
@@ -775,7 +842,9 @@ export function MapWorkspace({
     [availableMapDatasets],
   )
   const activeTemporalSeriesGroups = themeTemporalSeriesMap[activeTheme.id]
-  const availableEvolutionThemes = DATA_THEMES.filter((theme) => themeTemporalSeriesMap[theme.id].length > 0)
+  const availableEvolutionThemes = DATA_THEMES.filter((theme) =>
+    themeTemporalSeriesMap[theme.id].some((group) => group.items.length >= 2),
+  )
   const activeTemporalSeriesGroup = activeTemporalSeriesGroups.find((group) => group.key === selectedTemporalSeriesKey)
     ?? activeTemporalSeriesGroups[0]
   const activeTemporalSeries = activeTemporalSeriesGroup?.items ?? EMPTY_TEMPORAL_SERIES
@@ -815,13 +884,28 @@ export function MapWorkspace({
     (metric) => metric.metric_key !== activeSelectionResult?.summary?.primary_metric_key,
   )
   const terrainReliefMetric = activeSupportingMetrics.find((metric) => metric.metric_key === 'relief_m')
+  const populationDensityMetric = activeSupportingMetrics.find((metric) => metric.metric_key === 'population_density_mean_per_ha')
+  const scoreMedianMetric = activeSupportingMetrics.find((metric) => metric.metric_key.endsWith('_median'))
   const activeSecondaryMetric = activeMetricUnit === 'm TAW'
     ? terrainReliefMetric ? `${terrainReliefMetric.metric_value.toLocaleString('nl-BE', { maximumFractionDigits: 2 })} m reliëf` : null
-    : selectedAreaSquareMetres && selectedAreaSquareMetres > 0
+    : activeMetricUnit === 'inwoners'
-      ? activeMetricUnit === 'ha'
+      ? populationDensityMetric ? selectionMetricLabel(populationDensityMetric) : null
+      : activeMetricUnit.startsWith('score')
+        ? scoreMedianMetric ? selectionMetricLabel(scoreMedianMetric) : null
+        : selectedAreaSquareMetres && selectedAreaSquareMetres > 0
+          ? activeMetricUnit === 'ha'
       ? `${((activeMetricValue * 10_000) / selectedAreaSquareMetres * 100).toLocaleString('nl-BE', { maximumFractionDigits: 1 })}% dekking`
       : `${(activeMetricValue / (selectedAreaSquareMetres / 1_000_000)).toLocaleString('nl-BE', { maximumFractionDigits: 1 })} ${activeMetricUnit} / km2`
       : null
+  const activeSecondaryLabel = activeMetricUnit === 'ha'
+    ? 'Aandeel selectie'
+    : activeMetricUnit === 'm TAW'
+      ? 'Reliëf'
+      : activeMetricUnit === 'inwoners'
+        ? 'Gemiddelde dichtheid'
+        : activeMetricUnit.startsWith('score')
+          ? 'Mediaan'
+          : 'Dichtheid'
   const selectedResultProperties = useMemo(() => {
     const keys = new Map<string, Set<string>>()
     for (const feature of activeSelectionResult?.geojson.features ?? []) {
@@ -1185,7 +1269,7 @@ export function MapWorkspace({
                 const dataset = themeDatasetMap[theme.id]
                 const temporalGroups = themeTemporalSeriesMap[theme.id]
                 const temporalGroup = temporalGroups[0]
-                const evolutionAvailable = temporalGroups.length > 0
+                const evolutionAvailable = temporalGroups.some((group) => group.items.length >= 2)
                 const available = Boolean(dataset) && (analysisMode === 'current' || evolutionAvailable)
                 const active = activeThemeId === theme.id
                 const firstObservation = temporalGroup?.items[0]?.observed_at
@@ -1382,7 +1466,12 @@ export function MapWorkspace({
               />
               <div className="geo-map-legend" aria-label="Kaartlegende">
                 <span><i className="geo-legend-area" /> Werkgebied</span>
-                {activeImageOverlay ? <span><i className="geo-legend-imagery" /> {activeImageOverlay.label}</span> : null}
+                {thematicRasterImageOverlay ? (
+                  <span className="geo-legend-thematic">
+                    <i className={`geo-legend-ramp geo-legend-ramp-${activeTheme.id}`} />
+                    <small>{thematicLegendMin} → {thematicLegendMax}</small>
+                  </span>
+                ) : activeImageOverlay ? <span><i className="geo-legend-imagery" /> {activeImageOverlay.label}</span> : null}
                 {analysisOverlayActive ? (
                   <>
                     <span><i className="geo-legend-layer geo-legend-layer-buildings" /> AI-kandidaten</span>
@@ -1585,7 +1674,7 @@ export function MapWorkspace({
                         <strong>{activeSelectionResult ? resultMetricLabel(activeSelectionResult) : 'Geen resultaat'}</strong>
                       </div>
                       <div>
-                        <span>{activeMetricUnit === 'ha' ? 'Aandeel selectie' : activeMetricUnit === 'm TAW' ? 'Reliëf' : 'Dichtheid'}</span>
+                        <span>{activeSecondaryLabel}</span>
                         <strong>{activeSecondaryMetric ?? (selectedDensity === null ? 'n.v.t.' : `${selectedDensity.toLocaleString('nl-BE', { maximumFractionDigits: 1 })} / km2`)}</strong>
                       </div>
                     </div>
@@ -1658,7 +1747,7 @@ export function MapWorkspace({
 
                 {analysisMode === 'current' ? (
                   <div className="geo-result-actions">
-                    <button className="secondary-action" disabled={!activeSelectionResult || activeTheme.id === 'elevation' || activeTheme.id === 'flood_hazard'} type="button" onClick={downloadActiveThemeSelection}>Download GeoJSON</button>
+                    <button className="secondary-action" disabled={!activeSelectionResult || activeThemeDataset?.dataset_type === 'raster'} type="button" onClick={downloadActiveThemeSelection}>Download GeoJSON</button>
                     <button className="secondary-action" disabled={!activeSelectionResult} type="button" onClick={copyActiveThemeSelection}>Kopieer gegevens</button>
                   </div>
                 ) : null}

frontend/src/hooks/useMapSelectionExtract.ts

@@ -4,6 +4,7 @@ import { formatError } from '../lib/formatError'
 import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types'
 import { terrainSelectionToMapSelection } from '../lib/terrainSelection'
 import { floodHazardSelectionToMapSelection } from '../lib/floodHazardSelection'
+import { thematicRasterSelectionToMapSelection } from '../lib/thematicRaster'
 
 interface MapSelectionExtractOptions {
   selectedProjectId: string | null
@@ -42,8 +43,9 @@ export function useMapSelectionExtract({
     }
     const terrainDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'digitaal_vlaanderen_dhmv'
     const floodHazardDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'vmm_flood_hazard'
-    if (!isVectorDatasetType(selectedDataset.dataset_type) && !terrainDataset && !floodHazardDataset) {
+    const thematicRasterDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'department_omgeving_thematic_raster'
-      setMapSelectionError('Gebiedsanalyse ondersteunt een vectorlaag, DHMV-hoogtemodel of beheerd VMM-overstromingsscenario.')
+    if (!isVectorDatasetType(selectedDataset.dataset_type) && !terrainDataset && !floodHazardDataset && !thematicRasterDataset) {
+      setMapSelectionError('Gebiedsanalyse ondersteunt een vectorlaag of een beheerd thematisch raster.')
       return null
     }
 
@@ -63,6 +65,11 @@ export function useMapSelectionExtract({
               bbox: { ...bbox, crs: 'EPSG:4326' },
               area_id: areaId,
             }))
+          : thematicRasterDataset
+            ? thematicRasterSelectionToMapSelection(await datasetsApi.selectThematicRaster(selectedProjectId, selectedDataset.id, {
+                bbox: { ...bbox, crs: 'EPSG:4326' },
+                area_id: areaId,
+              }))
         : await datasetsApi.selectVectorFeatures(selectedProjectId, selectedDataset.id, {
             bbox: { ...bbox, crs: 'EPSG:4326' },
             area_id: areaId,

frontend/src/hooks/useMapThemeSelectionInsights.ts

@@ -4,6 +4,7 @@ import { datasetsApi } from '../services/api/datasets'
 import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types'
 import { terrainSelectionToMapSelection } from '../lib/terrainSelection'
 import { floodHazardSelectionToMapSelection } from '../lib/floodHazardSelection'
+import { thematicRasterSelectionToMapSelection } from '../lib/thematicRaster'
 
 export interface MapThemeQuery<TThemeId extends string> {
   themeId: TThemeId
@@ -66,6 +67,11 @@ export function useMapThemeSelectionInsights<TThemeId extends string>(
                   bbox,
                   area_id: areaId,
                 }))
+              : dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster'
+                ? thematicRasterSelectionToMapSelection(await datasetsApi.selectThematicRaster(selectedProjectId, dataset.id, {
+                    bbox,
+                    area_id: areaId,
+                  }))
             : await datasetsApi.selectVectorFeatures(selectedProjectId, dataset.id, {
                 bbox,
                 area_id: areaId,

frontend/src/lib/datasetDisplay.ts

@@ -9,8 +9,13 @@ const DATASET_LABEL_BY_LAYER: Record<string, string> = {
   forest: 'Bos en groen',
   nature_value: 'Natuurwaarde',
   agriculture: 'Landbouwgebruikspercelen',
+  soil: 'Digitale bodemkaart',
   flood_hazard: 'Overstromingsgevaar',
   elevation: 'Hoogte en reliëf',
+  space_occupation: 'Ruimtebeslag',
+  open_space: 'Open ruimte',
+  accessibility: 'Knooppuntwaarde',
+  services: 'Voorzieningenniveau',
   building_registry: 'Gebouwenregister',
   regional_boundary: 'Grens vervoerregio Kempen',
   municipality_boundaries: 'Gemeentegrenzen Kempen',
@@ -29,6 +34,8 @@ const DATASET_SOURCE_LABELS: Record<string, string> = {
   vrbg: 'Digitaal Vlaanderen',
   waterinfo: 'Waterinfo Vlaanderen',
   vmm_flood_hazard: 'Vlaamse Milieumaatschappij',
+  department_omgeving_thematic_raster: 'Departement Omgeving',
+  dov_soil_map: 'Databank Ondergrond Vlaanderen',
 }
 
 export function getDatasetSourceDisplayName(dataset: DatasetCreateResponse): string {
@@ -45,6 +52,10 @@ export function getDatasetDisplayName(dataset: DatasetCreateResponse): string {
     const productName = dataset.source_metadata?.['product_display_name']
     return typeof productName === 'string' && productName.trim() ? productName : 'VMM-overstromingsscenario'
   }
+  if (dataset.source_name === 'department_omgeving_thematic_raster') {
+    const productName = dataset.source_metadata?.['product_display_name']
+    return typeof productName === 'string' && productName.trim() ? productName : 'Officieel Vlaams themaraster'
+  }
   const layer = (dataset.reference_layer_name ?? dataset.source_metadata?.layer_name ?? dataset.source_metadata?.layer_type ?? '')
     .toString()
     .toLowerCase()

frontend/src/lib/sourcePortfolio.ts

@@ -35,6 +35,9 @@ const sourceNameIs = (dataset: DatasetCreateResponse, sourceName: string): boole
 const referenceLayerIs = (dataset: DatasetCreateResponse, layerName: string): boolean =>
   String(dataset.reference_layer_name ?? dataset.source_metadata?.['theme'] ?? '').toLowerCase() === layerName
 
+const thematicProductIs = (dataset: DatasetCreateResponse, productKey: string): boolean =>
+  sourceNameIs(dataset, 'department_omgeving_thematic_raster') && dataset.source_metadata?.['product_key'] === productKey
+
 export const SOURCE_DOMAINS: SourceDomainDefinition[] = [
   {
     key: 'space',
@@ -127,7 +130,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     metricExamples: 'hectare ruimtebeslag, aandeel, vergelijking met open ruimte',
     priority: 'next',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/ruimtebeslag-vlaanderen-toestand-2025',
-    matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_space_occupation'),
+    matches: (dataset) => thematicProductIs(dataset, 'space_occupation_2025'),
   },
   {
     key: 'settlement_typology',
@@ -137,7 +140,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     coverage: 'Toestand 2022',
     value: 'Begrijpelijke morfologie voor kern, lint en verspreide bebouwing.',
     metricExamples: 'oppervlakte en aandeel per morfologisch type',
-    priority: 'planned',
+    priority: 'next',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/kernen-linten-verspreide-bebouwing-in-vlaanderen-kernen-toestand-2022',
     matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_settlement_typology'),
   },
@@ -187,7 +190,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     metricExamples: 'hectare open ruimte, aandeel en fragmentatie',
     priority: 'planned',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/open-ruimte-vlaanderen-toestand-2022',
-    matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_open_space'),
+    matches: (dataset) => thematicProductIs(dataset, 'open_space_2022'),
   },
   {
     key: 'dhmv',
@@ -244,10 +247,10 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     owner: 'Departement Omgeving',
     coverage: '1 ha raster, toestand 2022',
     value: 'Modelmatige bereikbaarheid via collectief vervoer.',
-    metricExamples: 'gemiddelde score en aandeel per bereikbaarheidsklasse',
+    metricExamples: 'gemiddelde, mediaan en percentielen van de bronindex',
     priority: 'next',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/knooppuntwaarde-per-ha-toestand-2022',
-    matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_node_value'),
+    matches: (dataset) => thematicProductIs(dataset, 'node_value_2022'),
   },
   {
     key: 'hoppin',
@@ -295,7 +298,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     metricExamples: 'geschat aantal inwoners en inwoners per hectare',
     priority: 'next',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/inwonersdichtheid-per-ha-vlaanderen-toestand-2019',
-    matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_population_density'),
+    matches: (dataset) => thematicProductIs(dataset, 'population_density_2019'),
   },
   {
     key: 'service_level',
@@ -304,10 +307,10 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [
     owner: 'Departement Omgeving',
     coverage: '1 ha raster, toestand 2022',
     value: 'Nabijheid van dagelijkse voorzieningen in één brongetrouwe score.',
-    metricExamples: 'gemiddelde score en aandeel per voorzieningsklasse',
+    metricExamples: 'gemiddelde, mediaan en percentielen van de 0-1-score',
     priority: 'next',
     url: 'https://www.vlaanderen.be/datavindplaats/catalogus/totaal-voorzieningenniveau-toestand-2022',
-    matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_service_level'),
+    matches: (dataset) => thematicProductIs(dataset, 'service_level_2022'),
   },
   {
     key: 'municipality_indicators',

frontend/src/lib/thematicRaster.ts

@@ -0,0 +1,27 @@
+import type { ThematicRasterSelectionResponse, VectorSelectionResponse } from '../types'
+
+export function thematicRasterImageUrl(projectId: string, datasetId: string): string {
+  return `/api/v1/projects/${projectId}/datasets/${datasetId}/raster/thematic/image`
+}
+
+export function thematicRasterSelectionToMapSelection(result: ThematicRasterSelectionResponse): VectorSelectionResponse {
+  return {
+    selection_bbox: result.selection_bbox,
+    selection_area_id: result.selection_area_id,
+    feature_count: result.valid_cell_count,
+    total_feature_count: result.valid_cell_count,
+    limit: 0,
+    truncated: false,
+    geojson: { type: 'FeatureCollection', features: [] },
+    summary: {
+      ...result.summary,
+      feature_count: result.valid_cell_count,
+      is_estimate: true,
+      warning: result.limitation_message,
+      metrics: result.summary.metrics.map((metric) => ({
+        ...metric,
+        is_estimate: metric.is_estimate ?? true,
+      })),
+    },
+  }
+}

frontend/src/services/api/datasets.ts

@@ -24,6 +24,9 @@ import type {
   FloodHazardSelectionResponse,
   DhmvProductRead,
   TerrainSelectionResponse,
+  ThematicRasterAcquireRequest,
+  ThematicRasterProductRead,
+  ThematicRasterSelectionResponse,
 } from '../../types'
 
 const DATASET_PAGE_SIZE = 200
@@ -141,6 +144,16 @@ export const datasetsApi = {
     payload: { bbox: VectorSelectionRequest['bbox']; area_id?: string },
   ): Promise<FloodHazardSelectionResponse> =>
     apiPost<FloodHazardSelectionResponse>(`/api/v1/projects/${projectId}/datasets/${datasetId}/raster/flood-hazard/select`, payload),
+  acquireThematicRaster: (projectId: string, payload: ThematicRasterAcquireRequest): Promise<JobRead> =>
+    apiPost<JobRead>(`/api/v1/projects/${projectId}/datasets/thematic-raster/acquire`, payload),
+  listThematicRasterProducts: (projectId: string): Promise<{ items: ThematicRasterProductRead[]; total: number }> =>
+    apiGet<{ items: ThematicRasterProductRead[]; total: number }>(`/api/v1/projects/${projectId}/datasets/thematic-raster/products`),
+  selectThematicRaster: (
+    projectId: string,
+    datasetId: string,
+    payload: { bbox: VectorSelectionRequest['bbox']; area_id?: string },
+  ): Promise<ThematicRasterSelectionResponse> =>
+    apiPost<ThematicRasterSelectionResponse>(`/api/v1/projects/${projectId}/datasets/${datasetId}/raster/thematic/select`, payload),
   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

@@ -5751,13 +5751,18 @@ section {
 }
 
 .geo-theme-symbol-buildings { background: #d45f3d; }
+.geo-theme-symbol-space_occupation { background: #be3e33; }
+.geo-theme-symbol-open_space { background: #267a46; }
 .geo-theme-symbol-population { background: #7559a6; }
 .geo-theme-symbol-forest { background: #347950; }
 .geo-theme-symbol-nature_value { background: #9a4f64; }
 .geo-theme-symbol-agriculture { background: #7b8f32; }
+.geo-theme-symbol-soil { background: #9a7040; }
 .geo-theme-symbol-water { background: #2676a8; }
 .geo-theme-symbol-flood_hazard { background: #1597c2; }
 .geo-theme-symbol-elevation { background: #a57a4b; }
+.geo-theme-symbol-accessibility { background: #0f766e; }
+.geo-theme-symbol-services { background: #b66d16; }
 .geo-theme-symbol-roads { background: #6b7280; }
 .geo-theme-symbol-parcels { background: #a7792f; }
 
@@ -5932,6 +5937,16 @@ section {
   background: rgba(117, 89, 166, 0.24);
 }
 
+.geo-map-legend .geo-legend-layer-space_occupation {
+  border-color: #8f2c24;
+  background: rgba(190, 62, 51, 0.28);
+}
+
+.geo-map-legend .geo-legend-layer-open_space {
+  border-color: #175c32;
+  background: rgba(38, 122, 70, 0.26);
+}
+
 .geo-map-legend .geo-legend-layer-forest {
   border-color: #225f3b;
   background: rgba(52, 121, 80, 0.24);
@@ -5942,6 +5957,11 @@ section {
   background: rgba(38, 118, 168, 0.24);
 }
 
+.geo-map-legend .geo-legend-layer-soil {
+  border-color: #6f4c27;
+  background: rgba(154, 112, 64, 0.24);
+}
+
 .geo-map-legend .geo-legend-layer-flood_hazard {
   border-color: #075985;
   background: rgba(21, 151, 194, 0.28);
@@ -5952,6 +5972,16 @@ section {
   background: rgba(165, 122, 75, 0.26);
 }
 
+.geo-map-legend .geo-legend-layer-accessibility {
+  border-color: #115e59;
+  background: rgba(15, 118, 110, 0.26);
+}
+
+.geo-map-legend .geo-legend-layer-services {
+  border-color: #854d0e;
+  background: rgba(182, 109, 22, 0.26);
+}
+
 .geo-map-legend .geo-legend-layer-roads {
   border-color: #4b5563;
   background: rgba(107, 114, 128, 0.24);
@@ -5972,6 +6002,32 @@ section {
   background: linear-gradient(135deg, #7a9b68 0 33%, #d1b37a 33% 66%, #8eb5cb 66%);
 }
 
+.geo-map-legend .geo-legend-thematic {
+  display: inline-flex;
+  align-items: center;
+  gap: 0.42rem;
+}
+
+.geo-map-legend .geo-legend-thematic small {
+  color: #42504c;
+  font-size: 0.66rem;
+  font-weight: 700;
+}
+
+.geo-map-legend .geo-legend-ramp {
+  display: block;
+  width: 3.4rem;
+  height: 0.62rem;
+  border: 1px solid rgba(23, 39, 34, 0.25);
+  border-radius: 2px;
+}
+
+.geo-legend-ramp-space_occupation { background: linear-gradient(90deg, #fbe7d3, #be3e33); }
+.geo-legend-ramp-open_space { background: linear-gradient(90deg, #ddeedb, #267a46); }
+.geo-legend-ramp-population { background: linear-gradient(90deg, #eee7f6, #673a97); }
+.geo-legend-ramp-accessibility { background: linear-gradient(90deg, #e9f1f4, #0f766e); }
+.geo-legend-ramp-services { background: linear-gradient(90deg, #fff4bf, #b66d16); }
+
 .geo-map-legend .geo-legend-added {
   border-color: #15803d;
   background: rgba(22, 163, 74, 0.2);

frontend/src/types.ts

@@ -434,6 +434,57 @@ export interface FloodHazardSelectionResponse {
   generated_at: string
 }
 
+export interface ThematicRasterAcquireRequest {
+  bbox: VectorSelectionBBox
+  area_id?: string | null
+  product_key: string
+  force_refresh?: boolean
+}
+
+export interface ThematicRasterProductRead {
+  key: string
+  display_name: string
+  theme: 'space_occupation' | 'open_space' | 'population' | 'accessibility' | 'services'
+  metric_kind: 'binary_area' | 'population_density' | 'index_score' | 'normalized_score'
+  coverage_id: string
+  native_resolution_m: number
+  source_crs: 'EPSG:31370'
+  source_value_unit: string
+  observation_year: number
+  source_version: string
+  catalog_url: string
+  attribution: string
+  license_note: string
+  legend_min_label: string
+  legend_max_label: string
+  limitation_message: string
+}
+
+export interface ThematicRasterSelectionResponse {
+  dataset_id: string
+  product_key: string
+  theme: ThematicRasterProductRead['theme']
+  metric_kind: ThematicRasterProductRead['metric_kind']
+  selection_bbox: VectorSelectionBBox
+  selection_area_id?: string | null
+  selected_cell_count: number
+  valid_cell_count: number
+  coverage_ratio: number
+  resolution_m: number
+  observation_year: number
+  summary: {
+    metric_label: string
+    metric_value: number
+    metric_unit: string
+    aggregation_method: string
+    primary_metric_key: string
+    metrics: Array<VectorSelectionMetric & { is_estimate: boolean }>
+  }
+  unsupported_metrics: string[]
+  limitation_message: string
+  generated_at: string
+}
+
 export interface MapImageOverlay {
   url: string
   bbox: [number, number, number, number]

scripts/README.md

@@ -1613,6 +1613,33 @@ take a long time because every VMM WCS tile is bounded, rate-limited and
 validated. This is expected operator work; the app never fetches these rasters
 on page load or map click.
 
+## Cross-domain Mol profile
+
+Load the five official policy rasters for the exact Mol municipality Area and
+immediately verify each persisted selection result:
+
+```bash
+docker exec geointel python /app/scripts/provision_thematic_rasters.py
+```
+
+Plan the later complete Kempen rollout without source fetches or writes:
+
+```bash
+docker exec geointel python /app/scripts/provision_thematic_rasters.py \
+  --project-name "Kempen Regional Workbench" --all-municipalities --dry-run
+```
+
+Load the official DOV soil map for Mol:
+
+```bash
+docker exec geointel python /app/scripts/provision_mol_soil_map.py
+```
+
+`--fetch-only` builds the soil artifact and manifest without API import;
+`--force` is the only way to bypass an existing ready soil Dataset. Both
+operators use canonical APIs and persistent operator-evidence storage. They do
+not run on application startup.
+
 ## Tower deployment
 
 Push the local branch to Gitea, then rebuild the Unraid/Tower Docker runtime:

scripts/provision_mol_soil_map.py

@@ -0,0 +1,661 @@
+"""Provision the official DOV digital soil map for the municipality of Mol.
+
+The operator follows every bounded WFS page, retains checksummed source
+responses, clips soil polygons to the persisted Mol Area in EPSG:31370 and
+imports the result through GeoIntel's canonical dataset upload route. It does
+not write directly to vector_features and it does not treat the historical
+1949-1971 field survey as a current drainage observation.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import os
+import sys
+from collections import defaultdict
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any, Iterable
+
+import requests
+from pyproj import Transformer
+from requests.adapters import HTTPAdapter
+from shapely.geometry import MultiPolygon, Polygon, mapping, shape
+from shapely.ops import transform as transform_geometry
+from shapely.ops import unary_union
+from shapely.validation import make_valid
+from urllib3.util.retry import Retry
+
+
+WFS_URL = "https://www.dov.vlaanderen.be/geoserver/wfs"
+CATALOG_URL = (
+    "https://www.vlaanderen.be/datavindplaats/catalogus/"
+    "digitale-bodemkaart-van-het-vlaams-gewest-bodemtypes"
+)
+TYPE_NAME = "bodemkaart:bodemtypes"
+SOURCE_NAME = "dov_soil_map"
+SOURCE_VERSION = "Digitale uitgave juni 2017"
+SURVEY_PERIOD = "1949-1971"
+OBSERVED_AT = "1971-12-31T23:59:59Z"
+VALID_FROM = "1949-01-01T00:00:00Z"
+VALID_TO = OBSERVED_AT
+ATTRIBUTION = "Databank Ondergrond Vlaanderen - Digitale bodemkaart: bodemtypes"
+DEFAULT_API_URL = "http://127.0.0.1:8000"
+DEFAULT_OUTPUT_DIR = "/app/storage/operator-evidence/dov-soil-map/mol"
+DEFAULT_PROJECT_NAME = "Mol Municipality Workbench"
+DEFAULT_AREA_FRAGMENT = "Gemeente Mol"
+DATASET_FILENAME = "dov_soil_map_mol.geojson"
+MANIFEST_FILENAME = "dov_soil_map_mol.manifest.json"
+SCHEMA_VERSION = 1
+
+TO_LAMBERT72 = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True)
+TO_WGS84 = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True)
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Provision the official DOV soil map for Mol.")
+    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-fragment", default=DEFAULT_AREA_FRAGMENT)
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path(os.environ.get("GEOINTEL_SOIL_MAP_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)),
+    )
+    parser.add_argument("--page-limit", type=int, default=500)
+    parser.add_argument("--max-features", type=int, default=20_000)
+    parser.add_argument("--request-timeout", type=int, default=180)
+    parser.add_argument("--import-timeout", type=int, default=1800)
+    parser.add_argument("--force", action="store_true")
+    parser.add_argument("--fetch-only", action="store_true")
+    return parser.parse_args()
+
+
+def utc_now() -> str:
+    return datetime.now(timezone.utc).isoformat()
+
+
+def sha256_bytes(value: bytes) -> str:
+    return hashlib.sha256(value).hexdigest()
+
+
+def sha256_file(path: Path) -> str:
+    digest = hashlib.sha256()
+    with path.open("rb") as handle:
+        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
+            digest.update(chunk)
+    return digest.hexdigest()
+
+
+def write_bytes_atomic(path: Path, value: bytes) -> None:
+    path.parent.mkdir(parents=True, exist_ok=True)
+    temporary = path.with_suffix(path.suffix + ".tmp")
+    temporary.write_bytes(value)
+    temporary.replace(path)
+
+
+def write_json_atomic(path: Path, value: Any, *, pretty: bool = False) -> None:
+    encoded = json.dumps(
+        value,
+        ensure_ascii=False,
+        indent=2 if pretty else None,
+        separators=None if pretty else (",", ":"),
+    ).encode("utf-8")
+    write_bytes_atomic(path, encoded)
+
+
+def source_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-DOV-Soil-Mol-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 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 = list(page.get("items") or [])
+        page_total = int(page.get("total") or 0)
+        if total is None:
+            total = page_total
+        elif page_total != total:
+            raise RuntimeError("GeoIntel pagination total changed while locating the Mol workspace")
+        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 pagination returned {len(items)} of {total} items")
+    return items
+
+
+def polygonal_geometry(geometry):
+    if geometry is None or geometry.is_empty:
+        return None
+    if not geometry.is_valid:
+        geometry = make_valid(geometry)
+    polygons: list[Polygon] = []
+
+    def collect(candidate) -> None:
+        if candidate is None or candidate.is_empty:
+            return
+        if isinstance(candidate, Polygon):
+            polygons.append(candidate)
+        elif isinstance(candidate, MultiPolygon):
+            polygons.extend(part for part in candidate.geoms if not part.is_empty)
+        elif hasattr(candidate, "geoms"):
+            for part in candidate.geoms:
+                collect(part)
+
+    collect(geometry)
+    if not polygons:
+        return None
+    result = unary_union(polygons)
+    if not result.is_valid:
+        result = make_valid(result)
+    return result if not result.is_empty and result.is_valid else None
+
+
+def locate_workspace(
+    session: requests.Session,
+    base_url: str,
+    project_name: str,
+    area_fragment: str,
+    timeout: int,
+) -> tuple[str, str, Any, list[dict[str, Any]]]:
+    projects = paginated_items(session, f"{base_url}/api/v1/projects", timeout=timeout)
+    project = next((item for item in projects if item.get("name") == project_name), None)
+    if not project:
+        raise RuntimeError(f"Project {project_name!r} is missing")
+    project_id = str(project["id"])
+    areas = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", timeout=timeout)
+    area = next(
+        (item for item in areas if area_fragment.casefold() in str(item.get("name") or "").casefold()),
+        None,
+    )
+    if not area or not area.get("geometry"):
+        raise RuntimeError(f"Persisted Mol Area containing {area_fragment!r} is missing")
+    boundary_wgs84 = polygonal_geometry(shape(area["geometry"]))
+    if boundary_wgs84 is None:
+        raise RuntimeError("Persisted Mol Area is not valid polygonal geometry")
+    datasets = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/datasets", timeout=timeout)
+    return project_id, str(area["id"]), boundary_wgs84, datasets
+
+
+def iter_wfs_pages(
+    session: requests.Session,
+    bbox_lambert72: tuple[float, float, float, float],
+    *,
+    page_limit: int,
+    timeout: int,
+) -> Iterable[tuple[dict[str, Any], str, bytes]]:
+    start_index = 0
+    expected_total: int | None = None
+    while True:
+        params = {
+            "service": "WFS",
+            "version": "2.0.0",
+            "request": "GetFeature",
+            "typeNames": TYPE_NAME,
+            "srsName": "EPSG:4326",
+            "bbox": ",".join(f"{value:.3f}" for value in bbox_lambert72) + ",EPSG:31370",
+            "count": str(page_limit),
+            "startIndex": str(start_index),
+            "sortBy": "gid",
+            "outputFormat": "application/json",
+        }
+        response = session.get(WFS_URL, params=params, timeout=timeout)
+        response.raise_for_status()
+        payload = response.json()
+        if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection":
+            raise RuntimeError("DOV WFS returned an invalid FeatureCollection")
+        features = list(payload.get("features") or [])
+        matched = int(payload.get("numberMatched") or payload.get("totalFeatures") or 0)
+        if expected_total is None:
+            expected_total = matched
+        elif matched != expected_total:
+            raise RuntimeError("DOV WFS numberMatched changed during pagination")
+        yield payload, response.url, response.content
+        returned = int(payload.get("numberReturned") or len(features))
+        if returned != len(features):
+            raise RuntimeError("DOV WFS numberReturned does not match its feature payload")
+        start_index += returned
+        if returned == 0 or start_index >= expected_total:
+            if start_index != expected_total:
+                raise RuntimeError(f"DOV WFS returned {start_index} of {expected_total} matched features")
+            break
+
+
+def normalize_feature(feature: dict[str, Any], boundary_lambert72) -> tuple[dict[str, Any] | None, bool]:
+    geometry_payload = feature.get("geometry")
+    if not geometry_payload:
+        return None, False
+    source_wgs84 = polygonal_geometry(shape(geometry_payload))
+    if source_wgs84 is None:
+        return None, False
+    source_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, source_wgs84))
+    if source_lambert72 is None or not source_lambert72.intersects(boundary_lambert72):
+        return None, False
+    was_clipped = not source_lambert72.within(boundary_lambert72)
+    clipped_lambert72 = polygonal_geometry(source_lambert72.intersection(boundary_lambert72))
+    if clipped_lambert72 is None or clipped_lambert72.area <= 0:
+        return None, was_clipped
+    clipped_wgs84 = polygonal_geometry(transform_geometry(TO_WGS84.transform, clipped_lambert72))
+    if clipped_wgs84 is None:
+        return None, was_clipped
+
+    raw = dict(feature.get("properties") or {})
+    gid = raw.get("gid")
+    map_polygon_id = raw.get("id_kaartvlak")
+    source_id = str(feature.get("id") or f"{TYPE_NAME}:{gid or map_polygon_id}")
+    properties = {
+        "source_name": SOURCE_NAME,
+        "source_collection": TYPE_NAME,
+        "source_feature_id": source_id,
+        "source_gid": gid,
+        "source_map_polygon_id": map_polygon_id,
+        "reference_layer_name": "soil",
+        "theme": "soil",
+        "authority_level": "authoritative_historical_baseline",
+        "coverage_scope": "municipality",
+        "municipality": "Mol",
+        "nis_code": "13025",
+        "source_version": SOURCE_VERSION,
+        "survey_period": SURVEY_PERIOD,
+        "soil_type_code": raw.get("Bodemtype"),
+        "unified_soil_type_code": raw.get("Unibodemtype"),
+        "soil_series_code": raw.get("Bodemserie"),
+        "soil_series_description": raw.get("Beknopte_omschrijving_bodemserie"),
+        "soil_generalized_legend": raw.get("Gegeneraliseerde_legende"),
+        "soil_texture_class_code": raw.get("Textuurklasse_code"),
+        "soil_texture_class": raw.get("Textuurklasse"),
+        "soil_drainage_class_code": raw.get("Drainageklasse_code"),
+        "soil_drainage_class": raw.get("Drainageklasse"),
+        "soil_profile_group_code": raw.get("Profielontwikkelingsgroep_code"),
+        "soil_profile_group": raw.get("Profielontwikkelingsgroep"),
+        "soil_substrate_code": raw.get("Substraat_code"),
+        "soil_substrate": raw.get("Substraat_Vlaanderen") or raw.get("Substraat_legende"),
+        "soil_region": raw.get("Streek"),
+        "classification_type": raw.get("Type_classificatie"),
+        "soil_map_title": raw.get("Eenduidige_legende_titel"),
+        "clipped_area_ha": round(float(clipped_lambert72.area) / 10_000.0, 8),
+        "attribution": ATTRIBUTION,
+        "historical_drainage_limitation": (
+            "Drainage class derives from field data collected between 1949 and 1971 and may differ today."
+        ),
+    }
+    return {
+        "type": "Feature",
+        "id": source_id,
+        "geometry": mapping(clipped_wgs84),
+        "properties": properties,
+    }, was_clipped
+
+
+def prepare_artifact(
+    session: requests.Session,
+    boundary_wgs84,
+    output_dir: Path,
+    *,
+    page_limit: int,
+    max_features: int,
+    timeout: int,
+) -> tuple[Path, Path, dict[str, Any]]:
+    output_dir.mkdir(parents=True, exist_ok=True)
+    raw_dir = output_dir / "raw"
+    raw_dir.mkdir(parents=True, exist_ok=True)
+    boundary_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, boundary_wgs84))
+    if boundary_lambert72 is None:
+        raise RuntimeError("Mol boundary could not be transformed to EPSG:31370")
+
+    retained: list[dict[str, Any]] = []
+    raw_pages: list[dict[str, Any]] = []
+    source_urls: list[str] = []
+    seen_ids: set[str] = set()
+    raw_feature_count = 0
+    duplicate_count = 0
+    rejected_count = 0
+    clipped_count = 0
+    area_by_legend: dict[str, float] = defaultdict(float)
+    area_by_texture: dict[str, float] = defaultdict(float)
+    area_by_drainage: dict[str, float] = defaultdict(float)
+
+    for page_number, (payload, source_url, raw_bytes) in enumerate(
+        iter_wfs_pages(
+            session,
+            boundary_lambert72.bounds,
+            page_limit=page_limit,
+            timeout=timeout,
+        ),
+        start=1,
+    ):
+        page_path = raw_dir / f"dov_soil_map_page_{page_number:05d}.json"
+        write_bytes_atomic(page_path, raw_bytes)
+        features = list(payload.get("features") or [])
+        raw_feature_count += len(features)
+        if raw_feature_count > max_features:
+            raise RuntimeError(
+                f"DOV WFS exceeded the {max_features} feature safety limit; refusing a truncated import"
+            )
+        raw_pages.append(
+            {
+                "path": str(page_path.relative_to(output_dir)),
+                "sha256": sha256_bytes(raw_bytes),
+                "size_bytes": len(raw_bytes),
+                "feature_count": len(features),
+                "source_url": source_url,
+            }
+        )
+        source_urls.append(source_url)
+        for feature in features:
+            raw = dict(feature.get("properties") or {})
+            source_id = str(feature.get("id") or f"{TYPE_NAME}:{raw.get('gid')}")
+            if source_id in seen_ids:
+                duplicate_count += 1
+                continue
+            seen_ids.add(source_id)
+            normalized, was_clipped = normalize_feature(feature, boundary_lambert72)
+            if normalized is None:
+                rejected_count += 1
+                continue
+            if was_clipped:
+                clipped_count += 1
+            retained.append(normalized)
+            properties = normalized["properties"]
+            area = float(properties["clipped_area_ha"])
+            area_by_legend[str(properties.get("soil_generalized_legend") or "Onbekend")] += area
+            area_by_texture[str(properties.get("soil_texture_class") or "Onbekend")] += area
+            area_by_drainage[str(properties.get("soil_drainage_class") or "Onbekend")] += area
+
+    if not retained:
+        raise RuntimeError("DOV WFS returned no valid soil polygons inside the persisted Mol Area")
+    generated_at = utc_now()
+    artifact = {
+        "type": "FeatureCollection",
+        "name": "Digitale bodemkaart - Gemeente Mol",
+        "features": retained,
+        "source": ATTRIBUTION,
+        "source_version": SOURCE_VERSION,
+        "survey_period": SURVEY_PERIOD,
+        "catalog_url": CATALOG_URL,
+        "generated_at": generated_at,
+    }
+    artifact_path = output_dir / DATASET_FILENAME
+    write_json_atomic(artifact_path, artifact)
+    manifest = {
+        "schema_version": SCHEMA_VERSION,
+        "source_version": SOURCE_VERSION,
+        "survey_period": SURVEY_PERIOD,
+        "source_type_name": TYPE_NAME,
+        "wfs_url": WFS_URL,
+        "catalog_url": CATALOG_URL,
+        "attribution": ATTRIBUTION,
+        "generated_at": generated_at,
+        "crs_source_service": "EPSG:31370",
+        "crs_response_and_persisted": "EPSG:4326",
+        "crs_clip_and_area_measurement": "EPSG:31370",
+        "boundary_sha256": sha256_bytes(json.dumps(mapping(boundary_wgs84), sort_keys=True).encode("utf-8")),
+        "boundary_bbox_wgs84": list(boundary_wgs84.bounds),
+        "boundary_bbox_epsg31370": list(boundary_lambert72.bounds),
+        "page_limit": page_limit,
+        "page_count": len(raw_pages),
+        "raw_source_feature_count": raw_feature_count,
+        "feature_count": len(retained),
+        "duplicate_count": duplicate_count,
+        "rejected_or_outside_count": rejected_count,
+        "clipped_feature_count": clipped_count,
+        "reference_truncated": False,
+        "raw_pages": raw_pages,
+        "source_urls": source_urls,
+        "area_by_generalized_legend_ha": {key: round(value, 6) for key, value in sorted(area_by_legend.items())},
+        "area_by_texture_ha": {key: round(value, 6) for key, value in sorted(area_by_texture.items())},
+        "area_by_drainage_ha": {key: round(value, 6) for key, value in sorted(area_by_drainage.items())},
+        "artifact_path": str(artifact_path),
+        "artifact_sha256": sha256_file(artifact_path),
+        "artifact_size_bytes": artifact_path.stat().st_size,
+        "limitations": [
+            "The map is based on field data collected between 1949 and 1971.",
+            "Current drainage, land use and local soil disturbance may differ from the mapped class.",
+            "The 1:20,000 source is contextual evidence and not a parcel-scale soil investigation.",
+        ],
+    }
+    manifest_path = output_dir / MANIFEST_FILENAME
+    write_json_atomic(manifest_path, manifest, pretty=True)
+    return artifact_path, manifest_path, manifest
+
+
+def selection_metrics() -> list[dict[str, Any]]:
+    return [
+        {
+            "metric_key": "soil_dry_sand_area",
+            "method": "intersection_area",
+            "label": "Gekarteerd als droog zand",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "soil_generalized_legend",
+            "filter_values": ["Droog zand", "Zeer droog zand"],
+        },
+        {
+            "metric_key": "soil_moist_sand_area",
+            "method": "intersection_area",
+            "label": "Gekarteerd als vochtig zand",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "soil_generalized_legend",
+            "filter_values": ["Vochtig zand"],
+        },
+        {
+            "metric_key": "soil_wet_sand_area",
+            "method": "intersection_area",
+            "label": "Gekarteerd als nat zand",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "soil_generalized_legend",
+            "filter_values": ["Nat zand", "Zeer nat zand"],
+        },
+        {
+            "metric_key": "soil_anthropogenic_area",
+            "method": "intersection_area",
+            "label": "Antropogene bodemklasse",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "soil_generalized_legend",
+            "filter_values": ["Antropogeen"],
+        },
+    ]
+
+
+def upload_artifact(
+    session: requests.Session,
+    *,
+    base_url: str,
+    project_id: str,
+    area_id: str,
+    artifact_path: Path,
+    manifest_path: Path,
+    manifest: dict[str, Any],
+    timeout: int,
+) -> dict[str, Any]:
+    limitation = (
+        "Historische bodemkartering op schaal 1:20.000 op basis van veldwerk 1949-1971; "
+        "de huidige drainage en lokale bodemtoestand kunnen afwijken."
+    )
+    source_metadata = {
+        "provider": SOURCE_NAME,
+        "theme": "soil",
+        "layer_type": "soil",
+        "source_collection": TYPE_NAME,
+        "source_crs": "EPSG:31370",
+        "persisted_crs": "EPSG:4326",
+        "authority_level": "authoritative_historical_baseline",
+        "coverage_scope": "municipality",
+        "municipality": "Mol",
+        "nis_code": "13025",
+        "feature_count": manifest["feature_count"],
+        "geometry_clipped_to_area": True,
+        "semantic_metrics": False,
+        "survey_period": SURVEY_PERIOD,
+        "source_scale": "1:20,000",
+        "attribution": ATTRIBUTION,
+        "catalog_url": CATALOG_URL,
+        "license_note": "DOV standard attribution and public GDI reuse conditions apply.",
+        "limitation_message": limitation,
+        "selection_aggregation": {
+            "metric_key": "soil_mapped_area",
+            "method": "intersection_area",
+            "label": "Bodemkaartoppervlakte",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "warning": limitation,
+        },
+        "selection_metrics": selection_metrics(),
+    }
+    provenance_metadata = {
+        "operator_tool": "provision_mol_soil_map.py",
+        "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
+        "source_type_name": TYPE_NAME,
+        "wfs_url": WFS_URL,
+        "catalog_url": CATALOG_URL,
+        "manifest_path": str(manifest_path),
+        "artifact_sha256": manifest["artifact_sha256"],
+        "raw_page_checksums": {page["path"]: page["sha256"] for page in manifest["raw_pages"]},
+        "source_urls": manifest["source_urls"],
+        "reference_truncated": False,
+        "generated_at": manifest["generated_at"],
+        "limitations": manifest["limitations"],
+    }
+    with artifact_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": SOURCE_NAME,
+                "reference_layer_name": "soil",
+                "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": "dov:digital-soil-map:mol",
+                "observed_at": OBSERVED_AT,
+                "valid_from": VALID_FROM,
+                "valid_to": VALID_TO,
+                "temporal_granularity": "period",
+                "source_version": SOURCE_VERSION,
+            },
+            files={"file": (artifact_path.name, handle, "application/geo+json")},
+            timeout=timeout,
+        )
+    return response_data(response)
+
+
+def main() -> int:
+    args = parse_args()
+    if args.page_limit < 1 or args.page_limit > 2000 or args.max_features < args.page_limit:
+        print(json.dumps({"status": "error", "message": "Invalid page or feature safety limits"}), file=sys.stderr)
+        return 2
+    base_url = args.base_url.rstrip("/")
+    api_session = requests.Session()
+    try:
+        project_id, area_id, boundary, datasets = locate_workspace(
+            api_session,
+            base_url,
+            args.project_name,
+            args.area_fragment,
+            args.request_timeout,
+        )
+        existing = next(
+            (
+                item
+                for item in datasets
+                if item.get("source_name") == SOURCE_NAME
+                and str(item.get("area_id") or "") == area_id
+                and item.get("status") == "ready"
+            ),
+            None,
+        )
+        if existing and not args.force:
+            result = {
+                "status": "reused",
+                "project_id": project_id,
+                "area_id": area_id,
+                "dataset_id": existing["id"],
+                "feature_count": existing.get("feature_count"),
+            }
+        else:
+            artifact_path, manifest_path, manifest = prepare_artifact(
+                source_session(),
+                boundary,
+                args.output_dir,
+                page_limit=args.page_limit,
+                max_features=args.max_features,
+                timeout=args.request_timeout,
+            )
+            if args.fetch_only:
+                result = {
+                    "status": "prepared",
+                    "project_id": project_id,
+                    "area_id": area_id,
+                    "artifact_path": str(artifact_path),
+                    "feature_count": manifest["feature_count"],
+                }
+            else:
+                dataset = upload_artifact(
+                    api_session,
+                    base_url=base_url,
+                    project_id=project_id,
+                    area_id=area_id,
+                    artifact_path=artifact_path,
+                    manifest_path=manifest_path,
+                    manifest=manifest,
+                    timeout=args.import_timeout,
+                )
+                result = {
+                    "status": "imported",
+                    "project_id": project_id,
+                    "area_id": area_id,
+                    "dataset_id": dataset["id"],
+                    "feature_count": dataset.get("feature_count") or manifest["feature_count"],
+                    "artifact_path": str(artifact_path),
+                }
+        print(json.dumps(result, ensure_ascii=False, indent=2))
+        return 0
+    except (OSError, RuntimeError, requests.RequestException, ValueError) as exc:
+        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
+        return 1
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

scripts/provision_thematic_rasters.py

@@ -0,0 +1,182 @@
+"""Provision governed Flemish thematic rasters through the GeoIntel API.
+
+The safe default loads all five products for the official Mol municipality
+Area. Use --all-municipalities with an explicitly named regional project to
+load every persisted municipality Area. The operator never writes to PostGIS
+or storage directly and never accepts an arbitrary external service URL.
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import os
+from typing import Any
+
+import requests
+
+
+DEFAULT_API_URL = "http://127.0.0.1:8000"
+DEFAULT_PROJECT_NAME = "Mol Municipality Workbench"
+DEFAULT_AREA_FRAGMENT = "Gemeente Mol"
+DEFAULT_PRODUCTS = (
+    "space_occupation_2025",
+    "open_space_2022",
+    "population_density_2019",
+    "node_value_2022",
+    "service_level_2022",
+)
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Provision governed Flemish thematic raster products.")
+    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", default=DEFAULT_AREA_FRAGMENT, help="Case-insensitive Area name fragment.")
+    parser.add_argument("--products", default=",".join(DEFAULT_PRODUCTS))
+    parser.add_argument("--all-municipalities", action="store_true", help="Process every Area whose name starts with 'Gemeente '.")
+    parser.add_argument("--force-refresh", action="store_true")
+    parser.add_argument("--timeout", type=int, default=900)
+    parser.add_argument("--dry-run", action="store_true")
+    return parser.parse_args()
+
+
+def unwrap(response: requests.Response) -> Any:
+    try:
+        payload = response.json()
+    except ValueError as exc:
+        raise RuntimeError(f"GeoIntel returned non-JSON HTTP {response.status_code}: {response.text[:300]}") from exc
+    if not response.ok:
+        error = payload.get("error") if isinstance(payload, dict) else None
+        message = error.get("message") if isinstance(error, dict) else response.text[:300]
+        raise RuntimeError(f"GeoIntel HTTP {response.status_code}: {message}")
+    return payload.get("data") if isinstance(payload, dict) and "data" in payload else payload
+
+
+def paged_items(session: requests.Session, url: str, timeout: int) -> list[dict[str, Any]]:
+    items: list[dict[str, Any]] = []
+    offset = 0
+    while True:
+        separator = "&" if "?" in url else "?"
+        page = unwrap(session.get(f"{url}{separator}limit=200&offset={offset}", timeout=timeout))
+        rows = list(page.get("items") or [])
+        items.extend(rows)
+        total = int(page.get("total") or 0)
+        if not rows or len(items) >= total:
+            return items
+        offset += len(rows)
+
+
+def geometry_bbox(geometry: dict[str, Any]) -> dict[str, Any]:
+    points: list[tuple[float, float]] = []
+
+    def visit(value: Any) -> None:
+        if isinstance(value, list) and len(value) >= 2 and all(isinstance(item, (int, float)) for item in value[:2]):
+            points.append((float(value[0]), float(value[1])))
+            return
+        if isinstance(value, list):
+            for item in value:
+                visit(item)
+
+    visit(geometry.get("coordinates"))
+    if not points:
+        raise RuntimeError("Persisted Area geometry contains no coordinates")
+    return {
+        "min_x": min(point[0] for point in points),
+        "min_y": min(point[1] for point in points),
+        "max_x": max(point[0] for point in points),
+        "max_y": max(point[1] for point in points),
+        "crs": "EPSG:4326",
+    }
+
+
+def find_project(projects: list[dict[str, Any]], name: str) -> dict[str, Any]:
+    matches = [project for project in projects if str(project.get("name", "")).casefold() == name.casefold()]
+    if len(matches) != 1:
+        raise RuntimeError(f"Expected exactly one project named {name!r}, found {len(matches)}")
+    return matches[0]
+
+
+def select_areas(areas: list[dict[str, Any]], fragment: str, all_municipalities: bool) -> list[dict[str, Any]]:
+    if all_municipalities:
+        selected = [area for area in areas if str(area.get("name", "")).casefold().startswith("gemeente ")]
+    else:
+        selected = [area for area in areas if fragment.casefold() in str(area.get("name", "")).casefold()]
+    if not selected:
+        raise RuntimeError("No persisted Area matches the requested scope")
+    selected.sort(key=lambda item: str(item.get("name", "")).casefold())
+    return selected
+
+
+def main() -> int:
+    args = parse_args()
+    base_url = args.base_url.rstrip("/")
+    requested_products = [value.strip() for value in args.products.split(",") if value.strip()]
+    if not requested_products:
+        raise RuntimeError("Select at least one thematic raster product")
+    session = requests.Session()
+    session.headers.update({"User-Agent": "GeoIntel-Thematic-Raster-Operator/1.0"})
+
+    projects = paged_items(session, f"{base_url}/api/v1/projects", args.timeout)
+    project = find_project(projects, args.project_name)
+    project_id = str(project["id"])
+    areas = paged_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", args.timeout)
+    selected_areas = select_areas(areas, args.area, args.all_municipalities)
+    registry = unwrap(session.get(f"{base_url}/api/v1/projects/{project_id}/datasets/thematic-raster/products", timeout=args.timeout))
+    products = {str(item["key"]): item for item in registry.get("items") or []}
+    unknown = sorted(set(requested_products) - set(products))
+    if unknown:
+        raise RuntimeError(f"Products are not present in the canonical registry: {', '.join(unknown)}")
+
+    print(json.dumps({
+        "status": "planned" if args.dry_run else "running",
+        "project_id": project_id,
+        "project_name": project["name"],
+        "area_count": len(selected_areas),
+        "products": requested_products,
+    }, ensure_ascii=False))
+    if args.dry_run:
+        for area in selected_areas:
+            print(json.dumps({"area_id": area["id"], "area_name": area["name"], "bbox": geometry_bbox(area["geometry"])}, ensure_ascii=False))
+        return 0
+
+    results: list[dict[str, Any]] = []
+    for area in selected_areas:
+        bbox = geometry_bbox(area["geometry"])
+        for product_key in requested_products:
+            acquisition = unwrap(session.post(
+                f"{base_url}/api/v1/projects/{project_id}/datasets/thematic-raster/acquire",
+                json={
+                    "bbox": bbox,
+                    "area_id": area["id"],
+                    "product_key": product_key,
+                    "force_refresh": args.force_refresh,
+                },
+                timeout=args.timeout,
+            ))
+            if acquisition.get("status") != "success" or not acquisition.get("output_dataset_id"):
+                raise RuntimeError(f"Acquisition failed for {area['name']} / {product_key}: {acquisition}")
+            dataset_id = str(acquisition["output_dataset_id"])
+            analysis = unwrap(session.post(
+                f"{base_url}/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select",
+                json={"bbox": bbox, "area_id": area["id"]},
+                timeout=args.timeout,
+            ))
+            result = {
+                "area_id": area["id"],
+                "area_name": area["name"],
+                "product_key": product_key,
+                "dataset_id": dataset_id,
+                "reused": bool((acquisition.get("result_json") or {}).get("reused")),
+                "metric": analysis.get("summary"),
+                "coverage_ratio": analysis.get("coverage_ratio"),
+            }
+            results.append(result)
+            print(json.dumps(result, ensure_ascii=False))
+
+    print(json.dumps({"status": "complete", "dataset_count": len(results), "project_id": project_id}, ensure_ascii=False))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

scripts/run_readiness_check.sh

@@ -56,6 +56,8 @@ ${PYTHON_BIN} -m py_compile scripts/provision_buildings_addresses_register.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_dhmv.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_flood_hazards.py
 ${PYTHON_BIN} -m py_compile scripts/provision_regional_flood_hazards.py
+${PYTHON_BIN} -m py_compile scripts/provision_thematic_rasters.py
+${PYTHON_BIN} -m py_compile scripts/provision_mol_soil_map.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