feat: add governed Mol nature value layer

CHANGELOG.md

@@ -7,6 +7,20 @@
 
 # Changelog
 
+## Sprint 204 Governed BWK and Natura 2000 for Mol (2026-07-15)
+
+- Added an explicit operator for the official INBO BWK/Natura 2000 state 2025
+  WFS with complete link/start-index pagination, retained raw page checksums and
+  exact Mol clipping in EPSG:31370.
+- Preserved original BWK evaluation, mapping-unit, origin, habitat, percentage
+  and habitat-origin fields while adding readable, source-faithful properties.
+- Added filtered PostGIS selection metrics for each official BWK value class,
+  Natura 2000 habitat shares, regional biotopes and uncertain habitat gaps.
+- Added Nature Value as a seventh map theme and moved BWK/Natura 2000 from the
+  follow-up catalogue to loaded evidence only after a real Dataset exists.
+- Kept the 2025 edition honest: it is a map state rather than one uniform 2025
+  survey, and PHAB-based partial-selection hectares remain labelled estimates.
+
 ## Sprint 203 Governed hydrology and historical imagery (2026-07-15)
 
 - Added an explicit Waterinfo KiWIS operator for annual water-level and

backend/README.md

@@ -1158,3 +1158,18 @@ observations through the canonical dataset upload API. Every station has its
 own temporal-series key. Water levels and discharges remain Point measurements;
 they are never averaged across stations or presented as municipal water volume.
 Use `--fetch-only` to prepare and audit artifacts without persistence.
+
+## BWK/Natura 2000 state 2025
+
+Run the governed Mol operator after the regional workspace and Mol Area exist:
+
+```bash
+docker exec geointel python /app/scripts/provision_mol_bwk_natura2000.py
+```
+
+The command fetches the official INBO WFS, retains raw checksummed pages,
+clips in EPSG:31370 and imports through DatasetService. `--fetch-only` builds
+evidence without persistence. A conflicting checksum for an already persisted
+state-2025 Mol Dataset fails closed instead of creating a silent replacement.
+PostGIS selection summaries keep BWK value classes separate and label
+PHAB-derived habitat hectares as estimates.

backend/app/services/vector_feature_service.py

@@ -24,6 +24,7 @@ FULL_AREA_CLIPPED_OPERATOR_TOOLS = {
     "provision_regional_grb_buildings.py",
     "provision_regional_grb_context.py",
     "provision_waterinfo_station_history.py",
+    "provision_mol_bwk_natura2000.py",
 }
 
 
@@ -84,6 +85,7 @@ SEMANTIC_SELECTION_METRICS: dict[str, tuple[dict[str, Any], ...]] = {
             "warning": "GRB-percelen zijn een grafische referentie en vormen geen juridische grensopmeting.",
         },
     ),
+    "nature_value": (),
 }
 
 SEMANTIC_COUNT_LABELS = {
@@ -93,6 +95,7 @@ SEMANTIC_COUNT_LABELS = {
     "water": "Waterobjecten",
     "roads": "Wegsegmenten",
     "parcels": "Percelen",
+    "nature_value": "BWK-kaartvlakken",
 }
 
 
@@ -114,6 +117,10 @@ class VectorFeatureService:
             "waterways": "water",
             "road": "roads",
             "parcel": "parcels",
+            "nature": "nature_value",
+            "biodiversity": "nature_value",
+            "bwk": "nature_value",
+            "natura2000": "nature_value",
         }
         for candidate in candidates:
             if not isinstance(candidate, str) or not candidate.strip():
@@ -356,6 +363,17 @@ class VectorFeatureService:
             primary_config = semantic_metrics[0]
 
         metric_configs = [primary_config]
+        configured_metrics = source_metadata.get("selection_metrics")
+        if isinstance(configured_metrics, list):
+            existing_metric_keys = {str(primary_config.get("metric_key") or "")}
+            for configured_item in configured_metrics:
+                if not isinstance(configured_item, dict):
+                    continue
+                metric_key = str(configured_item.get("metric_key") or "").strip()
+                if not metric_key or metric_key in existing_metric_keys:
+                    continue
+                metric_configs.append(dict(configured_item))
+                existing_metric_keys.add(metric_key)
         for semantic_metric in semantic_metrics:
             signature = (semantic_metric["method"], semantic_metric["unit"])
             existing = {
@@ -419,6 +437,20 @@ class VectorFeatureService:
         metric_filter = selection_filter
         if dimension in {1, 2}:
             metric_filter += (func.ST_Dimension(VectorFeature.geometry) == int(dimension),)
+        filter_property = str(config.get("filter_property") or "").strip()
+        filter_values = config.get("filter_values")
+        if filter_property:
+            if not isinstance(filter_values, list) or not filter_values:
+                raise AppError(
+                    code="INVALID_SELECTION_AGGREGATION",
+                    message="Dataset selection metric filter requires one or more values",
+                    details={"dataset_id": str(dataset.id), "filter_property": filter_property},
+                    status_code=500,
+                )
+            normalized_filter_values = [str(value) for value in filter_values]
+            metric_filter += (
+                VectorFeature.properties_json.op("->>")(filter_property).in_(normalized_filter_values),
+            )
 
         if method == "intersection_area":
             measured_geometry = (
@@ -461,7 +493,7 @@ class VectorFeatureService:
             aggregate_function = func.avg if method == "mean" else func.sum
             aggregate_value = (
                 db.query(func.coalesce(aggregate_function(value_expression), 0.0))
-                .filter(*selection_filter)
+                .filter(*metric_filter)
                 .filter(VectorFeature.properties_json.op("->>")(property_name).isnot(None))
                 .scalar()
             )
@@ -469,16 +501,16 @@ class VectorFeatureService:
             if method == "area_weighted_sum" and not full_dataset_area:
                 partial_feature_count = (
                     db.query(func.count(VectorFeature.id))
-                    .filter(*selection_filter)
+                    .filter(*metric_filter)
                     .filter(coverage_ratio < 0.999999)
                     .scalar()
                 )
-                is_estimate = bool(partial_feature_count)
+                is_estimate = bool(config.get("is_estimate", False)) or bool(partial_feature_count)
                 if not is_estimate and config.get("warning_only_when_estimate", True):
                     warning = None
             elif method == "area_weighted_sum":
-                is_estimate = False
+                is_estimate = bool(config.get("is_estimate", False))
-                if config.get("warning_only_when_estimate", True):
+                if not is_estimate and config.get("warning_only_when_estimate", True):
                     warning = None
         elif method != "feature_count":
             raise AppError(

backend/tests/test_sprint204_bwk_natura2000.py

@@ -0,0 +1,270 @@
+from __future__ import annotations
+
+import importlib.util
+import json
+import sys
+from pathlib import Path
+from uuid import uuid4
+
+import pytest
+from shapely.geometry import box, shape
+from shapely.ops import transform as transform_geometry
+
+from app.models import Dataset
+from app.schemas.operations import VectorSelectionSummary
+from app.services.vector_feature_service import VectorFeatureService
+
+
+ROOT = Path(__file__).resolve().parents[2]
+BBOX = {"min_x": 5.0, "min_y": 51.1, "max_x": 5.2, "max_y": 51.3, "crs": "EPSG:4326"}
+
+
+def load_operator():
+    script_path = ROOT / "scripts" / "provision_mol_bwk_natura2000.py"
+    spec = importlib.util.spec_from_file_location("bwk_natura2000_operator", script_path)
+    assert spec is not None
+    assert spec.loader is not None
+    module = importlib.util.module_from_spec(spec)
+    sys.modules[spec.name] = module
+    spec.loader.exec_module(module)
+    return module
+
+
+class FakeResponse:
+    ok = True
+    status_code = 200
+    text = ""
+
+    def __init__(self, payload: dict, url: str):
+        self.payload = payload
+        self.url = url
+        self.content = json.dumps(payload).encode("utf-8")
+
+    def raise_for_status(self):
+        return None
+
+    def json(self):
+        return self.payload
+
+
+class FakeSession:
+    def __init__(self, responses: list[FakeResponse]):
+        self.responses = iter(responses)
+        self.calls: list[tuple[str, dict | None]] = []
+
+    def get(self, url, *, params=None, timeout):  # noqa: ANN001, ARG002
+        self.calls.append((url, params))
+        return next(self.responses)
+
+
+class ScalarQuery:
+    def __init__(self, value: float):
+        self.value = value
+
+    def filter(self, *args):  # noqa: ANN002, ARG002
+        return self
+
+    def scalar(self):
+        return self.value
+
+
+class SequenceScalarSession:
+    def __init__(self, values: list[float]):
+        self.values = iter(values)
+
+    def query(self, *args):  # noqa: ANN002, ARG002
+        return ScalarQuery(next(self.values))
+
+
+def test_official_bwk_contract_and_class_labels_are_fixed() -> None:
+    module = load_operator()
+
+    assert module.WFS_URL == "https://geo.api.vlaanderen.be/BWK/wfs"
+    assert module.TYPE_NAME == "BWK:Bwkhab"
+    assert module.SOURCE_VERSION == "2025"
+    assert module.ATTRIBUTION == "Bron: INBO"
+    assert module.EVALUATION_LABELS == {
+        "z": "Biologisch zeer waardevol",
+        "w": "Biologisch waardevol",
+        "m": "Biologisch minder waardevol",
+        "wz": "Complex van waardevolle en zeer waardevolle elementen",
+        "mwz": "Complex van minder waardevolle, waardevolle en zeer waardevolle elementen",
+        "mz": "Complex van minder waardevolle en zeer waardevolle elementen",
+        "mw": "Complex van minder waardevolle en waardevolle elementen",
+    }
+
+
+def test_wfs_pagination_follows_server_next_links() -> None:
+    module = load_operator()
+    page_one = {
+        "type": "FeatureCollection",
+        "features": [{"id": "one"}],
+        "numberReturned": 1,
+        "links": [{"rel": "next", "href": "https://geo.api.vlaanderen.be/BWK/wfs?STARTINDEX=1"}],
+    }
+    page_two = {"type": "FeatureCollection", "features": [], "numberReturned": 0, "links": []}
+    session = FakeSession(
+        [
+            FakeResponse(page_one, "https://geo.api.vlaanderen.be/BWK/wfs?first"),
+            FakeResponse(page_two, "https://geo.api.vlaanderen.be/BWK/wfs?STARTINDEX=1"),
+        ]
+    )
+
+    pages = list(module.iter_wfs_pages(session, (5.0, 51.1, 5.2, 51.3), page_limit=1000, timeout=30))
+
+    assert len(pages) == 2
+    assert session.calls[0][1]["sortBy"] == "UIDN"
+    assert session.calls[0][1]["srsName"] == "EPSG:4326"
+    assert session.calls[1] == ("https://geo.api.vlaanderen.be/BWK/wfs?STARTINDEX=1", None)
+
+
+def test_wfs_page_limit_without_next_link_uses_controlled_start_index_fallback() -> None:
+    module = load_operator()
+    session = FakeSession(
+        [
+            FakeResponse(
+                {"type": "FeatureCollection", "features": [{"id": "one"}], "numberReturned": 1},
+                "https://geo.api.vlaanderen.be/BWK/wfs",
+            ),
+            FakeResponse(
+                {"type": "FeatureCollection", "features": [], "numberReturned": 0},
+                "https://geo.api.vlaanderen.be/BWK/wfs?startIndex=1",
+            ),
+        ]
+    )
+
+    pages = list(module.iter_wfs_pages(session, (5.0, 51.1, 5.2, 51.3), page_limit=1, timeout=30))
+
+    assert len(pages) == 2
+    assert session.calls[1][1]["startIndex"] == "1"
+
+
+def test_feature_is_clipped_in_lambert72_and_keeps_bwk_habitat_provenance() -> None:
+    module = load_operator()
+    boundary_wgs84 = box(5.10, 51.20, 5.11, 51.21)
+    boundary_lambert72 = transform_geometry(module.TO_LAMBERT72.transform, boundary_wgs84)
+    feature = {
+        "type": "Feature",
+        "id": "Bwkhab.42",
+        "geometry": mapping_box(5.095, 51.195, 5.105, 51.205),
+        "properties": {
+            "UIDN": 42,
+            "EVAL": "wz",
+            "BWKLABEL": "qb + qs",
+            "EENH1": "qb",
+            "EENH2": "qs",
+            "HERK": "225",
+            "HAB1": "9190",
+            "PHAB1": 60,
+            "HAB2": "rbbppm",
+            "PHAB2": 30,
+            "HAB3": "gh",
+            "PHAB3": 10,
+            "HABLEGENDE": "phab",
+            "HERKHAB": "225",
+            "HERKPHAB": "a",
+        },
+    }
+
+    normalized, was_clipped = module.normalize_feature(feature, boundary_lambert72)
+
+    assert normalized is not None
+    assert was_clipped is True
+    normalized_geometry = shape(normalized["geometry"])
+    assert normalized_geometry.difference(boundary_wgs84.buffer(1e-7)).area < 1e-12
+    properties = normalized["properties"]
+    assert properties["bwk_evaluation_code"] == "wz"
+    assert properties["bwk_evaluation_label"].startswith("Complex van waardevolle")
+    assert properties["natura2000_codes"] == "9190"
+    assert properties["regional_biotope_codes"] == "rbbppm"
+    assert properties["natura2000_area_ha"] == pytest.approx(properties["clipped_area_ha"] * 0.6)
+    assert properties["regional_biotope_area_ha"] == pytest.approx(properties["clipped_area_ha"] * 0.3)
+    assert properties["habitat_share_origin_code"] == "a"
+
+
+def mapping_box(min_x: float, min_y: float, max_x: float, max_y: float) -> dict:
+    return {
+        "type": "Polygon",
+        "coordinates": [[
+            [min_x, min_y],
+            [max_x, min_y],
+            [max_x, max_y],
+            [min_x, max_y],
+            [min_x, min_y],
+        ]],
+    }
+
+
+def test_uncertain_habitat_status_is_not_presented_as_confirmed_habitat() -> None:
+    module = load_operator()
+
+    entries, natura_share, regional_share, uncertain_share = module.habitat_breakdown(
+        {"HAB1": "gh", "PHAB1": 100, "HABLEGENDE": "ohab"}
+    )
+
+    assert entries == [{"code": "gh", "share_percent": 100.0}]
+    assert natura_share == 0
+    assert regional_share == 0
+    assert uncertain_share == 100
+
+
+def test_nature_value_summary_returns_separate_official_classes_and_habitat_metrics() -> None:
+    module = load_operator()
+    dataset = Dataset(
+        id=uuid4(),
+        project_id=uuid4(),
+        name="bwk_natura2000_2025_mol.geojson",
+        dataset_type="vector",
+        dataset_role="reference",
+        source_name="inbo_bwk_natura2000",
+        reference_layer_name="nature_value",
+        source_metadata={
+            "theme": "nature_value",
+            "semantic_metrics": False,
+            "selection_aggregation": {
+                "metric_key": "bwk_mapped_area",
+                "method": "intersection_area",
+                "label": "BWK-gekarteerde oppervlakte",
+                "unit": "ha",
+                "geometry_dimension": 2,
+            },
+            "selection_metrics": module.selection_metrics(),
+        },
+    )
+
+    result = VectorFeatureService.summarize_features_by_bbox(
+        SequenceScalarSession([100_000.0, 10_000.0, 20_000.0, 30_000.0, 40_000.0, 5.5, 2.5, 1.5]),
+        dataset=dataset,
+        bbox=BBOX,
+        total_feature_count=125,
+        full_dataset_area=True,
+    )
+
+    assert result["primary_metric_key"] == "bwk_mapped_area"
+    assert result["metric_value"] == 10.0
+    metrics = {item["metric_key"]: item for item in result["metrics"]}
+    assert metrics["bwk_very_valuable_area"]["metric_value"] == 1.0
+    assert metrics["bwk_valuable_area"]["metric_value"] == 2.0
+    assert metrics["bwk_less_valuable_area"]["metric_value"] == 3.0
+    assert metrics["bwk_mixed_value_area"]["metric_value"] == 4.0
+    assert metrics["natura2000_area"]["metric_value"] == 5.5
+    assert metrics["natura2000_area"]["is_estimate"] is True
+    assert metrics["regional_biotope_area"]["metric_value"] == 2.5
+    assert metrics["uncertain_habitat_area"]["metric_value"] == 1.5
+    assert metrics["feature_count"]["metric_value"] == 125
+    VectorSelectionSummary(**result)
+
+
+def test_operator_is_packaged_readiness_checked_and_wired_to_map() -> None:
+    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")
+    service = (ROOT / "backend" / "app" / "services" / "vector_feature_service.py").read_text(encoding="utf-8")
+    map_workspace = (ROOT / "frontend" / "src" / "components" / "map" / "MapWorkspace.tsx").read_text(encoding="utf-8")
+    source_catalog = (ROOT / "frontend" / "src" / "components" / "datasets" / "SourceCatalogPanel.tsx").read_text(encoding="utf-8")
+
+    assert "COPY scripts/provision_mol_bwk_natura2000.py" in dockerfile
+    assert "py_compile scripts/provision_mol_bwk_natura2000.py" in readiness
+    assert '"provision_mol_bwk_natura2000.py"' in service
+    assert "id: 'nature_value'" in map_workspace
+    assert "Natuurwaarde" in map_workspace
+    assert "source.key === 'bwk'" in source_catalog

deploy/unraid/Dockerfile.all-in-one

@@ -78,6 +78,7 @@ COPY scripts/provision_mol_population_history.py /app/scripts/provision_mol_popu
 COPY scripts/provision_mol_historical_landuse.py /app/scripts/provision_mol_historical_landuse.py
 COPY scripts/provision_official_landuse_timeseries.py /app/scripts/provision_official_landuse_timeseries.py
 COPY scripts/provision_waterinfo_station_history.py /app/scripts/provision_waterinfo_station_history.py
+COPY scripts/provision_mol_bwk_natura2000.py /app/scripts/provision_mol_bwk_natura2000.py
 COPY scripts/provision_regional_timeseries.py /app/scripts/provision_regional_timeseries.py
 COPY scripts/geographic_scopes.py /app/scripts/geographic_scopes.py
 COPY scripts/provision_geographic_scope.py /app/scripts/provision_geographic_scope.py

docs/API_CONTRACTS.md

@@ -496,6 +496,11 @@ Rules:
 - Results are generated from persisted PostGIS `vector_features`, not from client-side map data.
 - `feature_count` is the number of GeoJSON features returned in the bounded preview. `total_feature_count` is the exact number of persisted rows intersecting the requested bbox or persisted Area geometry.
 - `summary` keeps one backwards-compatible primary metric and exposes all relevant measurements in `metrics`. Known themes use metric PostGIS calculations: building/forest/water/parcel surfaces in hectares, road and watercourse lengths in kilometres, population in inhabitants and intersecting feature counts as supporting evidence.
+- Governed datasets may declare additional `source_metadata.selection_metrics`.
+  Each metric may constrain one persisted feature-property to an explicit value
+  allowlist before the same PostGIS aggregation runs. The BWK/Natura 2000
+  dataset uses this only for official `EVAL` classes; it does not collapse mixed
+  classes into a made-up score.
 - Area and length calculations transform geometry to Belgian Lambert 72 (`EPSG:31370`); they are never calculated in geographic degrees.
 - Water volume is not inferred from 2D GRB geometry. It remains unavailable until a source provides reliable depth or bathymetry with compatible spatial coverage and provenance.
 - The response is capped by `limit` and returns `truncated=true` when `total_feature_count` exceeds the returned preview.
@@ -1582,6 +1587,9 @@ aggregation method, feature count, estimate status and an optional warning.
 Supported PostGIS aggregations are feature count, intersection area,
 intersection length, numeric sum and area-weighted numeric sum. Area and length
 are measured after transformation to EPSG:31370.
+Configured supplemental metrics can use `filter_property` plus
+`filter_values`. Filters are server-owned dataset metadata, not arbitrary
+client SQL or request expressions.
 
 ### PATCH `/api/v1/projects/{project_id}/datasets/{dataset_id}/temporal`
 

docs/DATA_SOURCES.md

@@ -227,6 +227,27 @@ until a governed operator import, provenance record and validation pass exist.
   LiDAR, suitable for elevation, slope and drainage. It does not provide water
   depth.
 
+## Governed BWK and Natura 2000 state 2025
+
+`scripts/provision_mol_bwk_natura2000.py` uses the official production WFS
+`https://geo.api.vlaanderen.be/BWK/wfs`, feature type `BWK:Bwkhab`. It follows
+server `next` links and uses WFS `startIndex` only when a full page omits that
+link. A short final page is required for completeness. Raw response pages,
+request URLs and SHA256 checksums are retained.
+
+The source bbox is reduced to the exact persisted Mol Area in Belgian Lambert
+72. Only valid polygonal intersections are emitted, transformed back to
+EPSG:4326 and imported through DatasetService with `source_name` equal to
+`inbo_bwk_natura2000` and `reference_layer_name` equal to `nature_value`.
+Required attribution is `Bron: INBO`.
+
+`EVAL` classes are reported separately. Natura 2000 and regionally important
+biotope hectares use official `PHAB` shares; uncertain `ohab` polygons remain
+knowledge gaps. The official report warns that parts of the 2025 edition,
+especially outside Habitats Directive areas, can still be based on older
+mapping. GeoIntel therefore preserves origin fields and never presents this
+single edition as a historical trend.
+
 ## Governed Waterinfo station history
 
 `scripts/provision_waterinfo_station_history.py` uses the public Waterinfo

docs/DATA_SPECIFICATION.md

@@ -213,6 +213,21 @@ series may not be merged into one area-wide value. Point water level and
 discharge may not be converted to water volume without governed compatible
 depth/profile data.
 
+### BWK and Natura 2000 polygons
+
+The governed state-2025 import uses `BWK:Bwkhab` polygons. Geometry is fetched
+in EPSG:4326, validated and clipped against the persisted Mol boundary in
+EPSG:31370, then persisted in EPSG:4326. Required retained fields include
+`EVAL`, `EENH1..8`, `V1..3`, `HERK`, `BWKLABEL`, `HAB1..5`, `PHAB1..5`,
+`HERKHAB`, `HERKPHAB` and `HABLEGENDE`.
+
+BWK value classes remain independent. Natura 2000 codes, regionally important
+biotope codes and uncertain `ohab` knowledge gaps are also independent. PHAB
+shares can be theoretical source allocations; derived hectares therefore keep
+an estimate warning, especially for partial rectangle intersections. State
+2025 is a product edition and must not be presented as a uniform 2025 field
+survey or an annual time series.
+
 ### Vector
 
 Ondersteund:

docs/STORAGE_ARCHITECTURE.md

@@ -109,6 +109,13 @@ product/layer, request/spatial hash, temporal validity and limitations are held
 in source/provenance metadata. Browser PNG rendering is derived on request and
 does not replace the stored GeoTIFF.
 
+BWK/Natura 2000 evidence lives under
+`storage/operator-evidence/bwk-natura2000-2025/mol/`. The `raw/` directory
+contains immutable WFS pages; the adjacent manifest records their URLs,
+checksums, boundary checksum, page/feature counts, clipping diagnostics and
+official class totals. The final GeoJSON is the retained source artifact while
+`datasets` and `vector_features` remain the canonical queryable PostGIS state.
+
 Offline demo export artifacts can be inspected and cleaned with:
 
 ```bash

frontend/README.md

@@ -470,6 +470,12 @@ The workbench shell clamps page-level horizontal overflow on mobile while keepin
 
 The QA/QC result list includes client-side search, status and check-type filters plus a latest-results cap so long-lived demo projects remain scan-friendly without changing the API response shape.
 
+The current-state Map explorer includes `Natuurwaarde` when the governed INBO
+BWK/Natura 2000 state-2025 Dataset is loaded. It shows exact mapped BWK area,
+separate official value classes, PHAB-derived Natura 2000 and regional-biotope
+hectares, uncertain habitat knowledge gaps and source provenance. This single
+edition is not exposed as a fabricated annual series in Evolution.
+
 ## Useful repository scripts
 
 - `bash scripts/frontend_install.sh`

frontend/src/components/datasets/SourceCatalogPanel.tsx

@@ -8,6 +8,7 @@ const THEME_LABELS: Record<string, string> = {
   buildings: 'Bebouwing',
   population: 'Bevolking',
   forest: 'Bos',
+  nature_value: 'Natuurwaarde',
   water: 'Water',
   roads: 'Wegen en transport',
   parcels: 'Percelen',
@@ -103,9 +104,11 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E
       dataset.source_name === 'digitaal_vlaanderen_orthophoto' &&
       String(dataset.source_metadata?.['product_key'] ?? 'most_recent') !== 'most_recent',
   )
+  const bwkDatasets = ready.filter((dataset) => dataset.source_name === 'inbo_bwk_natura2000')
   const pendingSources = AVAILABLE_SOURCES.filter((source) => {
     if (source.key === 'waterinfo') return waterinfoDatasets.length === 0
     if (source.key === 'historical_orthophoto') return historicalOrthophotos.length === 0
+    if (source.key === 'bwk') return bwkDatasets.length === 0
     return true
   })
   const themes = Object.keys(THEME_LABELS).map((theme) => {
@@ -165,7 +168,7 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E
         ))}
       </div>
 
-      {waterinfoDatasets.length > 0 || historicalOrthophotos.length > 0 ? (
+      {waterinfoDatasets.length > 0 || historicalOrthophotos.length > 0 || bwkDatasets.length > 0 ? (
         <div className="source-catalog-loaded" aria-label="Aanvullende ingeladen bronnen">
           {waterinfoDatasets.length > 0 ? (
             <article>
@@ -181,6 +184,13 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E
               <p>Officiële jaargangen en periodes zijn via de kaart beschikbaar zonder actuele GRB-validatie.</p>
             </article>
           ) : null}
+          {bwkDatasets.length > 0 ? (
+            <article>
+              <strong>BWK / Natura 2000</strong>
+              <span>Toestand 2025 · {bwkDatasets.reduce((total, dataset) => total + (dataset.feature_count ?? 0), 0).toLocaleString('nl-BE')} kaartvlakken</span>
+              <p>Biologische waardering en habitataandelen blijven afzonderlijke, brongetrouwe metingen met herkomstinformatie.</p>
+            </article>
+          ) : null}
         </div>
       ) : null}
 

frontend/src/components/map/MapWorkspace.tsx

@@ -13,7 +13,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' | 'water' | 'roads' | 'parcels'
+type DataThemeId = 'buildings' | 'population' | 'forest' | 'nature_value' | 'water' | 'roads' | 'parcels'
 
 interface DataTheme {
   id: DataThemeId
@@ -51,6 +51,13 @@ const DATA_THEMES: DataTheme[] = [
     description: 'Bos, natuur en groenbedekking uit een ingeladen vectorbron.',
     tokens: ['forest', 'forestry', 'woodland', 'bos', 'groen', 'vegetation'],
   },
+  {
+    id: 'nature_value',
+    label: 'Natuurwaarde',
+    shortLabel: 'BWK-oppervlakte',
+    description: 'Biologische waardering, Natura 2000-habitat en regionaal belangrijke biotopen uit de BWK.',
+    tokens: ['nature_value', 'nature value', 'natuurwaarde', 'bwk', 'natura2000', 'natura 2000', 'biodiversity'],
+  },
   {
     id: 'water',
     label: 'Water',
@@ -78,6 +85,7 @@ const DATA_THEME_MAP_STYLES: Record<DataThemeId, { fill: string; line: string }>
   buildings: { fill: '#d45f3d', line: '#9f3e24' },
   population: { fill: '#7559a6', line: '#5b3f88' },
   forest: { fill: '#347950', line: '#225f3b' },
+  nature_value: { fill: '#9a4f64', line: '#74364a' },
   water: { fill: '#2676a8', line: '#155b85' },
   roads: { fill: '#6b7280', line: '#4b5563' },
   parcels: { fill: '#a7792f', line: '#7d571f' },
@@ -111,6 +119,7 @@ function pickThemeDataset(datasets: DatasetCreateResponse[], theme: DataTheme):
       (dataset.reference_layer_name && theme.tokens.includes(dataset.reference_layer_name.toLowerCase()) ? 1_000_000 : 0) +
       (dataset.source_name === 'grb' ? 100_000 : 0) +
       (dataset.source_name === 'department_omgeving_land_use' ? 90_000 : 0) +
+      (dataset.source_name === 'inbo_bwk_natura2000' ? 95_000 : 0) +
       (dataset.dataset_role === 'reference' ? 10_000 : 0) +
       (dataset.observed_at ? new Date(dataset.observed_at).getTime() / 100_000_000 : 0) +
       (dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0)

frontend/src/styles/app.css

@@ -5753,6 +5753,7 @@ section {
 .geo-theme-symbol-buildings { background: #d45f3d; }
 .geo-theme-symbol-population { background: #7559a6; }
 .geo-theme-symbol-forest { background: #347950; }
+.geo-theme-symbol-nature_value { background: #9a4f64; }
 .geo-theme-symbol-water { background: #2676a8; }
 .geo-theme-symbol-roads { background: #6b7280; }
 .geo-theme-symbol-parcels { background: #a7792f; }

scripts/README.md

@@ -1438,6 +1438,25 @@ missing discharge series is reported without synthesizing values. Different
 stations remain separate temporal series and may not be treated as area-wide
 water level or volume.
 
+## Governed BWK/Natura 2000 state 2025 for Mol
+
+Build and validate raw WFS evidence without importing:
+
+```bash
+docker exec geointel python /app/scripts/provision_mol_bwk_natura2000.py --fetch-only
+```
+
+Import after reviewing the manifest:
+
+```bash
+docker exec geointel python /app/scripts/provision_mol_bwk_natura2000.py
+```
+
+The operator uses the official `BWK:Bwkhab` layer, follows complete WFS
+pagination, retains raw JSON/checksums, clips polygon geometry in EPSG:31370
+and writes only through the Dataset upload flow. `--force` refetches source
+evidence but never silently replaces a conflicting persisted 2025 artifact.
+
 ## Tower deployment
 
 Push the local branch to Gitea, then rebuild the Unraid/Tower Docker runtime:

scripts/provision_mol_bwk_natura2000.py

@@ -0,0 +1,815 @@
+"""Provision the official 2025 BWK/Natura 2000 dataset for Mol.
+
+The operator follows every WFS page, retains raw checksummed responses, clips
+polygon geometry to the persisted Mol Area in EPSG:31370 and imports the final
+GeoJSON through the canonical DatasetService upload route. It never writes
+directly to vector_features.
+"""
+
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import os
+import re
+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://geo.api.vlaanderen.be/BWK/wfs"
+CATALOG_URL = (
+    "https://www.vlaanderen.be/datavindplaats/catalogus/"
+    "biologische-waarderingskaart-en-natura-2000-habitatkaart-toestand-2025"
+)
+REPORT_URL = "https://doi.org/10.21436/inbor.129502912"
+TYPE_NAME = "BWK:Bwkhab"
+SOURCE_VERSION = "2025"
+PUBLICATION_DATE = "2025-12-10T00:00:00Z"
+ATTRIBUTION = "Bron: INBO"
+DEFAULT_API_URL = "http://127.0.0.1:8000"
+DEFAULT_OUTPUT_DIR = "/app/storage/operator-evidence/bwk-natura2000-2025/mol"
+DEFAULT_PROJECT_NAME = "Kempen Regional Workbench"
+DEFAULT_AREA_NAME = "Gemeente Mol - officiele grens"
+DATASET_FILENAME = "bwk_natura2000_2025_mol.geojson"
+MANIFEST_FILENAME = "bwk_natura2000_2025_mol.manifest.json"
+SCHEMA_VERSION = 1
+
+EVALUATION_LABELS = {
+    "z": "Biologisch zeer waardevol",
+    "w": "Biologisch waardevol",
+    "m": "Biologisch minder waardevol",
+    "wz": "Complex van waardevolle en zeer waardevolle elementen",
+    "mwz": "Complex van minder waardevolle, waardevolle en zeer waardevolle elementen",
+    "mz": "Complex van minder waardevolle en zeer waardevolle elementen",
+    "mw": "Complex van minder waardevolle en waardevolle elementen",
+}
+
+HABITAT_STATUS_LABELS = {
+    "gh": "Geen Natura 2000-habitat aangeduid",
+    "hab": "Volledig habitatwaardig kaartvlak",
+    "phab": "Gedeeltelijk habitatwaardig kaartvlak",
+    "ohab": "Onzeker habitat of kennislacune",
+}
+
+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 official BWK/Natura 2000 state 2025 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-name", default=DEFAULT_AREA_NAME)
+    parser.add_argument(
+        "--output-dir",
+        type=Path,
+        default=Path(os.environ.get("GEOINTEL_BWK_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)),
+    )
+    parser.add_argument("--page-limit", type=int, default=1000)
+    parser.add_argument("--max-features", type=int, default=30_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", help="Refetch source pages and rebuild the artifact.")
+    parser.add_argument("--fetch-only", action="store_true", help="Build evidence without importing a Dataset.")
+    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-BWK-Natura2000-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 reading the 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)
+            return
+        if isinstance(candidate, MultiPolygon):
+            polygons.extend(part for part in candidate.geoms if not part.is_empty)
+            return
+        if 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_name: 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)
+    normalized_target = area_name.casefold().replace("officiele", "officiële")
+    area = next(
+        (
+            item
+            for item in areas
+            if str(item.get("name") or "").casefold().replace("officiele", "officiële") == normalized_target
+        ),
+        None,
+    )
+    if not area:
+        area = next((item for item in areas if "gemeente mol" in str(item.get("name") or "").casefold()), None)
+    if not area or not area.get("geometry"):
+        raise RuntimeError(f"Persisted Mol Area {area_name!r} with geometry is missing")
+    boundary_wgs84 = polygonal_geometry(shape(area["geometry"]))
+    if boundary_wgs84 is None:
+        raise RuntimeError("Persisted Mol Area geometry is invalid or not polygonal")
+    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 next_page_url(payload: dict[str, Any]) -> str | None:
+    for link in payload.get("links") or []:
+        if str(link.get("rel") or "").lower() == "next" and link.get("href"):
+            return str(link["href"])
+    return None
+
+
+def iter_wfs_pages(
+    session: requests.Session,
+    bbox: tuple[float, float, float, float],
+    *,
+    page_limit: int,
+    timeout: int,
+) -> Iterable[tuple[dict[str, Any], str, bytes]]:
+    url: str | None = WFS_URL
+    base_params: dict[str, str] = {
+        "service": "WFS",
+        "version": "2.0.0",
+        "request": "GetFeature",
+        "typeNames": TYPE_NAME,
+        "srsName": "EPSG:4326",
+        "bbox": ",".join(f"{value:.8f}" for value in bbox) + ",EPSG:4326",
+        "count": str(page_limit),
+        "sortBy": "UIDN",
+        "outputFormat": "application/json",
+    }
+    params: dict[str, str] | None = dict(base_params)
+    seen: set[str] = set()
+    next_start_index = 0
+    while url:
+        request_key = requests.Request("GET", url, params=params).prepare().url or url
+        if request_key in seen:
+            raise RuntimeError("BWK WFS pagination loop detected")
+        seen.add(request_key)
+        response = session.get(url, params=params, timeout=timeout)
+        params = None
+        response.raise_for_status()
+        payload = response.json()
+        if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection":
+            raise RuntimeError("BWK WFS returned an invalid FeatureCollection")
+        yield payload, response.url, response.content
+        next_url = next_page_url(payload)
+        number_returned = int(payload.get("numberReturned") or len(payload.get("features") or []))
+        next_start_index += number_returned
+        if next_url:
+            url = next_url
+            params = None
+        elif number_returned >= page_limit:
+            url = WFS_URL
+            params = {**base_params, "startIndex": str(next_start_index)}
+        else:
+            url = None
+            params = None
+
+
+def habitat_breakdown(properties: dict[str, Any]) -> tuple[list[dict[str, Any]], float, float, float]:
+    entries: list[dict[str, Any]] = []
+    natura_share = 0.0
+    regional_share = 0.0
+    uncertain_share = 0.0
+    for index in range(1, 6):
+        code = str(properties.get(f"HAB{index}") or "").strip()
+        try:
+            share = max(0.0, min(100.0, float(properties.get(f"PHAB{index}") or 0.0)))
+        except (TypeError, ValueError):
+            share = 0.0
+        if not code or share <= 0:
+            continue
+        normalized = code.lower()
+        entries.append({"code": code, "share_percent": share})
+        if re.match(r"^\d", normalized):
+            natura_share += share
+        elif normalized.startswith("rbb"):
+            regional_share += share
+        elif normalized == "ohab":
+            uncertain_share += share
+    status = str(properties.get("HABLEGENDE") or "").strip().lower()
+    if status == "ohab" and uncertain_share <= 0:
+        uncertain_share = 100.0
+    return entries, min(natura_share, 100.0), min(regional_share, 100.0), min(uncertain_share, 100.0)
+
+
+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_properties = dict(feature.get("properties") or {})
+    source_id = str(feature.get("id") or raw_properties.get("UIDN") or raw_properties.get("OIDN") or "").strip()
+    if not source_id:
+        source_id = sha256_bytes(json.dumps(geometry_payload, sort_keys=True).encode("utf-8"))
+    stable_id = f"BWK:Bwkhab:{raw_properties.get('UIDN') or source_id}"
+    evaluation_code = str(raw_properties.get("EVAL") or "").strip().lower()
+    habitat_status = str(raw_properties.get("HABLEGENDE") or "").strip().lower()
+    habitats, natura_share, regional_share, uncertain_share = habitat_breakdown(raw_properties)
+    clipped_area_ha = float(clipped_lambert72.area) / 10_000.0
+    natura_codes = [entry["code"] for entry in habitats if re.match(r"^\d", str(entry["code"]))]
+    regional_codes = [entry["code"] for entry in habitats if str(entry["code"]).lower().startswith("rbb")]
+
+    properties = {
+        **raw_properties,
+        "source_name": "inbo_bwk_natura2000",
+        "source_collection": TYPE_NAME,
+        "source_feature_id": stable_id,
+        "reference_layer_name": "nature_value",
+        "theme": "nature_value",
+        "authority_level": "authoritative",
+        "coverage_scope": "municipality",
+        "municipality": "Mol",
+        "nis_code": "13025",
+        "source_version": SOURCE_VERSION,
+        "attribution": ATTRIBUTION,
+        "bwk_evaluation_code": evaluation_code or "unknown",
+        "bwk_evaluation_label": EVALUATION_LABELS.get(evaluation_code, "Onbekende of ontbrekende BWK-waardering"),
+        "bwk_label": str(raw_properties.get("BWKLABEL") or "").strip(),
+        "bwk_units": ", ".join(
+            str(raw_properties.get(f"EENH{index}") or "").strip()
+            for index in range(1, 9)
+            if str(raw_properties.get(f"EENH{index}") or "").strip()
+        ),
+        "bwk_origin_code": str(raw_properties.get("HERK") or "").strip(),
+        "habitat_status_code": habitat_status or "unknown",
+        "habitat_status_label": HABITAT_STATUS_LABELS.get(habitat_status, "Onbekende of ontbrekende habitatstatus"),
+        "natura2000_codes": ", ".join(natura_codes),
+        "regional_biotope_codes": ", ".join(regional_codes),
+        "habitat_entries": habitats,
+        "habitat_origin_code": str(raw_properties.get("HERKHAB") or "").strip(),
+        "habitat_share_origin_code": str(raw_properties.get("HERKPHAB") or "").strip(),
+        "clipped_area_ha": round(clipped_area_ha, 8),
+        "natura2000_share_percent": natura_share,
+        "regional_biotope_share_percent": regional_share,
+        "uncertain_habitat_share_percent": uncertain_share,
+        "natura2000_area_ha": round(clipped_area_ha * natura_share / 100.0, 8),
+        "regional_biotope_area_ha": round(clipped_area_ha * regional_share / 100.0, 8),
+        "uncertain_habitat_area_ha": round(clipped_area_ha * uncertain_share / 100.0, 8),
+        "habitat_share_spatial_limitation": (
+            "PHAB percentages apply to the full source polygon; a partial map selection scales them by intersected area."
+        ),
+    }
+    return {
+        "type": "Feature",
+        "id": stable_id,
+        "geometry": mapping(clipped_wgs84),
+        "properties": properties,
+    }, was_clipped
+
+
+def reusable_artifact(output_dir: Path) -> tuple[Path, Path, dict[str, Any]] | None:
+    artifact_path = output_dir / DATASET_FILENAME
+    manifest_path = output_dir / MANIFEST_FILENAME
+    if not artifact_path.is_file() or not manifest_path.is_file():
+        return None
+    try:
+        manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
+    except (OSError, ValueError):
+        return None
+    if manifest.get("schema_version") != SCHEMA_VERSION or manifest.get("source_version") != SOURCE_VERSION:
+        return None
+    if manifest.get("artifact_sha256") != sha256_file(artifact_path):
+        return None
+    for page in manifest.get("raw_pages") or []:
+        page_path = output_dir / str(page.get("path") or "")
+        if not page_path.is_file() or page.get("sha256") != sha256_file(page_path):
+            return None
+    return artifact_path, manifest_path, manifest
+
+
+def prepare_artifact(
+    session: requests.Session,
+    boundary_wgs84,
+    output_dir: Path,
+    *,
+    page_limit: int,
+    max_features: int,
+    timeout: int,
+    force: bool,
+) -> tuple[Path, Path, dict[str, Any]]:
+    if not force:
+        reusable = reusable_artifact(output_dir)
+        if reusable:
+            return reusable
+
+    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
+
+    for page_number, (payload, source_url, raw_bytes) in enumerate(
+        iter_wfs_pages(session, boundary_wgs84.bounds, page_limit=page_limit, timeout=timeout),
+        start=1,
+    ):
+        page_path = raw_dir / f"bwk_bwkhab_page_{page_number:05d}.json"
+        write_bytes_atomic(page_path, raw_bytes)
+        page_features = list(payload.get("features") or [])
+        raw_feature_count += len(page_features)
+        if raw_feature_count > max_features:
+            raise RuntimeError(
+                f"BWK WFS returned more than 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(page_features),
+                "source_url": source_url,
+            }
+        )
+        source_urls.append(source_url)
+        for source_feature in page_features:
+            source_identity = str(
+                source_feature.get("id")
+                or (source_feature.get("properties") or {}).get("UIDN")
+                or ""
+            )
+            if source_identity and source_identity in seen_ids:
+                duplicate_count += 1
+                continue
+            if source_identity:
+                seen_ids.add(source_identity)
+            normalized, was_clipped = normalize_feature(source_feature, boundary_lambert72)
+            if normalized is None:
+                rejected_count += 1
+                continue
+            if was_clipped:
+                clipped_count += 1
+            retained.append(normalized)
+
+    if not retained:
+        raise RuntimeError("BWK WFS returned no valid polygon features inside the persisted Mol Area")
+
+    evaluation_area: dict[str, float] = defaultdict(float)
+    habitat_status_area: dict[str, float] = defaultdict(float)
+    natura_area = 0.0
+    regional_area = 0.0
+    uncertain_area = 0.0
+    for feature in retained:
+        properties = feature["properties"]
+        area = float(properties["clipped_area_ha"])
+        evaluation_area[str(properties["bwk_evaluation_code"])] += area
+        habitat_status_area[str(properties["habitat_status_code"])] += area
+        natura_area += float(properties["natura2000_area_ha"])
+        regional_area += float(properties["regional_biotope_area_ha"])
+        uncertain_area += float(properties["uncertain_habitat_area_ha"])
+
+    generated_at = utc_now()
+    artifact = {
+        "type": "FeatureCollection",
+        "name": "BWK en Natura 2000 - toestand 2025 - Gemeente Mol",
+        "features": retained,
+        "source": "INBO BWK/Natura 2000 WFS",
+        "source_version": SOURCE_VERSION,
+        "attribution": ATTRIBUTION,
+        "catalog_url": CATALOG_URL,
+        "report_url": REPORT_URL,
+        "coverage_scope": "municipality",
+        "municipality": "Mol",
+        "nis_code": "13025",
+        "reference_truncated": False,
+        "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,
+        "source_type_name": TYPE_NAME,
+        "wfs_url": WFS_URL,
+        "catalog_url": CATALOG_URL,
+        "report_url": REPORT_URL,
+        "attribution": ATTRIBUTION,
+        "generated_at": generated_at,
+        "crs_source": "EPSG:4326",
+        "crs_clip": "EPSG:31370",
+        "crs_persisted": "EPSG:4326",
+        "boundary_sha256": sha256_bytes(json.dumps(mapping(boundary_wgs84), sort_keys=True).encode("utf-8")),
+        "boundary_bbox_wgs84": list(boundary_wgs84.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,
+        "evaluation_area_ha": {key: round(value, 6) for key, value in sorted(evaluation_area.items())},
+        "habitat_status_area_ha": {key: round(value, 6) for key, value in sorted(habitat_status_area.items())},
+        "natura2000_area_ha": round(natura_area, 6),
+        "regional_biotope_area_ha": round(regional_area, 6),
+        "uncertain_habitat_area_ha": round(uncertain_area, 6),
+        "artifact_path": str(artifact_path),
+        "artifact_sha256": sha256_file(artifact_path),
+        "artifact_size_bytes": artifact_path.stat().st_size,
+        "limitations": [
+            "The 2025 edition is the best available map state, not one uniform 2025 field survey.",
+            "PHAB percentages may be theoretical shares; partial selections scale shares by intersected polygon area.",
+            "The real field situation remains authoritative for policy and legal use.",
+        ],
+    }
+    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]]:
+    share_warning = (
+        "Oppervlakte op basis van de officiele PHAB-aandelen. Automatisch verdeelde aandelen kunnen lokaal "
+        "afwijken van de terreinsituatie."
+    )
+    return [
+        {
+            "metric_key": "bwk_very_valuable_area",
+            "method": "intersection_area",
+            "label": "Biologisch zeer waardevol",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "bwk_evaluation_code",
+            "filter_values": ["z"],
+        },
+        {
+            "metric_key": "bwk_valuable_area",
+            "method": "intersection_area",
+            "label": "Biologisch waardevol",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "bwk_evaluation_code",
+            "filter_values": ["w"],
+        },
+        {
+            "metric_key": "bwk_less_valuable_area",
+            "method": "intersection_area",
+            "label": "Biologisch minder waardevol",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "bwk_evaluation_code",
+            "filter_values": ["m"],
+        },
+        {
+            "metric_key": "bwk_mixed_value_area",
+            "method": "intersection_area",
+            "label": "Gemengde BWK-waardering",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "filter_property": "bwk_evaluation_code",
+            "filter_values": ["wz", "mwz", "mz", "mw"],
+        },
+        {
+            "metric_key": "natura2000_area",
+            "method": "area_weighted_sum",
+            "property": "natura2000_area_ha",
+            "label": "Natura 2000-habitat",
+            "unit": "ha",
+            "is_estimate": True,
+            "warning": share_warning,
+            "warning_only_when_estimate": False,
+        },
+        {
+            "metric_key": "regional_biotope_area",
+            "method": "area_weighted_sum",
+            "property": "regional_biotope_area_ha",
+            "label": "Regionaal belangrijk biotoop",
+            "unit": "ha",
+            "is_estimate": True,
+            "warning": share_warning,
+            "warning_only_when_estimate": False,
+        },
+        {
+            "metric_key": "uncertain_habitat_area",
+            "method": "area_weighted_sum",
+            "property": "uncertain_habitat_area_ha",
+            "label": "Onzeker habitat / kennislacune",
+            "unit": "ha",
+            "is_estimate": True,
+            "warning": "Dit is een maximale potentiele oppervlakte van kaartvlakken met de status ohab, geen bevestigd habitat.",
+            "warning_only_when_estimate": False,
+        },
+    ]
+
+
+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]:
+    source_metadata = {
+        "provider": "Instituut voor Natuur- en Bosonderzoek",
+        "theme": "nature_value",
+        "layer_name": "BWK/Natura 2000 toestand 2025",
+        "authority_level": "authoritative",
+        "coverage_scope": "municipality",
+        "municipality": "Mol",
+        "nis_code": "13025",
+        "feature_count": manifest["feature_count"],
+        "geometry_clipped_to_area": True,
+        "semantic_metrics": False,
+        "attribution": ATTRIBUTION,
+        "catalog_url": CATALOG_URL,
+        "report_url": REPORT_URL,
+        "selection_aggregation": {
+            "metric_key": "bwk_mapped_area",
+            "method": "intersection_area",
+            "label": "BWK-gekarteerde oppervlakte",
+            "unit": "ha",
+            "geometry_dimension": 2,
+            "warning": (
+                "Uitgave 2025 is de best beschikbare kaarttoestand, maar niet elk kaartvlak is in 2025 op terrein gekarteerd."
+            ),
+        },
+        "selection_metrics": selection_metrics(),
+    }
+    provenance_metadata = {
+        "operator_tool": "provision_mol_bwk_natura2000.py",
+        "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
+        "source_type_name": TYPE_NAME,
+        "wfs_url": WFS_URL,
+        "catalog_url": CATALOG_URL,
+        "report_url": REPORT_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": "inbo_bwk_natura2000",
+                "reference_layer_name": "nature_value",
+                "source_metadata_json": json.dumps(source_metadata, ensure_ascii=False),
+                "provenance_metadata_json": json.dumps(provenance_metadata, ensure_ascii=False),
+                "area_id": area_id,
+                "observed_at": PUBLICATION_DATE,
+                "temporal_granularity": "snapshot",
+                "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 > 5000 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("/")
+    try:
+        with requests.Session() as api_session:
+            project_id, area_id, boundary, existing = locate_workspace(
+                api_session,
+                base_url,
+                args.project_name,
+                args.area_name,
+                args.import_timeout,
+            )
+            with source_session() as official_session:
+                artifact_path, manifest_path, manifest = prepare_artifact(
+                    official_session,
+                    boundary,
+                    args.output_dir,
+                    page_limit=args.page_limit,
+                    max_features=args.max_features,
+                    timeout=args.request_timeout,
+                    force=args.force,
+                )
+            existing_dataset = next(
+                (
+                    item
+                    for item in existing
+                    if item.get("source_name") == "inbo_bwk_natura2000"
+                    and item.get("source_version") == SOURCE_VERSION
+                    and str(item.get("area_id") or "") == area_id
+                ),
+                None,
+            )
+            if existing_dataset:
+                persisted_checksum = str(existing_dataset.get("checksum_sha256") or "")
+                if persisted_checksum and persisted_checksum != manifest["artifact_sha256"]:
+                    raise RuntimeError(
+                        "A different BWK 2025 Mol artifact is already persisted; refusing a silent duplicate or replacement"
+                    )
+                result = {
+                    "status": "existing",
+                    "dataset_id": existing_dataset["id"],
+                    "feature_count": existing_dataset.get("feature_count") or manifest["feature_count"],
+                }
+            elif args.fetch_only:
+                result = {
+                    "status": "prepared",
+                    "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",
+                    "dataset_id": dataset["id"],
+                    "feature_count": dataset.get("feature_count"),
+                }
+    except (KeyError, OSError, RuntimeError, ValueError, requests.RequestException) as exc:
+        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
+        return 1
+
+    print(
+        json.dumps(
+            {
+                "status": "ok",
+                "project": args.project_name,
+                "area": args.area_name,
+                "source_version": SOURCE_VERSION,
+                "manifest_path": str(manifest_path),
+                "metrics": {
+                    "evaluation_area_ha": manifest["evaluation_area_ha"],
+                    "natura2000_area_ha": manifest["natura2000_area_ha"],
+                    "regional_biotope_area_ha": manifest["regional_biotope_area_ha"],
+                    "uncertain_habitat_area_ha": manifest["uncertain_habitat_area_ha"],
+                },
+                "result": result,
+            },
+            ensure_ascii=False,
+            indent=2,
+        )
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

scripts/run_readiness_check.sh

@@ -48,6 +48,7 @@ ${PYTHON_BIN} -m py_compile scripts/provision_mol_population_history.py
 ${PYTHON_BIN} -m py_compile scripts/provision_mol_historical_landuse.py
 ${PYTHON_BIN} -m py_compile scripts/provision_official_landuse_timeseries.py
 ${PYTHON_BIN} -m py_compile scripts/provision_waterinfo_station_history.py
+${PYTHON_BIN} -m py_compile scripts/provision_mol_bwk_natura2000.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