"""Provision an audited BWK/Natura 2000 snapshot for an approved region. The official WFS is queried per municipality. Every source response is retained as deterministic gzip evidence, geometry is clipped in EPSG:31370, and one regional GeoJSON snapshot is persisted through the canonical DatasetService API. """ 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_bwk_natura2000 as bwk 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/bwk-natura2000-2025/regional") SCHEMA_VERSION = 1 def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description="Provision regional official BWK/Natura 2000 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", bwk.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_BWK_OUTPUT_ROOT", DEFAULT_OUTPUT_ROOT)), ) parser.add_argument("--page-limit", type=int, default=1000) parser.add_argument("--max-features-per-partition", type=int, default=30_000) parser.add_argument("--max-total-features", type=int, default=300_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: return bwk.sha256_bytes(json.dumps(mapping(geometry), sort_keys=True, separators=(",", ":")).encode("utf-8")) 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 bwk.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 = bwk.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 / "bwk_natura2000_2025.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") != bwk.SOURCE_VERSION or manifest.get("nis_code") != member.nis_code or manifest.get("boundary_sha256") != boundary_hash or manifest.get("output_sha256") != bwk.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") != bwk.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") != bwk.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 = bwk.polygonal_geometry(transform_geometry(bwk.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 for page_number, (payload, source_url, raw_bytes) in enumerate( bwk.iter_wfs_pages(session, boundary_wgs84.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"BWK WFS returned more than {max_features} features for {member.name}; refusing a truncated partition" ) compressed = gzip.compress(raw_bytes, mtime=0) raw_path = raw_dir / f"bwk_bwkhab_page_{page_number:05d}.json.gz" bwk.write_bytes_atomic(raw_path, compressed) raw_pages.append( { "artifact_path": str(raw_path.relative_to(manifest_path.parent)), "artifact_sha256": bwk.sha256_bytes(compressed), "response_sha256": bwk.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: source_id = str(source_feature.get("id") or (source_feature.get("properties") or {}).get("UIDN") or "") if source_id and source_id in seen_source_ids: duplicate_count += 1 continue if source_id: seen_source_ids.add(source_id) normalized, was_clipped = bwk.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) artifact = { "type": "FeatureCollection", "name": f"BWK en Natura 2000 - toestand 2025 - {member.name}", "features": retained, "source": "INBO BWK/Natura 2000 WFS", "source_version": bwk.SOURCE_VERSION, "coverage_scope": scope.key, "municipality": member.name, "nis_code": member.nis_code, "reference_truncated": False, } bwk.write_json_atomic(output_path, artifact) manifest = { "schema_version": SCHEMA_VERSION, "source_version": bwk.SOURCE_VERSION, "scope_key": scope.key, "municipality": member.name, "nis_code": member.nis_code, "boundary_sha256": current_boundary_hash, "boundary_bbox_wgs84": list(boundary_wgs84.bounds), "page_limit": page_limit, "page_count": len(raw_pages), "raw_source_feature_count": raw_feature_count, "feature_count": len(retained), "duplicate_count": duplicate_count, "rejected_or_outside_count": rejected_count, "clipped_feature_count": clipped_count, "reference_truncated": False, "raw_pages": raw_pages, "output_path": str(output_path), "output_sha256": bwk.sha256_file(output_path), "generated_at": bwk.utc_now(), } bwk.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 / "bwk_natura2000_2025.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("BWK partition order/completeness does not match the approved geographic scope") output_path, manifest_path = snapshot_paths(output_root, scope) partition_identity = bwk.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") == bwk.sha256_file(output_path) ): return output_path, manifest_path, existing features: list[dict[str, Any]] = [] feature_ids: set[str] = set() evaluation_area: dict[str, float] = defaultdict(float) habitat_status_area: dict[str, float] = defaultdict(float) natura_area = 0.0 regional_area = 0.0 uncertain_area = 0.0 for 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"BWK 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 BWK 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) evaluation_area[str(properties.get("bwk_evaluation_code") or "unknown")] += area habitat_status_area[str(properties.get("habitat_status_code") or "unknown")] += area natura_area += float(properties.get("natura2000_area_ha") or 0.0) regional_area += float(properties.get("regional_biotope_area_ha") or 0.0) uncertain_area += float(properties.get("uncertain_habitat_area_ha") or 0.0) features.append(feature) if len(features) > max_total_features: raise RuntimeError( f"Regional BWK snapshot exceeds the {max_total_features} feature safety limit; refusing truncation" ) generated_at = bwk.utc_now() artifact = { "type": "FeatureCollection", "name": f"BWK en Natura 2000 - toestand 2025 - {scope.display_name}", "features": features, "source": "INBO BWK/Natura 2000 WFS", "source_version": bwk.SOURCE_VERSION, "attribution": bwk.ATTRIBUTION, "catalog_url": bwk.CATALOG_URL, "report_url": bwk.REPORT_URL, "coverage_scope": scope.key, "member_count": len(scope.members), "reference_truncated": False, "generated_at": generated_at, } bwk.write_json_atomic(output_path, artifact) empty_partitions = [item["nis_code"] for item in partitions if int(item["feature_count"]) == 0] manifest = { "schema_version": SCHEMA_VERSION, "source_version": bwk.SOURCE_VERSION, "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), "empty_partitions": empty_partitions, "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 ], "evaluation_area_ha": {key: round(value, 6) for key, value in sorted(evaluation_area.items())}, "habitat_status_area_ha": {key: round(value, 6) for key, value in sorted(habitat_status_area.items())}, "natura2000_area_ha": round(natura_area, 6), "regional_biotope_area_ha": round(regional_area, 6), "uncertain_habitat_area_ha": round(uncertain_area, 6), "reference_truncated": False, "output_sha256": bwk.sha256_file(output_path), "output_size_bytes": output_path.stat().st_size, "generated_at": generated_at, "limitations": [ "The 2025 edition is the best available map state, not one uniform 2025 field survey.", "PHAB percentages may be theoretical shares; partial selections scale shares by intersected polygon area.", "Municipality partitions split source polygons at administrative boundaries; source ids carry a NIS suffix.", "The real field situation remains authoritative for policy and legal use.", ], } bwk.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 = bwk.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 = bwk.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 = bwk.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]: source_metadata = { "provider": "Instituut voor Natuur- en Bosonderzoek", "theme": "nature_value", "layer_name": "BWK/Natura 2000 toestand 2025", "authority_level": "authoritative", "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, "attribution": bwk.ATTRIBUTION, "catalog_url": bwk.CATALOG_URL, "report_url": bwk.REPORT_URL, "selection_aggregation": { "metric_key": "bwk_mapped_area", "method": "intersection_area", "label": "BWK-gekarteerde oppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": "Uitgave 2025 is de best beschikbare kaarttoestand, maar niet elk kaartvlak is in 2025 op terrein gekarteerd.", }, "selection_metrics": bwk.selection_metrics(), } provenance_metadata = { "operator_tool": "provision_regional_bwk_natura2000.py", "operator_explicit_fetch": True, "geometry_clipped_to_area": True, "source_type_name": bwk.TYPE_NAME, "wfs_url": bwk.WFS_URL, "catalog_url": bwk.CATALOG_URL, "report_url": bwk.REPORT_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": "inbo_bwk_natura2000", "reference_layer_name": "nature_value", "source_metadata_json": json.dumps(source_metadata, ensure_ascii=False), "provenance_metadata_json": json.dumps(provenance_metadata, ensure_ascii=False), "area_id": area_id, "temporal_series_key": f"inbo-bwk-natura2000:{scope.key}", "observed_at": bwk.PUBLICATION_DATE, "temporal_granularity": "snapshot", "source_version": bwk.SOURCE_VERSION, }, files={"file": (path.name, handle, "application/geo+json")}, timeout=timeout, ) return bwk.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 BWK 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) source_session = bwk.source_session() source_session.headers.update({"User-Agent": "GeoIntel-BWK-Natura2000-Regional-Operator/1.0"}) with source_session: 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, ) result: dict[str, Any] if args.fetch_only: result = {"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") == "inbo_bwk_natura2000" and item.get("source_version") == bwk.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 BWK 2025 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": { "evaluation_area_ha": manifest["evaluation_area_ha"], "natura2000_area_ha": manifest["natura2000_area_ha"], "regional_biotope_area_ha": manifest["regional_biotope_area_ha"], "uncertain_habitat_area_ha": manifest["uncertain_habitat_area_ha"], }, "result": result, }, ensure_ascii=False, indent=2, ) ) return 0 if __name__ == "__main__": sys.exit(main())