"""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())