"""Provision partition-audited historical land use for an approved region. The official historical land-use WFS caps broad regional result counts. This operator therefore fetches and clips each approved municipality separately, retains the exact source responses as checksummed gzip artifacts, assembles one regional snapshot per theme/year and persists only through the dataset API. """ from __future__ import annotations import argparse import gzip import hashlib import json import os import sys from datetime import datetime, timezone from pathlib import Path from typing import Any import requests from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, mapping, shape from shapely.ops import unary_union from shapely.validation import make_valid from geographic_scopes import GEOGRAPHIC_SCOPES, GeographicScope, ScopeMember from provision_mol_historical_landuse import ( ATTRIBUTION, COLLECTIONS, THEMES, WFS_URL, ThemeDefinition, build_session, response_data, wfs_filter_xml, wfs_request_xml, ) DEFAULT_API_URL = "http://127.0.0.1:8000" DEFAULT_SCOPE_KEY = "kempen-transport-region" DEFAULT_SCOPE_OUTPUT_ROOT = Path("/app/storage/operator-data/geographic-scopes") DEFAULT_OUTPUT_ROOT = Path("/app/storage/operator-data/regional-historical-landuse") SOURCE_CATALOG_URL = ( "https://www.vlaanderen.be/datavindplaats/catalogus/" "digitalisatie-historisch-landgebruik-en-landgebruiksveranderingen-in-vlaanderen-1778-2022" ) SUPPORTED_THEME_KEYS = frozenset({"buildings", "water", "roads"}) GEOJSON_CRS = {"type": "name", "properties": {"name": "EPSG:4326"}} def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description="Provision regional historical buildings, water and roads.") 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", DEFAULT_API_URL)) parser.add_argument("--years", default="1778,1873,1969") parser.add_argument("--themes", default="buildings,water,roads") 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_HISTORICAL_OUTPUT_ROOT", DEFAULT_OUTPUT_ROOT)), ) parser.add_argument("--page-size", type=int, default=500) parser.add_argument("--max-features-per-partition", type=int, default=100_000) parser.add_argument("--max-total-features", type=int, default=500_000) parser.add_argument("--simplify-tolerance-degrees", type=float, default=0.00001) parser.add_argument("--request-timeout", type=int, default=300) parser.add_argument("--import-timeout", type=int, default=3600) parser.add_argument("--force", action="store_true") parser.add_argument("--fetch-only", action="store_true") return parser.parse_args() def utc_now() -> str: return datetime.now(timezone.utc).isoformat() def sha256_bytes(content: bytes) -> str: return hashlib.sha256(content).hexdigest() def sha256_file(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: for chunk in iter(lambda: handle.read(1024 * 1024), b""): digest.update(chunk) return digest.hexdigest() def atomic_write_bytes(path: Path, content: bytes) -> None: path.parent.mkdir(parents=True, exist_ok=True) temporary = path.with_suffix(f"{path.suffix}.partial") temporary.write_bytes(content) temporary.replace(path) def atomic_write_json(path: Path, payload: dict[str, Any], *, pretty: bool = False) -> None: content = json.dumps( payload, ensure_ascii=False, indent=2 if pretty else None, separators=None if pretty else (",", ":"), sort_keys=pretty, ).encode("utf-8") atomic_write_bytes(path, content) def normalize_polygonal(geometry): if geometry is None or geometry.is_empty: return None if not geometry.is_valid: geometry = make_valid(geometry) if isinstance(geometry, (Polygon, MultiPolygon)): return geometry if isinstance(geometry, GeometryCollection): polygons = [ part for part in geometry.geoms if isinstance(part, (Polygon, MultiPolygon)) and not part.is_empty ] if polygons: merged = unary_union(polygons) if isinstance(merged, (Polygon, MultiPolygon)) and not merged.is_empty: return merged return None def resolve_member_boundaries(scope: GeographicScope, scope_output_root: Path) -> Path: 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}; run provision_geographic_scope.py first" ) 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 "") if not members_path.is_file() or sha256_file(members_path) != manifest.get("municipalities_sha256"): raise RuntimeError("Official municipality-boundary artifact is missing or fails its scope checksum") return members_path 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 = normalize_polygonal(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 fetch_source_page( session: requests.Session, *, collection: str, filter_xml: str, page_size: int, start_index: int, timeout: int, ) -> tuple[bytes, list[dict[str, Any]]]: response = session.post( WFS_URL, data=wfs_request_xml( collection=collection, filter_xml=filter_xml, page_size=page_size, start_index=start_index, ).encode("utf-8"), headers={"Content-Type": "application/xml; charset=UTF-8", "Accept": "application/json"}, timeout=timeout, ) response.raise_for_status() raw_content = response.content payload = response.json() features = payload.get("features") if isinstance(payload, dict) else None if not isinstance(features, list): raise RuntimeError("Historical land-use WFS returned an invalid FeatureCollection") return raw_content, features def partition_paths(output_root: Path, year: int, theme: str, nis_code: str) -> tuple[Path, Path, Path]: directory = output_root / "partitions" / str(year) / theme output_path = directory / f"{nis_code}.geojson" return output_path, output_path.with_suffix(".manifest.json"), directory / f"{nis_code}.raw" def cached_partition( output_path: Path, manifest_path: Path, *, year: int, theme: str, nis_code: str, page_size: int, simplify_tolerance_degrees: float, boundary_geometry_sha256: str, ): if not output_path.is_file() or not manifest_path.is_file(): return None manifest = json.loads(manifest_path.read_text(encoding="utf-8")) if ( manifest.get("status") != "complete" or manifest.get("year") != year or manifest.get("theme") != theme or manifest.get("nis_code") != nis_code or int(manifest.get("page_size") or 0) != page_size or float(manifest.get("geometry_simplification_tolerance_degrees") or 0.0) != simplify_tolerance_degrees or manifest.get("boundary_geometry_sha256") != boundary_geometry_sha256 or sha256_file(output_path) != manifest.get("output_sha256") ): return None for page in manifest.get("raw_pages") or []: raw_path = manifest_path.parent / str(page.get("artifact_path") or "") if not raw_path.is_file() or sha256_file(raw_path) != page.get("artifact_sha256"): return None return manifest def prepare_partition( session: requests.Session, *, output_root: Path, year: int, definition: ThemeDefinition, scope_key: str, member: ScopeMember, boundary, page_size: int, max_features: int, simplify_tolerance_degrees: float, timeout: int, force: bool, ) -> dict[str, Any]: output_path, manifest_path, raw_dir = partition_paths(output_root, year, definition.key, member.nis_code) boundary_geometry_sha256 = sha256_bytes(boundary.wkb) if not force: cached = cached_partition( output_path, manifest_path, year=year, theme=definition.key, nis_code=member.nis_code, page_size=page_size, simplify_tolerance_degrees=simplify_tolerance_degrees, boundary_geometry_sha256=boundary_geometry_sha256, ) if cached: return cached collection = COLLECTIONS[year] filter_xml = wfs_filter_xml(definition, boundary.bounds) source_seen: set[str] = set() output_features: list[dict[str, Any]] = [] raw_pages: list[dict[str, Any]] = [] invalid_geometry_count = 0 outside_boundary_count = 0 start_index = 0 while True: raw_content, page = fetch_source_page( session, collection=collection, filter_xml=filter_xml, page_size=page_size, start_index=start_index, timeout=timeout, ) compressed = gzip.compress(raw_content, compresslevel=6, mtime=0) raw_path = raw_dir / f"page_{start_index:09d}.geojson.gz" atomic_write_bytes(raw_path, compressed) raw_pages.append( { "start_index": start_index, "feature_count": len(page), "response_sha256": sha256_bytes(raw_content), "response_size_bytes": len(raw_content), "artifact_path": str(raw_path.relative_to(manifest_path.parent)), "artifact_sha256": sha256_bytes(compressed), "artifact_size_bytes": len(compressed), } ) new_source_ids = 0 for raw_feature in page: source_feature_id = str(raw_feature.get("id") or "") if not source_feature_id or source_feature_id in source_seen: continue source_seen.add(source_feature_id) new_source_ids += 1 properties = dict(raw_feature.get("properties") or {}) landuse_class = str(properties.get("KLASSE") or "").strip().lower() if not definition.matches(landuse_class): continue try: source_geometry = normalize_polygonal(shape(raw_feature.get("geometry"))) except (AttributeError, TypeError, ValueError): source_geometry = None if source_geometry is None: invalid_geometry_count += 1 continue clipped = normalize_polygonal(source_geometry.intersection(boundary)) if clipped is None: outside_boundary_count += 1 continue if simplify_tolerance_degrees > 0: clipped = normalize_polygonal( clipped.simplify(simplify_tolerance_degrees, preserve_topology=True) ) if clipped is None: invalid_geometry_count += 1 continue partition_feature_id = f"{source_feature_id}:{member.nis_code}" properties.update( { "source_name": "historical_landuse", "source_feature_id": partition_feature_id, "original_source_feature_id": source_feature_id, "reference_layer_name": definition.key, "authority_level": "authoritative", "coverage_scope": scope_key, "municipality": member.name, "nis_code": member.nis_code, "observation_year": year, "historical_landuse_class": landuse_class, "attribution": ATTRIBUTION, } ) output_features.append( { "type": "Feature", "id": partition_feature_id, "geometry": mapping(clipped), "properties": properties, } ) if len(output_features) > max_features: raise RuntimeError( f"Historical {definition.key} {year} exceeds {max_features} clipped features in {member.name}" ) if len(source_seen) > max_features: raise RuntimeError( f"Historical {definition.key} {year} exceeds {max_features} source features in {member.name}" ) if len(page) < page_size: break if new_source_ids == 0: raise RuntimeError( f"Historical WFS pagination stalled for {definition.key} {year} in {member.name}" ) start_index += len(page) payload = { "type": "FeatureCollection", "name": f"{definition.label} - {member.name} {year}", "crs": GEOJSON_CRS, "features": output_features, "observation_year": year, "municipality": member.name, "nis_code": member.nis_code, "attribution": ATTRIBUTION, } atomic_write_json(output_path, payload) manifest = { "status": "complete", "year": year, "theme": definition.key, "collection": collection, "municipality": member.name, "nis_code": member.nis_code, "boundary_bbox": [float(value) for value in boundary.bounds], "boundary_geometry_sha256": boundary_geometry_sha256, "page_size": page_size, "source_feature_count": len(source_seen), "feature_count": len(output_features), "invalid_geometry_count": invalid_geometry_count, "outside_boundary_count": outside_boundary_count, "geometry_simplification_tolerance_degrees": simplify_tolerance_degrees, "raw_pages": raw_pages, "output_path": str(output_path), "output_sha256": sha256_file(output_path), "output_size_bytes": output_path.stat().st_size, "generated_at": utc_now(), } atomic_write_json(manifest_path, manifest, pretty=True) return manifest def snapshot_paths(output_root: Path, scope: GeographicScope, year: int, theme: str) -> tuple[Path, Path]: directory = output_root / "snapshots" stem = f"{scope.key.replace('-', '_')}_historical_{theme}_{year}" output_path = directory / f"{stem}.geojson" return output_path, directory / f"{stem}.manifest.json" def assemble_snapshot( *, output_root: Path, scope: GeographicScope, year: int, definition: ThemeDefinition, partitions: list[dict[str, Any]], max_total_features: int, ) -> tuple[Path, dict[str, Any]]: if len(partitions) != len(scope.members): raise RuntimeError(f"Expected {len(scope.members)} partition manifests, received {len(partitions)}") output_path, manifest_path = snapshot_paths(output_root, scope, year, definition.key) partition_identity = 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(): existing = json.loads(manifest_path.read_text(encoding="utf-8")) if ( existing.get("partition_identity_sha256") == partition_identity and sha256_file(output_path) == existing.get("output_sha256") ): return output_path, existing output_path.parent.mkdir(parents=True, exist_ok=True) temporary = output_path.with_suffix(f"{output_path.suffix}.partial") feature_count = 0 seen_ids: set[str] = set() with temporary.open("w", encoding="utf-8") as output: header = { "type": "FeatureCollection", "name": f"{definition.label} - {scope.display_name} {year}", "crs": GEOJSON_CRS, "scope_key": scope.key, "member_count": len(scope.members), "observation_year": year, "attribution": ATTRIBUTION, } output.write(json.dumps(header, ensure_ascii=False, separators=(",", ":"))[:-1]) output.write(',"features":[') first = True for partition in partitions: partition_path = Path(str(partition["output_path"])) payload = json.loads(partition_path.read_text(encoding="utf-8")) features = payload.get("features") if isinstance(payload, dict) else None if not isinstance(features, list) or len(features) != int(partition["feature_count"]): raise RuntimeError(f"Partition feature count drift for NIS {partition['nis_code']}") for feature in features: feature_id = str(feature.get("id") or "") if not feature_id or feature_id in seen_ids: raise RuntimeError(f"Duplicate or missing partition feature id in {definition.key} {year}") seen_ids.add(feature_id) if not first: output.write(",") output.write(json.dumps(feature, ensure_ascii=False, separators=(",", ":"))) first = False feature_count += 1 if feature_count > max_total_features: raise RuntimeError( f"Regional historical {definition.key} {year} exceeds {max_total_features} features" ) output.write("]}") if feature_count == 0: temporary.unlink(missing_ok=True) raise RuntimeError(f"Regional historical {definition.key} {year} contains no clipped features") temporary.replace(output_path) empty_partitions = [item["nis_code"] for item in partitions if int(item["feature_count"]) == 0] manifest = { "status": "complete", "scope_key": scope.key, "scope_type": scope.scope_type, "member_count": len(scope.members), "year": year, "theme": definition.key, "collection": COLLECTIONS[year], "feature_count": feature_count, "coverage_complete": True, "empty_partitions": empty_partitions, "partition_identity_sha256": partition_identity, "partitions": [ { "municipality": item["municipality"], "nis_code": item["nis_code"], "source_feature_count": item["source_feature_count"], "feature_count": item["feature_count"], "raw_page_count": len(item.get("raw_pages") or []), "output_sha256": item["output_sha256"], } for item in partitions ], "output_path": str(output_path), "output_sha256": sha256_file(output_path), "output_size_bytes": output_path.stat().st_size, "generated_at": utc_now(), } atomic_write_json(manifest_path, manifest, pretty=True) return output_path, manifest def locate_workspace(session: requests.Session, base_url: str, scope: GeographicScope, timeout: int): projects = response_data(session.get(f"{base_url}/api/v1/projects", params={"limit": 200}, timeout=timeout)) project = next((item for item in projects.get("items") or [] 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 = response_data( session.get(f"{base_url}/api/v1/projects/{project_id}/areas", params={"limit": 200}, timeout=timeout) ) area = next((item for item in areas.get("items") or [] if item.get("name") == scope.area_name), None) if not area: raise RuntimeError(f"Official scope Area {scope.area_name!r} is missing") datasets = response_data( session.get(f"{base_url}/api/v1/projects/{project_id}/datasets", params={"limit": 200}, timeout=timeout) ) return project_id, str(area["id"]), list(datasets.get("items") or []) def upload_snapshot( session: requests.Session, *, base_url: str, project_id: str, area_id: str, scope: GeographicScope, year: int, definition: ThemeDefinition, path: Path, manifest: dict[str, Any], simplify_tolerance_degrees: float, timeout: int, ) -> dict[str, Any]: observed_at = f"{year}-01-01T00:00:00Z" source_metadata = { "provider": "Digitaal Vlaanderen", "collection": COLLECTIONS[year], "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"]), "geometry_clipped_to_area": True, "empty_partition_nis_codes": manifest["empty_partitions"], "attribution": ATTRIBUTION, "source_catalog_url": SOURCE_CATALOG_URL, "identity_stable": False, "semantic_metrics": False, "geometry_simplification_tolerance_degrees": simplify_tolerance_degrees, "selection_aggregation": { "metric_key": f"{definition.key}_area", "method": "intersection_area", "label": f"Oppervlakte {definition.label.lower()}", "unit": "ha", "is_estimate": False, "warning": ( "Historische kaartklassen, bronkaarten en karteermethodes verschillen per bronjaar; " "interpreteer evoluties binnen die methodologische context." ), }, } provenance_metadata = { "operator_tool": "provision_regional_historical_landuse.py", "operator_explicit_fetch": True, "wfs_url": WFS_URL, "collection": COLLECTIONS[year], "partition_count": len(manifest["partitions"]), "partition_identity_sha256": manifest["partition_identity_sha256"], "combined_output_sha256": manifest["output_sha256"], "raw_source_responses_retained": True, "geometry_clipped_to_area": True, "geometry_simplification_tolerance_degrees": simplify_tolerance_degrees, "generated_at": manifest["generated_at"], } series_key = f"digitaal-vlaanderen:historical-landuse:{definition.key}:{scope.key}" 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": "historical_landuse", "reference_layer_name": definition.key, "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": series_key, "observed_at": observed_at, "valid_from": observed_at, "temporal_granularity": "year", "source_version": str(year), }, files={"file": (path.name, handle, "application/geo+json")}, timeout=timeout, ) return response_data(response) def requested_configuration(args: argparse.Namespace) -> tuple[list[int], list[ThemeDefinition]]: try: years = sorted({int(value.strip()) for value in args.years.split(",") if value.strip()}) except ValueError as exc: raise ValueError("Years must be comma-separated integers") from exc theme_keys = {value.strip().lower() for value in args.themes.split(",") if value.strip()} definitions = [item for item in THEMES if item.key in theme_keys and item.key in SUPPORTED_THEME_KEYS] unsupported_years = sorted(set(years) - set(COLLECTIONS)) unsupported_themes = sorted(theme_keys - SUPPORTED_THEME_KEYS) if not years or not definitions or unsupported_years or unsupported_themes: raise ValueError(f"Unsupported years={unsupported_years}, themes={unsupported_themes}") return years, definitions def main() -> int: args = parse_args() scope = GEOGRAPHIC_SCOPES[args.scope] try: if args.page_size <= 0 or args.max_features_per_partition <= 0 or args.max_total_features <= 0: raise ValueError("Page and feature safety limits must be greater than zero") years, definitions = requested_configuration(args) members_path = resolve_member_boundaries(scope, args.scope_output_root) boundaries = load_member_boundaries(members_path, scope) output_root = args.output_root / scope.key prepared: list[tuple[int, ThemeDefinition, Path, dict[str, Any]]] = [] with build_session() as source_session: for year in years: for definition in definitions: partitions = [ prepare_partition( source_session, output_root=output_root, year=year, definition=definition, scope_key=scope.key, member=member, boundary=boundaries[member.nis_code][1], page_size=args.page_size, max_features=args.max_features_per_partition, simplify_tolerance_degrees=args.simplify_tolerance_degrees, timeout=args.request_timeout, force=args.force, ) for member in scope.members ] path, manifest = assemble_snapshot( output_root=output_root, scope=scope, year=year, definition=definition, partitions=partitions, max_total_features=args.max_total_features, ) prepared.append((year, definition, path, manifest)) results: list[dict[str, Any]] = [] if args.fetch_only: results = [ { "year": year, "theme": definition.key, "feature_count": manifest["feature_count"], "partition_count": len(manifest["partitions"]), "empty_partitions": manifest["empty_partitions"], "path": str(path), "status": "prepared", } for year, definition, path, manifest in prepared ] else: base_url = args.base_url.rstrip("/") with requests.Session() as api_session: project_id, area_id, existing = locate_workspace(api_session, base_url, scope, args.import_timeout) for year, definition, path, manifest in prepared: series_key = f"digitaal-vlaanderen:historical-landuse:{definition.key}:{scope.key}" dataset = next( ( item for item in existing if item.get("temporal_series_key") == series_key and str(item.get("observed_at") or "").startswith(str(year)) ), None, ) if dataset: if int(dataset.get("feature_count") or 0) != int(manifest["feature_count"]): raise RuntimeError( f"Persisted {definition.key} {year} count differs from the audited source artifact" ) status = "existing" else: dataset = upload_snapshot( api_session, base_url=base_url, project_id=project_id, area_id=area_id, scope=scope, year=year, definition=definition, path=path, manifest=manifest, simplify_tolerance_degrees=args.simplify_tolerance_degrees, timeout=args.import_timeout, ) existing.append(dataset) status = "imported" results.append( { "year": year, "theme": definition.key, "dataset_id": dataset["id"], "feature_count": dataset.get("feature_count"), "partition_count": len(manifest["partitions"]), "empty_partitions": manifest["empty_partitions"], "status": status, } ) except (OSError, RuntimeError, ValueError, KeyError, requests.RequestException, json.JSONDecodeError) as exc: print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr) return 1 print( json.dumps( { "status": "ok", "mode": "fetch_only" if args.fetch_only else "synchronized", "scope": scope.key, "display_name": scope.display_name, "member_count": len(scope.members), "snapshots": results, }, ensure_ascii=False, indent=2, ) ) return 0 if __name__ == "__main__": sys.exit(main())