Add regional historical land-use operator

scripts/provision_regional_historical_landuse.py

@@ -0,0 +1,756 @@
+"""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"]),
+        "empty_partition_nis_codes": manifest["empty_partitions"],
+        "attribution": ATTRIBUTION,
+        "source_catalog_url": SOURCE_CATALOG_URL,
+        "identity_stable": False,
+        "geometry_simplification_tolerance_degrees": simplify_tolerance_degrees,
+        "selection_aggregation": {
+            "method": "intersection_area",
+            "label": "Oppervlakte",
+            "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_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())