Files
geointel/scripts/provision_regional_bwk_natura2000.py
Jens faeb58ef6d
GeoIntel release gates / Compile, test, contracts and builds (push) Successful in 1m49s
GeoIntel release gates / Python and npm vulnerability policy (push) Successful in 21s
GeoIntel release gates / Production AI image, SBOM and container scan (push) Successful in 5m39s
GeoIntel release gates / Deploy exact gated revision to Unraid (push) Failing after 58m43s
Initial public release
2026-08-31 21:56:53 +02:00

588 lines
26 KiB
Python

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