Add regional DOV soil operator

backend/app/services/vector_feature_service.py

@@ -30,6 +30,7 @@ FULL_AREA_CLIPPED_OPERATOR_TOOLS = {
     "provision_agricultural_parcel_history.py",
     "provision_buildings_addresses_register.py",
     "provision_mol_soil_map.py",
+    "provision_regional_soil_map.py",
 }
 
 SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS = {
@@ -40,6 +41,7 @@ SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS = {
 
 PRECLIPPED_MUNICIPALITY_PARTITION_OPERATOR_TOOLS = {
     "provision_regional_bwk_natura2000.py",
+    "provision_regional_soil_map.py",
 }
 
 

backend/tests/test_sprint217_regional_dov_soil_map.py

@@ -0,0 +1,287 @@
+from __future__ import annotations
+
+import gzip
+import importlib.util
+import json
+from pathlib import Path
+import sys
+
+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]
+SCRIPTS = ROOT / "scripts"
+if str(SCRIPTS) not in sys.path:
+    sys.path.insert(0, str(SCRIPTS))
+
+
+def load_script():
+    path = SCRIPTS / "provision_regional_soil_map.py"
+    spec = importlib.util.spec_from_file_location("test_provision_regional_soil_map", path)
+    assert spec is not None and 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:
+    status_code = 200
+    ok = True
+    text = ""
+
+    def __init__(self, payload):
+        self.payload = payload
+
+    def json(self):
+        return self.payload
+
+
+class FakeApiSession:
+    def __init__(self, payload):
+        self.payload = payload
+        self.calls = []
+
+    def post(self, url, **kwargs):
+        self.calls.append((url, kwargs))
+        return FakeResponse({"data": self.payload})
+
+
+def source_feature(feature_id: str, geometry) -> dict:
+    return {
+        "type": "Feature",
+        "id": feature_id,
+        "geometry": mapping(geometry),
+        "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_normalized_regional_features_keep_member_identity_and_unique_ids() -> None:
+    module = load_script()
+    boundary = box(5.0, 51.0, 5.1, 51.1)
+    boundary_lambert72 = transform_geometry(module.soil.TO_LAMBERT72.transform, boundary)
+    feature = source_feature("bodemtypes.1", box(4.98, 51.02, 5.05, 51.08))
+
+    left, _ = module.soil.normalize_feature(
+        feature,
+        boundary_lambert72,
+        municipality="Mol",
+        nis_code="13025",
+        coverage_scope="test-region",
+        feature_id_suffix="13025",
+    )
+    right, _ = module.soil.normalize_feature(
+        feature,
+        boundary_lambert72,
+        municipality="Balen",
+        nis_code="13003",
+        coverage_scope="test-region",
+        feature_id_suffix="13003",
+    )
+
+    assert left is not None and right is not None
+    assert left["id"] == "bodemtypes.1:13025"
+    assert right["id"] == "bodemtypes.1:13003"
+    assert left["properties"]["municipality"] == "Mol"
+    assert left["properties"]["coverage_scope"] == "test-region"
+    assert shape(left["geometry"]).within(boundary.buffer(1e-7))
+
+
+def test_partition_retains_gzipped_source_and_is_checksum_reusable(tmp_path: Path, monkeypatch) -> None:
+    module = load_script()
+    scope = module.GeographicScope(
+        key="test-region",
+        display_name="Test region",
+        project_name="Test",
+        project_region="Test",
+        area_name="Test area",
+        authority_name="Test",
+        authority_url="https://example.test",
+        scope_type="test",
+        limitation_message="Test",
+        members=(module.ScopeMember("Mol", "13025"),),
+    )
+    boundary = box(5.0, 51.0, 5.1, 51.1)
+    payload = {
+        "type": "FeatureCollection",
+        "features": [source_feature("bodemtypes.1", box(4.98, 51.02, 5.05, 51.08))],
+    }
+    raw_bytes = json.dumps(payload, separators=(",", ":")).encode("utf-8")
+    monkeypatch.setattr(
+        module.soil,
+        "iter_wfs_pages",
+        lambda *_args, **_kwargs: iter([(payload, "https://example.test/page", raw_bytes)]),
+    )
+
+    manifest = module.prepare_partition(
+        object(),
+        output_root=tmp_path,
+        scope=scope,
+        member=scope.members[0],
+        boundary_wgs84=boundary,
+        page_limit=100,
+        max_features=1000,
+        timeout=30,
+        force=False,
+    )
+    output = json.loads(Path(manifest["output_path"]).read_text(encoding="utf-8"))
+    raw_path = Path(manifest["manifest_path"]).parent / manifest["raw_pages"][0]["artifact_path"]
+
+    assert manifest["feature_count"] == 1
+    assert output["features"][0]["id"] == "bodemtypes.1:13025"
+    assert gzip.decompress(raw_path.read_bytes()) == raw_bytes
+
+    cached = module.prepare_partition(
+        object(),
+        output_root=tmp_path,
+        scope=scope,
+        member=scope.members[0],
+        boundary_wgs84=boundary,
+        page_limit=100,
+        max_features=1000,
+        timeout=30,
+        force=False,
+    )
+    assert cached["output_sha256"] == manifest["output_sha256"]
+
+
+def test_snapshot_assembles_all_partitions_with_governed_area_metrics(tmp_path: Path) -> None:
+    module = load_script()
+    scope = module.GeographicScope(
+        key="test-region",
+        display_name="Test region",
+        project_name="Test",
+        project_region="Test",
+        area_name="Test area",
+        authority_name="Test",
+        authority_url="https://example.test",
+        scope_type="test",
+        limitation_message="Test",
+        members=(module.ScopeMember("Left", "10001"), module.ScopeMember("Right", "10002")),
+    )
+    partitions = []
+    for index, member in enumerate(scope.members):
+        output_path, manifest_path, _raw_dir = module.partition_paths(tmp_path / scope.key, member.nis_code)
+        feature = source_feature(f"bodemtypes.{index}", box(5.0 + index * 0.1, 51.0, 5.05 + index * 0.1, 51.05))
+        feature["id"] = f"bodemtypes.{index}:{member.nis_code}"
+        feature["properties"].update(
+            {
+                "clipped_area_ha": 1.0 + index,
+                "soil_generalized_legend": "Nat zand",
+                "soil_texture_class": "zand",
+                "soil_drainage_class": "nat",
+            }
+        )
+        module.soil.write_json_atomic(output_path, {"type": "FeatureCollection", "features": [feature]})
+        partitions.append(
+            {
+                "municipality": member.name,
+                "nis_code": member.nis_code,
+                "feature_count": 1,
+                "raw_source_feature_count": 1,
+                "page_count": 1,
+                "output_path": str(output_path),
+                "output_sha256": module.soil.sha256_file(output_path),
+                "manifest_path": str(manifest_path),
+            }
+        )
+
+    output_path, _manifest_path, manifest = module.assemble_snapshot(
+        output_root=tmp_path,
+        scope=scope,
+        partitions=partitions,
+        member_boundaries_sha256="boundaries-hash",
+        max_total_features=10,
+    )
+    output = json.loads(output_path.read_text(encoding="utf-8"))
+
+    assert manifest["coverage_complete"] is True
+    assert manifest["feature_count"] == 2
+    assert manifest["area_by_generalized_legend_ha"]["Nat zand"] == 3.0
+    assert manifest["area_by_texture_ha"]["zand"] == 3.0
+    assert len({feature["id"] for feature in output["features"]}) == 2
+
+
+def test_upload_contract_is_regional_historical_and_canonical(tmp_path: Path) -> None:
+    module = load_script()
+    scope = module.GEOGRAPHIC_SCOPES["kempen-transport-region"]
+    path = tmp_path / "soil.geojson"
+    path.write_text('{"type":"FeatureCollection","features":[]}', encoding="utf-8")
+    manifest_path = tmp_path / "manifest.json"
+    manifest = {
+        "coverage_complete": True,
+        "feature_count": 42,
+        "output_sha256": "output-hash",
+        "partition_identity_sha256": "partition-hash",
+        "partitions": [{} for _ in scope.members],
+        "generated_at": "2026-07-16T00:00:00+00:00",
+        "limitations": ["historical"],
+    }
+    session = FakeApiSession({"id": "dataset-id", "feature_count": 42})
+
+    result = module.upload_snapshot(
+        session,
+        base_url="http://backend:8000",
+        project_id="project-id",
+        area_id="area-id",
+        scope=scope,
+        path=path,
+        manifest_path=manifest_path,
+        manifest=manifest,
+        timeout=30,
+    )
+    data = session.calls[0][1]["data"]
+    source_metadata = json.loads(data["source_metadata_json"])
+    provenance = json.loads(data["provenance_metadata_json"])
+
+    assert result["id"] == "dataset-id"
+    assert data["area_id"] == "area-id"
+    assert data["temporal_series_key"] == "dov:digital-soil-map:kempen-transport-region"
+    assert data["valid_from"] == module.soil.VALID_FROM
+    assert data["valid_to"] == module.soil.VALID_TO
+    assert source_metadata["coverage_scope"] == "kempen-transport-region"
+    assert source_metadata["member_count"] == 28
+    assert source_metadata["authority_level"] == "authoritative_historical_baseline"
+    assert provenance["operator_tool"] == "provision_regional_soil_map.py"
+    assert provenance["raw_source_responses_retained"] is True
+
+
+def test_regional_operator_is_packaged_release_checked_and_query_optimized() -> 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")
+
+    assert "COPY scripts/provision_regional_soil_map.py" in dockerfile
+    assert "py_compile scripts/provision_regional_soil_map.py" in readiness
+    assert '"provision_regional_soil_map.py"' in service
+
+    dataset = Dataset(
+        name="regional-soil.geojson",
+        dataset_type="vector",
+        status="ready",
+        source_metadata={"partitioned_source_audit": True, "geometry_clipped_to_area": True},
+        provenance_metadata={"operator_tool": "provision_regional_soil_map.py"},
+    )
+    assert VectorFeatureService.preclipped_partition_filter(
+        dataset, "Gemeente Mol - officiele grens"
+    ) == ("municipality", "Mol")
+    assert VectorFeatureService.preclipped_partition_filter(dataset, "Vervoerregio Kempen") is None

deploy/unraid/Dockerfile.all-in-one

@@ -79,6 +79,7 @@ COPY scripts/provision_mol_flood_hazards.py /app/scripts/provision_mol_flood_haz
 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_regional_soil_map.py /app/scripts/provision_regional_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

scripts/provision_mol_soil_map.py

@@ -257,7 +257,15 @@ def iter_wfs_pages(
             break
 
 
-def normalize_feature(feature: dict[str, Any], boundary_lambert72) -> tuple[dict[str, Any] | None, bool]:
+def normalize_feature(
+    feature: dict[str, Any],
+    boundary_lambert72,
+    *,
+    municipality: str = "Mol",
+    nis_code: str = "13025",
+    coverage_scope: str = "municipality",
+    feature_id_suffix: str | None = None,
+) -> tuple[dict[str, Any] | None, bool]:
     geometry_payload = feature.get("geometry")
     if not geometry_payload:
         return None, False
@@ -279,18 +287,19 @@ def normalize_feature(feature: dict[str, Any], boundary_lambert72) -> tuple[dict
     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}")
+    persisted_id = f"{source_id}:{feature_id_suffix}" if feature_id_suffix else source_id
     properties = {
         "source_name": SOURCE_NAME,
         "source_collection": TYPE_NAME,
-        "source_feature_id": source_id,
+        "source_feature_id": persisted_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",
+        "coverage_scope": coverage_scope,
+        "municipality": municipality,
+        "nis_code": nis_code,
         "source_version": SOURCE_VERSION,
         "survey_period": SURVEY_PERIOD,
         "soil_type_code": raw.get("Bodemtype"),
@@ -317,7 +326,7 @@ def normalize_feature(feature: dict[str, Any], boundary_lambert72) -> tuple[dict
     }
     return {
         "type": "Feature",
-        "id": source_id,
+        "id": persisted_id,
         "geometry": mapping(clipped_wgs84),
         "properties": properties,
     }, was_clipped

scripts/provision_regional_soil_map.py

@@ -0,0 +1,619 @@
+"""Provision the official DOV digital soil map for an approved region.
+
+The DOV WFS is queried and audited per municipality. Source responses are
+retained as deterministic gzip artifacts, polygons are clipped in EPSG:31370,
+and one complete regional GeoJSON snapshot is imported through DatasetService.
+The 1949-1971 survey remains a historical baseline, never a current drainage
+or parcel-scale soil investigation.
+"""
+
+from __future__ import annotations
+
+import argparse
+from collections import defaultdict
+import gzip
+import json
+import os
+from pathlib import Path
+import sys
+from typing import Any
+
+import requests
+from shapely.geometry import mapping, shape
+from shapely.ops import transform as transform_geometry
+
+from geographic_scopes import GEOGRAPHIC_SCOPES, GeographicScope, ScopeMember
+import provision_mol_soil_map as soil
+
+
+DEFAULT_SCOPE_KEY = "kempen-transport-region"
+DEFAULT_SCOPE_OUTPUT_ROOT = Path("/app/storage/operator-data/geographic-scopes")
+DEFAULT_OUTPUT_ROOT = Path("/app/storage/operator-evidence/dov-soil-map/regional")
+SCHEMA_VERSION = 1
+
+
+def parse_args() -> argparse.Namespace:
+    parser = argparse.ArgumentParser(description="Provision regional official DOV soil reference data.")
+    parser.add_argument("--scope", choices=sorted(GEOGRAPHIC_SCOPES), default=DEFAULT_SCOPE_KEY)
+    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", soil.DEFAULT_API_URL))
+    parser.add_argument(
+        "--scope-output-root",
+        type=Path,
+        default=Path(os.environ.get("GEOINTEL_SCOPE_OUTPUT_ROOT", DEFAULT_SCOPE_OUTPUT_ROOT)),
+    )
+    parser.add_argument(
+        "--output-root",
+        type=Path,
+        default=Path(os.environ.get("GEOINTEL_REGIONAL_SOIL_OUTPUT_ROOT", DEFAULT_OUTPUT_ROOT)),
+    )
+    parser.add_argument("--page-limit", type=int, default=1000)
+    parser.add_argument("--max-features-per-partition", type=int, default=20_000)
+    parser.add_argument("--max-total-features", type=int, default=200_000)
+    parser.add_argument("--request-timeout", type=int, default=300)
+    parser.add_argument("--import-timeout", type=int, default=3600)
+    parser.add_argument("--fetch-only", action="store_true")
+    parser.add_argument("--force", action="store_true")
+    return parser.parse_args()
+
+
+def boundary_sha256(geometry) -> str:
+    encoded = json.dumps(mapping(geometry), sort_keys=True, separators=(",", ":")).encode("utf-8")
+    return soil.sha256_bytes(encoded)
+
+
+def resolve_member_boundaries(scope: GeographicScope, scope_output_root: Path) -> tuple[Path, str]:
+    scope_dir = scope_output_root / scope.key
+    manifest_path = scope_dir / f"{scope.key.replace('-', '_')}_scope_manifest.json"
+    if not manifest_path.is_file():
+        raise RuntimeError(f"Official scope manifest is missing at {manifest_path}")
+    manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
+    if (
+        manifest.get("status") != "complete"
+        or manifest.get("scope_key") != scope.key
+        or int(manifest.get("member_count") or 0) != len(scope.members)
+    ):
+        raise RuntimeError(f"Official scope manifest at {manifest_path} is incomplete or inconsistent")
+    members_path = scope_dir / str(manifest.get("municipalities_filename") or "")
+    expected_sha256 = str(manifest.get("municipalities_sha256") or "")
+    if not members_path.is_file() or not expected_sha256 or soil.sha256_file(members_path) != expected_sha256:
+        raise RuntimeError("Official municipality-boundary artifact is missing or fails its scope checksum")
+    return members_path, expected_sha256
+
+
+def load_member_boundaries(path: Path, scope: GeographicScope) -> dict[str, tuple[ScopeMember, Any]]:
+    payload = json.loads(path.read_text(encoding="utf-8"))
+    features = payload.get("features") if isinstance(payload, dict) else None
+    if not isinstance(features, list):
+        raise RuntimeError("Municipality-boundary artifact is not a GeoJSON FeatureCollection")
+    expected = {member.nis_code: member for member in scope.members}
+    selected: dict[str, tuple[ScopeMember, Any]] = {}
+    for feature in features:
+        properties = feature.get("properties") or {}
+        nis_code = str(properties.get("nis_code") or properties.get("NISCODE") or "")
+        if nis_code not in expected:
+            continue
+        if nis_code in selected:
+            raise RuntimeError(f"Municipality-boundary artifact contains duplicate NIS code {nis_code}")
+        geometry = soil.polygonal_geometry(shape(feature.get("geometry")))
+        if geometry is None:
+            raise RuntimeError(f"Municipality boundary for {expected[nis_code].name} is invalid")
+        selected[nis_code] = (expected[nis_code], geometry)
+    missing = [member.name for member in scope.members if member.nis_code not in selected]
+    if missing:
+        raise RuntimeError(f"Municipality-boundary artifact is missing: {', '.join(missing)}")
+    return {member.nis_code: selected[member.nis_code] for member in scope.members}
+
+
+def partition_paths(output_root: Path, nis_code: str) -> tuple[Path, Path, Path]:
+    partition_dir = output_root / "partitions" / nis_code
+    return partition_dir / "dov_soil_map.geojson", partition_dir / "manifest.json", partition_dir / "raw"
+
+
+def cached_partition(
+    output_root: Path,
+    member: ScopeMember,
+    boundary_hash: str,
+) -> dict[str, Any] | None:
+    output_path, manifest_path, _raw_dir = partition_paths(output_root, member.nis_code)
+    if not output_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") != soil.SOURCE_VERSION
+        or manifest.get("nis_code") != member.nis_code
+        or manifest.get("boundary_sha256") != boundary_hash
+        or manifest.get("output_sha256") != soil.sha256_file(output_path)
+    ):
+        return None
+    for page in manifest.get("raw_pages") or []:
+        page_path = manifest_path.parent / str(page.get("artifact_path") or "")
+        if not page_path.is_file() or page.get("artifact_sha256") != soil.sha256_file(page_path):
+            return None
+        try:
+            raw_bytes = gzip.decompress(page_path.read_bytes())
+        except (OSError, EOFError):
+            return None
+        if page.get("response_sha256") != soil.sha256_bytes(raw_bytes):
+            return None
+    manifest["output_path"] = str(output_path)
+    manifest["manifest_path"] = str(manifest_path)
+    return manifest
+
+
+def prepare_partition(
+    session,
+    *,
+    output_root: Path,
+    scope: GeographicScope,
+    member: ScopeMember,
+    boundary_wgs84,
+    page_limit: int,
+    max_features: int,
+    timeout: int,
+    force: bool,
+) -> dict[str, Any]:
+    output_path, manifest_path, raw_dir = partition_paths(output_root, member.nis_code)
+    current_boundary_hash = boundary_sha256(boundary_wgs84)
+    if not force:
+        reusable = cached_partition(output_root, member, current_boundary_hash)
+        if reusable:
+            return reusable
+
+    raw_dir.mkdir(parents=True, exist_ok=True)
+    boundary_lambert72 = soil.polygonal_geometry(transform_geometry(soil.TO_LAMBERT72.transform, boundary_wgs84))
+    if boundary_lambert72 is None:
+        raise RuntimeError(f"Boundary for {member.name} could not be transformed to EPSG:31370")
+
+    retained: list[dict[str, Any]] = []
+    raw_pages: list[dict[str, Any]] = []
+    seen_source_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(
+        soil.iter_wfs_pages(
+            session,
+            boundary_lambert72.bounds,
+            page_limit=page_limit,
+            timeout=timeout,
+        ),
+        start=1,
+    ):
+        page_features = list(payload.get("features") or [])
+        raw_feature_count += len(page_features)
+        if raw_feature_count > max_features:
+            raise RuntimeError(
+                f"DOV WFS returned more than {max_features} features for {member.name}; refusing truncation"
+            )
+        compressed = gzip.compress(raw_bytes, mtime=0)
+        raw_path = raw_dir / f"dov_soil_map_page_{page_number:05d}.json.gz"
+        soil.write_bytes_atomic(raw_path, compressed)
+        raw_pages.append(
+            {
+                "artifact_path": str(raw_path.relative_to(manifest_path.parent)),
+                "artifact_sha256": soil.sha256_bytes(compressed),
+                "response_sha256": soil.sha256_bytes(raw_bytes),
+                "response_size_bytes": len(raw_bytes),
+                "feature_count": len(page_features),
+                "source_url": source_url,
+            }
+        )
+        for source_feature in page_features:
+            raw = dict(source_feature.get("properties") or {})
+            source_id = str(source_feature.get("id") or f"{soil.TYPE_NAME}:{raw.get('gid')}")
+            if source_id in seen_source_ids:
+                duplicate_count += 1
+                continue
+            seen_source_ids.add(source_id)
+            normalized, was_clipped = soil.normalize_feature(
+                source_feature,
+                boundary_lambert72,
+                municipality=member.name,
+                nis_code=member.nis_code,
+                coverage_scope=scope.key,
+                feature_id_suffix=member.nis_code,
+            )
+            if normalized is None:
+                rejected_count += 1
+                continue
+            clipped_count += int(was_clipped)
+            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(f"DOV WFS returned no valid soil polygons inside {member.name}")
+    artifact = {
+        "type": "FeatureCollection",
+        "name": f"Digitale bodemkaart - {member.name}",
+        "features": retained,
+        "source": soil.ATTRIBUTION,
+        "source_version": soil.SOURCE_VERSION,
+        "survey_period": soil.SURVEY_PERIOD,
+        "coverage_scope": scope.key,
+        "municipality": member.name,
+        "nis_code": member.nis_code,
+        "reference_truncated": False,
+    }
+    soil.write_json_atomic(output_path, artifact)
+    manifest = {
+        "schema_version": SCHEMA_VERSION,
+        "source_version": soil.SOURCE_VERSION,
+        "survey_period": soil.SURVEY_PERIOD,
+        "scope_key": scope.key,
+        "municipality": member.name,
+        "nis_code": member.nis_code,
+        "boundary_sha256": current_boundary_hash,
+        "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,
+        "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())},
+        "output_path": str(output_path),
+        "output_sha256": soil.sha256_file(output_path),
+        "generated_at": soil.utc_now(),
+    }
+    soil.write_json_atomic(manifest_path, manifest, pretty=True)
+    manifest["manifest_path"] = str(manifest_path)
+    return manifest
+
+
+def snapshot_paths(output_root: Path, scope: GeographicScope) -> tuple[Path, Path]:
+    snapshot_dir = output_root / "snapshots" / scope.key
+    return snapshot_dir / "dov_soil_map.geojson", snapshot_dir / "manifest.json"
+
+
+def assemble_snapshot(
+    *,
+    output_root: Path,
+    scope: GeographicScope,
+    partitions: list[dict[str, Any]],
+    member_boundaries_sha256: str,
+    max_total_features: int,
+) -> tuple[Path, Path, dict[str, Any]]:
+    if [item.get("nis_code") for item in partitions] != list(scope.nis_codes):
+        raise RuntimeError("Soil partition order/completeness does not match the approved geographic scope")
+    output_path, manifest_path = snapshot_paths(output_root, scope)
+    partition_identity = soil.sha256_bytes(
+        json.dumps(
+            [(item["nis_code"], item["output_sha256"]) for item in partitions],
+            separators=(",", ":"),
+        ).encode("utf-8")
+    )
+    if output_path.is_file() and manifest_path.is_file():
+        try:
+            existing = json.loads(manifest_path.read_text(encoding="utf-8"))
+        except (OSError, ValueError):
+            existing = {}
+        if (
+            existing.get("partition_identity_sha256") == partition_identity
+            and existing.get("output_sha256") == soil.sha256_file(output_path)
+        ):
+            return output_path, manifest_path, existing
+
+    features: list[dict[str, Any]] = []
+    feature_ids: set[str] = set()
+    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 partition in partitions:
+        payload = json.loads(Path(str(partition["output_path"])).read_text(encoding="utf-8"))
+        partition_features = payload.get("features") if isinstance(payload, dict) else None
+        if not isinstance(partition_features, list) or len(partition_features) != int(partition["feature_count"]):
+            raise RuntimeError(f"Soil partition feature-count drift for NIS {partition['nis_code']}")
+        for feature in partition_features:
+            feature_id = str(feature.get("id") or "")
+            if not feature_id or feature_id in feature_ids:
+                raise RuntimeError(f"Duplicate or missing regional soil feature id {feature_id!r}")
+            feature_ids.add(feature_id)
+            properties = feature.get("properties") or {}
+            area = float(properties.get("clipped_area_ha") or 0.0)
+            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
+            features.append(feature)
+            if len(features) > max_total_features:
+                raise RuntimeError(
+                    f"Regional soil snapshot exceeds the {max_total_features} feature safety limit; refusing truncation"
+                )
+
+    generated_at = soil.utc_now()
+    artifact = {
+        "type": "FeatureCollection",
+        "name": f"Digitale bodemkaart - {scope.display_name}",
+        "features": features,
+        "source": soil.ATTRIBUTION,
+        "source_version": soil.SOURCE_VERSION,
+        "survey_period": soil.SURVEY_PERIOD,
+        "catalog_url": soil.CATALOG_URL,
+        "coverage_scope": scope.key,
+        "member_count": len(scope.members),
+        "reference_truncated": False,
+        "generated_at": generated_at,
+    }
+    soil.write_json_atomic(output_path, artifact)
+    manifest = {
+        "schema_version": SCHEMA_VERSION,
+        "source_version": soil.SOURCE_VERSION,
+        "survey_period": soil.SURVEY_PERIOD,
+        "scope_key": scope.key,
+        "scope_display_name": scope.display_name,
+        "member_count": len(scope.members),
+        "member_nis_codes": list(scope.nis_codes),
+        "member_boundaries_sha256": member_boundaries_sha256,
+        "coverage_complete": len(partitions) == len(scope.members),
+        "feature_count": len(features),
+        "raw_source_feature_count": sum(int(item["raw_source_feature_count"]) for item in partitions),
+        "raw_response_count": sum(int(item["page_count"]) for item in partitions),
+        "partition_identity_sha256": partition_identity,
+        "partitions": [
+            {
+                "municipality": item["municipality"],
+                "nis_code": item["nis_code"],
+                "feature_count": item["feature_count"],
+                "raw_source_feature_count": item["raw_source_feature_count"],
+                "output_sha256": item["output_sha256"],
+                "manifest_path": item["manifest_path"],
+            }
+            for item in partitions
+        ],
+        "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())},
+        "reference_truncated": False,
+        "output_sha256": soil.sha256_file(output_path),
+        "output_size_bytes": output_path.stat().st_size,
+        "generated_at": generated_at,
+        "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.",
+            "Municipality partitions split source polygons at administrative boundaries; source ids carry a NIS suffix.",
+        ],
+    }
+    soil.write_json_atomic(manifest_path, manifest, pretty=True)
+    return output_path, manifest_path, manifest
+
+
+def locate_workspace(session, base_url: str, scope: GeographicScope, timeout: int):
+    projects = soil.paginated_items(session, f"{base_url}/api/v1/projects", timeout=timeout)
+    project = next((item for item in projects if item.get("name") == scope.project_name), None)
+    if not project:
+        raise RuntimeError(f"Project {scope.project_name!r} is missing")
+    project_id = str(project["id"])
+    areas = soil.paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", timeout=timeout)
+    area = next((item for item in areas if item.get("name") == scope.area_name), None)
+    if not area:
+        raise RuntimeError(f"Official scope Area {scope.area_name!r} is missing")
+    datasets = soil.paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/datasets", timeout=timeout)
+    return project_id, str(area["id"]), datasets
+
+
+def upload_snapshot(
+    session,
+    *,
+    base_url: str,
+    project_id: str,
+    area_id: str,
+    scope: GeographicScope,
+    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": soil.SOURCE_NAME,
+        "theme": "soil",
+        "layer_type": "soil",
+        "source_collection": soil.TYPE_NAME,
+        "source_crs": "EPSG:31370",
+        "persisted_crs": "EPSG:4326",
+        "authority_level": "authoritative_historical_baseline",
+        "coverage_scope": scope.key,
+        "scope_type": scope.scope_type,
+        "scope_display_name": scope.display_name,
+        "member_count": len(scope.members),
+        "member_nis_codes": list(scope.nis_codes),
+        "partitioned_source_audit": True,
+        "coverage_complete": bool(manifest["coverage_complete"]),
+        "feature_count": manifest["feature_count"],
+        "geometry_clipped_to_area": True,
+        "semantic_metrics": False,
+        "survey_period": soil.SURVEY_PERIOD,
+        "source_scale": "1:20,000",
+        "attribution": soil.ATTRIBUTION,
+        "catalog_url": soil.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": soil.selection_metrics(),
+    }
+    provenance_metadata = {
+        "operator_tool": "provision_regional_soil_map.py",
+        "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
+        "source_type_name": soil.TYPE_NAME,
+        "wfs_url": soil.WFS_URL,
+        "catalog_url": soil.CATALOG_URL,
+        "manifest_path": str(manifest_path),
+        "combined_output_sha256": manifest["output_sha256"],
+        "partition_count": len(manifest["partitions"]),
+        "partition_identity_sha256": manifest["partition_identity_sha256"],
+        "raw_source_responses_retained": True,
+        "reference_truncated": False,
+        "generated_at": manifest["generated_at"],
+        "limitations": manifest["limitations"],
+    }
+    with path.open("rb") as handle:
+        response = session.post(
+            f"{base_url}/api/v1/projects/{project_id}/datasets/upload",
+            data={
+                "dataset_type": "vector",
+                "source": "operator_official_import",
+                "dataset_role": "reference",
+                "source_name": soil.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": f"dov:digital-soil-map:{scope.key}",
+                "observed_at": soil.OBSERVED_AT,
+                "valid_from": soil.VALID_FROM,
+                "valid_to": soil.VALID_TO,
+                "temporal_granularity": "period",
+                "source_version": soil.SOURCE_VERSION,
+            },
+            files={"file": (path.name, handle, "application/geo+json")},
+            timeout=timeout,
+        )
+    return soil.response_data(response)
+
+
+def main() -> int:
+    args = parse_args()
+    try:
+        if (
+            args.page_limit < 1
+            or args.page_limit > 5000
+            or args.max_features_per_partition < args.page_limit
+            or args.max_total_features < args.max_features_per_partition
+        ):
+            raise ValueError("Invalid soil-map page or feature safety limits")
+        scope = GEOGRAPHIC_SCOPES[args.scope]
+        members_path, members_sha256 = resolve_member_boundaries(scope, args.scope_output_root)
+        boundaries = load_member_boundaries(members_path, scope)
+        with soil.source_session() as source_session:
+            source_session.headers.update({"User-Agent": "GeoIntel-DOV-Soil-Regional-Operator/1.0"})
+            partitions = [
+                prepare_partition(
+                    source_session,
+                    output_root=args.output_root / scope.key,
+                    scope=scope,
+                    member=member,
+                    boundary_wgs84=boundary,
+                    page_limit=args.page_limit,
+                    max_features=args.max_features_per_partition,
+                    timeout=args.request_timeout,
+                    force=args.force,
+                )
+                for member, boundary in boundaries.values()
+            ]
+        path, manifest_path, manifest = assemble_snapshot(
+            output_root=args.output_root,
+            scope=scope,
+            partitions=partitions,
+            member_boundaries_sha256=members_sha256,
+            max_total_features=args.max_total_features,
+        )
+        if args.fetch_only:
+            result: dict[str, Any] = {
+                "status": "prepared",
+                "artifact_path": str(path),
+                "feature_count": manifest["feature_count"],
+            }
+        else:
+            with requests.Session() as api_session:
+                project_id, area_id, datasets = locate_workspace(
+                    api_session,
+                    args.base_url.rstrip("/"),
+                    scope,
+                    args.import_timeout,
+                )
+                existing = next(
+                    (
+                        item
+                        for item in datasets
+                        if item.get("source_name") == soil.SOURCE_NAME
+                        and item.get("source_version") == soil.SOURCE_VERSION
+                        and str(item.get("area_id") or "") == area_id
+                    ),
+                    None,
+                )
+                if existing:
+                    if str(existing.get("checksum_sha256") or "") != manifest["output_sha256"]:
+                        raise RuntimeError("A different regional DOV soil artifact already exists; refusing replacement")
+                    result = {
+                        "status": "existing",
+                        "dataset_id": existing["id"],
+                        "feature_count": existing.get("feature_count"),
+                    }
+                else:
+                    dataset = upload_snapshot(
+                        api_session,
+                        base_url=args.base_url.rstrip("/"),
+                        project_id=project_id,
+                        area_id=area_id,
+                        scope=scope,
+                        path=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",
+                "scope": scope.key,
+                "member_count": len(scope.members),
+                "partition_count": len(partitions),
+                "feature_count": manifest["feature_count"],
+                "raw_source_feature_count": manifest["raw_source_feature_count"],
+                "raw_response_count": manifest["raw_response_count"],
+                "manifest_path": str(manifest_path),
+                "metrics": {
+                    "area_by_generalized_legend_ha": manifest["area_by_generalized_legend_ha"],
+                    "area_by_texture_ha": manifest["area_by_texture_ha"],
+                    "area_by_drainage_ha": manifest["area_by_drainage_ha"],
+                },
+                "result": result,
+            },
+            ensure_ascii=False,
+            indent=2,
+        )
+    )
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())

scripts/run_readiness_check.sh

@@ -58,6 +58,7 @@ ${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_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