diff --git a/backend/app/services/vector_feature_service.py b/backend/app/services/vector_feature_service.py index 3a119ddb..4c23d4b6 100644 --- a/backend/app/services/vector_feature_service.py +++ b/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", } diff --git a/backend/tests/test_sprint217_regional_dov_soil_map.py b/backend/tests/test_sprint217_regional_dov_soil_map.py new file mode 100644 index 00000000..da342022 --- /dev/null +++ b/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 diff --git a/deploy/unraid/Dockerfile.all-in-one b/deploy/unraid/Dockerfile.all-in-one index 0b190c90..e7f2a329 100644 --- a/deploy/unraid/Dockerfile.all-in-one +++ b/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 diff --git a/scripts/provision_mol_soil_map.py b/scripts/provision_mol_soil_map.py index a738920d..d07bae60 100644 --- a/scripts/provision_mol_soil_map.py +++ b/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 diff --git a/scripts/provision_regional_soil_map.py b/scripts/provision_regional_soil_map.py new file mode 100644 index 00000000..faeb51ee --- /dev/null +++ b/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()) diff --git a/scripts/run_readiness_check.sh b/scripts/run_readiness_check.sh index 77e2e15c..5e007fa6 100755 --- a/scripts/run_readiness_check.sh +++ b/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