geointel/scripts/provision_official_landuse_timeseries.py

"""Provision official modern Flemish land-use snapshots for an explicit area.

The operator downloads categorical 10 metre GeoTIFF subsets from the public
Departement Omgeving MercatorNet WCS, clips them to a supplied boundary and
polygonizes only explicitly supported classes. Raw rasters and checksum
manifests remain provenance artifacts. Vector output is imported through the
normal GeoIntel dataset API and is never written directly to PostGIS.

Mol is the safe default. Other scopes must provide their own boundary, project,
area name and identity explicitly.
"""

from __future__ import annotations

import argparse
import hashlib
import json
import math
import os
import sys
from dataclasses import dataclass
from datetime import datetime, timezone
from pathlib import Path
from typing import Any

import numpy as np
import rasterio
import requests
from pyproj import Transformer
from rasterio.features import geometry_mask, shapes
from rasterio.merge import merge
from requests.adapters import HTTPAdapter
from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, mapping, shape
from shapely.ops import transform, unary_union
from shapely.validation import make_valid
from urllib3.util.retry import Retry


DEFAULT_API_URL = "http://127.0.0.1:8000"
DEFAULT_PROJECT_NAME = "Mol Municipality Workbench"
DEFAULT_AREA_NAME = "Gemeente Mol"
DEFAULT_MUNICIPALITY_NAME = "Mol"
DEFAULT_NIS_CODE = "13025"
DEFAULT_SCOPE_KEY = "mol"
DEFAULT_BOUNDARY_PATH = Path("/app/storage/operator-data/mol-municipality/mol_municipality_boundary.geojson")
DEFAULT_OUTPUT_DIR = Path("/app/storage/operator-data/official-landuse/mol")

WCS_URL = "https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs"
WCS_VERSION = "1.0.0"
SOURCE_CRS = "EPSG:31370"
OUTPUT_CRS = "EPSG:4326"
SOURCE_RESOLUTION_METRES = 10.0
SUPPORTED_YEARS = (2013, 2016, 2019, 2022, 2025)
ATTRIBUTION = "Bron: Landgebruik Vlaanderen, Departement Omgeving"
SERIES_LABEL = "Moderne landgebruikskaart (10 m)"

CATALOGUE_URLS = {
    year: f"https://www.vlaanderen.be/datavindplaats/catalogus/landgebruik-vlaanderen-toestand-{year}"
    for year in SUPPORTED_YEARS
}

LAND_USE_CLASSES = {
    1: "Huizen en tuinen",
    2: "Industrie",
    3: "Commerciele doeleinden",
    4: "Diensten",
    5: "Transportinfrastructuur",
    6: "Recreatie",
    7: "Landbouwgebouwen en -infrastructuur",
    8: "Overige bebouwde terreinen",
    9: "Overige onbebouwde terreinen",
    10: "Actieve groeves",
    11: "Luchthavens",
    12: "Bos",
    13: "Akker",
    14: "Grasland in landbouwgebruik",
    15: "Struikgewas",
    16: "Braakliggend en duinen",
    17: "Water",
    18: "Moeras",
    19: "Overige graslanden",
}


@dataclass(frozen=True)
class ThemeDefinition:
    key: str
    label: str
    class_ids: tuple[int, ...]
    reference_layer_name: str
    metric_label: str


@dataclass(frozen=True)
class PreparedSnapshot:
    year: int
    theme: ThemeDefinition
    raster_path: Path
    vector_path: Path
    manifest_path: Path
    feature_count: int
    raster_sha256: str
    vector_sha256: str


THEMES = (
    ThemeDefinition("forest", "Bos", (12,), "forest", "Bosoppervlakte"),
    ThemeDefinition("water", "Water", (17,), "water", "Wateroppervlakte"),
    ThemeDefinition(
        "built",
        "Bebouwde functies",
        (1, 2, 3, 4, 6, 7, 8),
        "buildings",
        "Oppervlakte bebouwde functies",
    ),
    ThemeDefinition("transport", "Transportinfrastructuur", (5,), "roads", "Oppervlakte transportinfrastructuur"),
)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description="Provision official 2013-2025 Flemish land-use snapshots.")
    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL))
    parser.add_argument("--project-name", default=DEFAULT_PROJECT_NAME)
    parser.add_argument("--area-name", default=DEFAULT_AREA_NAME, help="Case-insensitive fragment identifying the persisted Area.")
    parser.add_argument("--municipality-name", default=DEFAULT_MUNICIPALITY_NAME)
    parser.add_argument("--nis-code", default=DEFAULT_NIS_CODE)
    parser.add_argument("--scope-key", default=DEFAULT_SCOPE_KEY)
    parser.add_argument("--years", default=",".join(str(year) for year in SUPPORTED_YEARS))
    parser.add_argument("--themes", default="forest,water,built,transport")
    parser.add_argument(
        "--boundary-path",
        type=Path,
        default=Path(os.environ.get("OFFICIAL_LANDUSE_BOUNDARY_PATH", DEFAULT_BOUNDARY_PATH)),
    )
    parser.add_argument(
        "--partition-boundaries-path",
        type=Path,
        default=None,
        help="Optional FeatureCollection of non-overlapping member boundaries used for bounded WCS requests.",
    )
    parser.add_argument(
        "--output-dir",
        type=Path,
        default=Path(os.environ.get("OFFICIAL_LANDUSE_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)),
    )
    parser.add_argument("--request-timeout", type=int, default=240)
    parser.add_argument("--import-timeout", type=int, default=1800)
    parser.add_argument("--max-features", type=int, default=100000)
    parser.add_argument("--force", action="store_true", help="Refetch and rebuild local source artifacts; persisted datasets stay immutable.")
    parser.add_argument("--fetch-only", action="store_true", help="Prepare and verify artifacts without changing GeoIntel persistence.")
    return parser.parse_args()


def utc_now() -> str:
    return datetime.now(timezone.utc).isoformat()


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 write_json_atomic(path: Path, payload: dict[str, Any], *, pretty: bool = False) -> None:
    temporary = path.with_suffix(f"{path.suffix}.partial")
    temporary.write_text(
        json.dumps(
            payload,
            ensure_ascii=False,
            indent=2 if pretty else None,
            separators=None if pretty else (",", ":"),
            sort_keys=pretty,
        ),
        encoding="utf-8",
    )
    temporary.replace(path)


def build_session() -> requests.Session:
    retry = Retry(
        total=5,
        connect=5,
        read=5,
        status=5,
        backoff_factor=1.0,
        status_forcelist=(429, 500, 502, 503, 504),
        allowed_methods=frozenset({"GET"}),
        raise_on_status=True,
    )
    session = requests.Session()
    session.headers.update({"User-Agent": "GeoIntel-Official-Landuse-Operator/1.0"})
    adapter = HTTPAdapter(max_retries=retry)
    session.mount("https://", adapter)
    session.mount("http://", adapter)
    return session


def coverage_id(year: int) -> str:
    return f"lu:lu_landgebruik_vlaa_{year}_v3"


def series_key(theme: ThemeDefinition, scope_key: str) -> str:
    return f"department-omgeving:land-use:{theme.key}:{scope_key.strip().lower()}"


def scope_identity(scope_display_name: str, raw_nis_codes: str) -> dict[str, Any]:
    nis_codes = [value.strip() for value in raw_nis_codes.split(",") if value.strip()]
    return {
        "scope_display_name": scope_display_name,
        "member_nis_codes": nis_codes,
        "municipality": scope_display_name if len(nis_codes) == 1 else None,
        "nis_code": nis_codes[0] if len(nis_codes) == 1 else None,
    }


def load_boundary(path: Path):
    if not path.exists():
        raise RuntimeError(f"Boundary artifact is missing at {path}")
    payload = json.loads(path.read_text(encoding="utf-8"))
    features = payload.get("features") or []
    if len(features) != 1:
        raise RuntimeError("Boundary artifact must contain exactly one feature")
    boundary = normalize_polygonal(shape(features[0].get("geometry")))
    if boundary is None:
        raise RuntimeError("Boundary artifact is empty, invalid or non-polygonal")
    return boundary


def load_partition_boundaries(path: Path | None) -> list[tuple[str, Any]]:
    if path is None:
        return []
    if not path.is_file():
        raise RuntimeError(f"Partition boundary artifact is missing at {path}")
    payload = json.loads(path.read_text(encoding="utf-8"))
    features = payload.get("features") or []
    if not features:
        raise RuntimeError("Partition boundary artifact contains no features")
    partitions: list[tuple[str, Any]] = []
    seen: set[str] = set()
    for index, feature in enumerate(features):
        properties = feature.get("properties") or {}
        key = str(properties.get("nis_code") or properties.get("NISCODE") or feature.get("id") or index).strip()
        if not key or key in seen:
            raise RuntimeError(f"Partition boundary artifact contains an empty or duplicate key: {key!r}")
        geometry = normalize_polygonal(shape(feature.get("geometry")))
        if geometry is None:
            raise RuntimeError(f"Partition boundary {key} is empty, invalid or non-polygonal")
        seen.add(key)
        partitions.append((key, geometry))
    return partitions


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 not polygons:
            return None
        merged = unary_union(polygons)
        return merged if isinstance(merged, (Polygon, MultiPolygon)) and not merged.is_empty else None
    return None


def metric_boundary(boundary):
    transformer = Transformer.from_crs(OUTPUT_CRS, SOURCE_CRS, always_xy=True)
    projected = normalize_polygonal(transform(transformer.transform, boundary))
    if projected is None:
        raise RuntimeError("Boundary could not be projected to EPSG:31370")
    return projected


def snapped_bounds(bounds: tuple[float, float, float, float]) -> tuple[float, float, float, float]:
    min_x, min_y, max_x, max_y = bounds
    resolution = SOURCE_RESOLUTION_METRES
    return (
        math.floor(min_x / resolution) * resolution,
        math.floor(min_y / resolution) * resolution,
        math.ceil(max_x / resolution) * resolution,
        math.ceil(max_y / resolution) * resolution,
    )


def build_wcs_params(year: int, bounds: tuple[float, float, float, float]) -> dict[str, str]:
    min_x, min_y, max_x, max_y = snapped_bounds(bounds)
    return {
        "SERVICE": "WCS",
        "VERSION": WCS_VERSION,
        "REQUEST": "GetCoverage",
        "COVERAGE": coverage_id(year),
        "CRS": SOURCE_CRS,
        "BBOX": f"{min_x:.3f},{min_y:.3f},{max_x:.3f},{max_y:.3f}",
        "RESX": str(int(SOURCE_RESOLUTION_METRES)),
        "RESY": str(int(SOURCE_RESOLUTION_METRES)),
        "FORMAT": "image/tiff",
        "RESPONSE_CRS": SOURCE_CRS,
    }


def verify_coverages(session: requests.Session, years: list[int], timeout: int) -> None:
    response = session.get(
        WCS_URL,
        params={"SERVICE": "WCS", "VERSION": WCS_VERSION, "REQUEST": "GetCapabilities"},
        timeout=timeout,
    )
    response.raise_for_status()
    missing = [coverage_id(year) for year in years if coverage_id(year) not in response.text]
    if missing:
        raise RuntimeError(f"Official WCS is missing expected coverages: {', '.join(missing)}")


def validate_raster(path: Path) -> dict[str, Any]:
    try:
        with rasterio.open(path) as dataset:
            epsg = dataset.crs.to_epsg() if dataset.crs else None
            resolution = (abs(float(dataset.res[0])), abs(float(dataset.res[1])))
            if epsg != 31370:
                raise RuntimeError(f"Expected EPSG:31370 source raster, received {dataset.crs}")
            if dataset.count != 1:
                raise RuntimeError(f"Expected one categorical raster band, received {dataset.count}")
            if any(abs(value - SOURCE_RESOLUTION_METRES) > 0.01 for value in resolution):
                raise RuntimeError(f"Expected 10 metre source resolution, received {resolution}")
            if not np.issubdtype(np.dtype(dataset.dtypes[0]), np.integer):
                raise RuntimeError(f"Expected integer land-use classes, received {dataset.dtypes[0]}")
            return {
                "width": dataset.width,
                "height": dataset.height,
                "dtype": dataset.dtypes[0],
                "nodata": dataset.nodata,
                "crs": SOURCE_CRS,
                "resolution_metres": SOURCE_RESOLUTION_METRES,
                "bounds": list(dataset.bounds),
            }
    except rasterio.errors.RasterioError as exc:
        raise RuntimeError(f"Official WCS response is not a readable GeoTIFF: {exc}") from exc


def download_raster(
    session: requests.Session,
    *,
    year: int,
    bounds: tuple[float, float, float, float],
    path: Path,
    timeout: int,
) -> dict[str, Any]:
    response = session.get(WCS_URL, params=build_wcs_params(year, bounds), timeout=timeout, stream=True)
    response.raise_for_status()
    content_type = str(response.headers.get("content-type") or "").lower()
    if "tiff" not in content_type:
        preview = response.content[:500].decode("utf-8", errors="replace")
        raise RuntimeError(f"Official WCS returned {content_type or 'unknown content'} instead of GeoTIFF: {preview}")
    temporary = path.with_suffix(f"{path.suffix}.partial")
    try:
        with temporary.open("wb") as handle:
            for chunk in response.iter_content(chunk_size=1024 * 1024):
                if chunk:
                    handle.write(chunk)
        profile = validate_raster(temporary)
        temporary.replace(path)
    except Exception:
        temporary.unlink(missing_ok=True)
        raise
    return {**profile, "request_url": response.url, "retrieved_at": utc_now()}


def merge_partition_rasters(paths: list[Path], destination: Path) -> dict[str, Any]:
    if not paths:
        raise RuntimeError("At least one partition raster is required for a regional mosaic")
    sources = [rasterio.open(path) for path in paths]
    temporary = destination.with_suffix(f"{destination.suffix}.partial")
    try:
        mosaic, transform_matrix = merge(sources)
        profile = sources[0].profile.copy()
        profile.pop("blockxsize", None)
        profile.pop("blockysize", None)
        tile_width = min(512, (mosaic.shape[2] // 16) * 16)
        tile_height = min(512, (mosaic.shape[1] // 16) * 16)
        profile.update(
            driver="GTiff",
            width=mosaic.shape[2],
            height=mosaic.shape[1],
            transform=transform_matrix,
            count=mosaic.shape[0],
            compress="deflate",
            tiled=tile_width >= 16 and tile_height >= 16,
            BIGTIFF="IF_SAFER",
        )
        if tile_width >= 16 and tile_height >= 16:
            profile.update(blockxsize=tile_width, blockysize=tile_height)
        with rasterio.open(temporary, "w", **profile) as target:
            target.write(mosaic)
        validated = validate_raster(temporary)
        temporary.replace(destination)
        return validated
    except Exception:
        temporary.unlink(missing_ok=True)
        raise
    finally:
        for source in sources:
            source.close()


def download_partitioned_raster(
    session: requests.Session,
    *,
    year: int,
    partitions: list[tuple[str, Any]],
    output_dir: Path,
    destination: Path,
    timeout: int,
    force: bool,
) -> dict[str, Any]:
    partition_dir = output_dir / "source-partitions" / str(year)
    partition_dir.mkdir(parents=True, exist_ok=True)
    paths: list[Path] = []
    details: list[dict[str, Any]] = []
    for key, boundary in partitions:
        path = partition_dir / f"land_use_{year}_{key}.tif"
        boundary_metric = metric_boundary(boundary)
        if force or not path.is_file():
            profile = download_raster(
                session,
                year=year,
                bounds=boundary_metric.bounds,
                path=path,
                timeout=timeout,
            )
            status = "downloaded"
        else:
            profile = validate_raster(path)
            profile["request_url"] = requests.Request(
                "GET", WCS_URL, params=build_wcs_params(year, boundary_metric.bounds)
            ).prepare().url
            profile["retrieved_at"] = None
            status = "reused"
        paths.append(path)
        details.append(
            {
                "partition_key": key,
                "path": str(path),
                "sha256": sha256_file(path),
                "status": status,
                "request_url": profile.get("request_url"),
                "width": profile["width"],
                "height": profile["height"],
            }
        )

    mosaic_profile = merge_partition_rasters(paths, destination)
    return {
        **mosaic_profile,
        "request_mode": "partitioned_scope_members",
        "partition_count": len(details),
        "partitions": details,
        "retrieved_at": utc_now(),
    }


def polygonize_snapshot(
    *,
    raster_path: Path,
    boundary,
    year: int,
    theme: ThemeDefinition,
    municipality_name: str,
    nis_code: str,
    scope_key: str,
    max_features: int,
) -> tuple[dict[str, Any], dict[str, Any]]:
    identity = scope_identity(municipality_name, nis_code)
    boundary_metric = metric_boundary(boundary)
    to_output = Transformer.from_crs(SOURCE_CRS, OUTPUT_CRS, always_xy=True)
    raster_sha256 = sha256_file(raster_path)

    with rasterio.open(raster_path) as dataset:
        validate_raster(raster_path)
        values = dataset.read(1)
        inside = geometry_mask(
            [mapping(boundary_metric)],
            out_shape=values.shape,
            transform=dataset.transform,
            invert=True,
            # Exact clipping happens after polygonization, so retain every
            # classified source cell that overlaps the requested boundary.
            all_touched=True,
        )
        valid_inside = inside.copy()
        nodata_count = 0
        if dataset.nodata is not None:
            nodata = np.isclose(values, dataset.nodata)
            nodata_count = int(np.count_nonzero(inside & nodata))
            valid_inside &= ~nodata
        available_values, available_counts = np.unique(values[valid_inside], return_counts=True)
        class_histogram = {str(int(value)): int(count) for value, count in zip(available_values, available_counts)}
        unknown_classes = sorted(int(value) for value in available_values if int(value) not in LAND_USE_CLASSES)
        if unknown_classes:
            raise RuntimeError(f"Land-use raster {year} contains undocumented classes: {unknown_classes}")

        class_mask = np.isin(values, np.asarray(theme.class_ids)) & valid_inside
        source_pixel_count = int(np.count_nonzero(class_mask))
        if source_pixel_count == 0:
            raise RuntimeError(f"Land-use raster {year} contains no {theme.label} cells inside the boundary")

        features: list[dict[str, Any]] = []
        polygon_area_m2 = 0.0
        for raw_geometry, raw_value in shapes(
            class_mask.astype("uint8"),
            mask=class_mask,
            transform=dataset.transform,
            connectivity=4,
        ):
            if int(raw_value) != 1:
                continue
            geometry_metric = normalize_polygonal(shape(raw_geometry).intersection(boundary_metric))
            if geometry_metric is None or geometry_metric.area <= 0:
                continue
            geometry_output = normalize_polygonal(transform(to_output.transform, geometry_metric))
            if geometry_output is None:
                continue
            if len(features) >= max_features:
                raise RuntimeError(
                    f"Land-use {theme.key} {year} exceeds the {max_features} feature safety limit; refusing truncation"
                )
            feature_hash = hashlib.sha256(
                f"{year}:{theme.key}:".encode("utf-8") + geometry_metric.wkb
            ).hexdigest()[:24]
            feature_id = f"land-use-{year}-{theme.key}-{feature_hash}"
            area_m2 = float(geometry_metric.area)
            polygon_area_m2 += area_m2
            features.append(
                {
                    "type": "Feature",
                    "id": feature_id,
                    "geometry": mapping(geometry_output),
                    "properties": {
                        "source_name": "department_omgeving_land_use",
                        "source_feature_id": feature_id,
                        "reference_layer_name": theme.reference_layer_name,
                        "layer_type": theme.reference_layer_name,
                        "authority_level": "authoritative",
                        "coverage_scope": scope_key,
                        **identity,
                        "observation_year": year,
                        "land_use_class_ids": list(theme.class_ids),
                        "land_use_class_names": [LAND_USE_CLASSES[class_id] for class_id in theme.class_ids],
                        "source_resolution_m": SOURCE_RESOLUTION_METRES,
                        "polygon_area_m2": round(area_m2, 3),
                        "source_coverage_id": coverage_id(year),
                        "source_raster_sha256": raster_sha256,
                        "attribution": ATTRIBUTION,
                    },
                }
            )

    payload = {
        "type": "FeatureCollection",
        "name": f"{theme.label} - {municipality_name} {year}",
        "crs": {"type": "name", "properties": {"name": OUTPUT_CRS}},
        **identity,
        "scope_key": scope_key,
        "observation_year": year,
        "source_coverage_id": coverage_id(year),
        "source_resolution_m": SOURCE_RESOLUTION_METRES,
        "attribution": ATTRIBUTION,
        "features": features,
    }
    stats = {
        "feature_count": len(features),
        "source_pixel_count": source_pixel_count,
        "source_pixel_area_m2": source_pixel_count * SOURCE_RESOLUTION_METRES**2,
        "polygon_area_m2": polygon_area_m2,
        "nodata_pixels_inside_boundary": nodata_count,
        "class_histogram": class_histogram,
        "raster_sha256": raster_sha256,
    }
    return payload, stats


def existing_snapshot(
    *,
    year: int,
    theme: ThemeDefinition,
    raster_path: Path,
    vector_path: Path,
    manifest_path: Path,
) -> PreparedSnapshot | None:
    if not raster_path.exists() or not vector_path.exists() or not manifest_path.exists():
        return None
    try:
        manifest = json.loads(manifest_path.read_text(encoding="utf-8"))
        raster_sha256 = sha256_file(raster_path)
        vector_sha256 = sha256_file(vector_path)
        if (
            manifest.get("year") != year
            or manifest.get("theme") != theme.key
            or manifest.get("coverage_id") != coverage_id(year)
            or manifest.get("class_ids") != list(theme.class_ids)
            or manifest.get("raster_sha256") != raster_sha256
            or manifest.get("vector_sha256") != vector_sha256
        ):
            return None
        validate_raster(raster_path)
        feature_count = int(manifest["feature_count"])
        if feature_count <= 0:
            return None
        return PreparedSnapshot(
            year=year,
            theme=theme,
            raster_path=raster_path,
            vector_path=vector_path,
            manifest_path=manifest_path,
            feature_count=feature_count,
            raster_sha256=raster_sha256,
            vector_sha256=vector_sha256,
        )
    except (OSError, ValueError, KeyError, RuntimeError, rasterio.errors.RasterioError):
        return None


def prepare_snapshot(
    session: requests.Session,
    *,
    args: argparse.Namespace,
    boundary,
    boundary_metric,
    partition_boundaries: list[tuple[str, Any]],
    year: int,
    theme: ThemeDefinition,
    refresh_source: bool = False,
) -> PreparedSnapshot:
    stem = f"{args.scope_key}_land_use_{theme.key}_{year}"
    legacy_forest_raster = args.output_dir / f"{args.scope_key}_land_use_forest_{year}.tif"
    shared_raster = args.output_dir / f"{args.scope_key}_land_use_source_{year}.tif"
    raster_path = legacy_forest_raster if not args.force and legacy_forest_raster.exists() else shared_raster
    vector_path = args.output_dir / f"{stem}.geojson"
    manifest_path = args.output_dir / f"{stem}.manifest.json"
    if not args.force:
        prepared = existing_snapshot(
            year=year,
            theme=theme,
            raster_path=raster_path,
            vector_path=vector_path,
            manifest_path=manifest_path,
        )
        if prepared:
            return prepared

    raster_profile: dict[str, Any]
    if refresh_source or not raster_path.exists():
        if partition_boundaries:
            raster_profile = download_partitioned_raster(
                session,
                year=year,
                partitions=partition_boundaries,
                output_dir=args.output_dir,
                destination=raster_path,
                timeout=args.request_timeout,
                force=args.force,
            )
        else:
            raster_profile = download_raster(
                session,
                year=year,
                bounds=boundary_metric.bounds,
                path=raster_path,
                timeout=args.request_timeout,
            )
    else:
        raster_profile = validate_raster(raster_path)
        raster_profile["request_url"] = requests.Request(
            "GET", WCS_URL, params=build_wcs_params(year, boundary_metric.bounds)
        ).prepare().url
        raster_profile["retrieved_at"] = None

    payload, stats = polygonize_snapshot(
        raster_path=raster_path,
        boundary=boundary,
        year=year,
        theme=theme,
        municipality_name=args.municipality_name,
        nis_code=args.nis_code,
        scope_key=args.scope_key,
        max_features=args.max_features,
    )
    write_json_atomic(vector_path, payload)
    vector_sha256 = sha256_file(vector_path)
    identity = scope_identity(args.municipality_name, args.nis_code)
    manifest = {
        "schema_version": 1,
        "year": year,
        "theme": theme.key,
        "class_ids": list(theme.class_ids),
        "class_names": [LAND_USE_CLASSES[class_id] for class_id in theme.class_ids],
        "coverage_id": coverage_id(year),
        "catalogue_url": CATALOGUE_URLS[year],
        "wcs_url": WCS_URL,
        "wcs_version": WCS_VERSION,
        "wcs_request_url": raster_profile.get("request_url"),
        "source_crs": SOURCE_CRS,
        "output_crs": OUTPUT_CRS,
        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
        "boundary_path": str(args.boundary_path),
        "partition_boundaries_path": (
            str(getattr(args, "partition_boundaries_path", None))
            if getattr(args, "partition_boundaries_path", None)
            else None
        ),
        **identity,
        "scope_key": args.scope_key,
        "raster_path": str(raster_path),
        "vector_path": str(vector_path),
        "raster_profile": raster_profile,
        "raster_sha256": stats["raster_sha256"],
        "vector_sha256": vector_sha256,
        "feature_count": stats["feature_count"],
        "source_pixel_count": stats["source_pixel_count"],
        "source_pixel_area_m2": stats["source_pixel_area_m2"],
        "polygon_area_m2": stats["polygon_area_m2"],
        "nodata_pixels_inside_boundary": stats["nodata_pixels_inside_boundary"],
        "class_histogram": stats["class_histogram"],
        "generated_at": utc_now(),
    }
    write_json_atomic(manifest_path, manifest, pretty=True)
    return PreparedSnapshot(
        year=year,
        theme=theme,
        raster_path=raster_path,
        vector_path=vector_path,
        manifest_path=manifest_path,
        feature_count=int(stats["feature_count"]),
        raster_sha256=str(stats["raster_sha256"]),
        vector_sha256=vector_sha256,
    )


def response_data(response: requests.Response) -> Any:
    try:
        payload = response.json()
    except ValueError as exc:
        raise RuntimeError(f"GeoIntel API returned non-JSON ({response.status_code}): {response.text[:300]}") from exc
    if not response.ok:
        raise RuntimeError(f"GeoIntel API failed ({response.status_code}): {json.dumps(payload, ensure_ascii=False)[:800]}")
    if not isinstance(payload, dict) or "data" not in payload:
        raise RuntimeError("GeoIntel API response does not use the canonical data envelope")
    return payload["data"]


def list_paginated_items(session: requests.Session, url: str, *, timeout: int) -> list[dict[str, Any]]:
    items: list[dict[str, Any]] = []
    offset = 0
    total: int | None = None
    while total is None or offset < total:
        page = response_data(session.get(url, params={"limit": 200, "offset": offset}, timeout=timeout))
        page_items = page.get("items") if isinstance(page, dict) else None
        if not isinstance(page_items, list):
            raise RuntimeError(f"GeoIntel list response for {url} has no items array")
        if total is None:
            total = int(page.get("total", len(page_items)))
        items.extend(page_items)
        if not page_items:
            break
        offset += len(page_items)
    if total is not None and len(items) != total:
        raise RuntimeError(f"GeoIntel list response for {url} returned {len(items)} of {total} items")
    return items


def locate_workspace(session: requests.Session, base_url: str, args: argparse.Namespace):
    projects = list_paginated_items(session, f"{base_url}/api/v1/projects", timeout=args.import_timeout)
    project = next((item for item in projects if item.get("name") == args.project_name), None)
    if not project:
        raise RuntimeError(f"Project {args.project_name!r} is missing")
    project_id = str(project["id"])
    areas = list_paginated_items(
        session,
        f"{base_url}/api/v1/projects/{project_id}/areas",
        timeout=args.import_timeout,
    )
    area_fragment = args.area_name.strip().casefold()
    matches = [item for item in areas if area_fragment in str(item.get("name") or "").casefold()]
    if len(matches) != 1:
        raise RuntimeError(f"Expected one Area matching {args.area_name!r}, received {len(matches)}")
    datasets = list_paginated_items(
        session,
        f"{base_url}/api/v1/projects/{project_id}/datasets",
        timeout=args.import_timeout,
    )
    return project_id, str(matches[0]["id"]), datasets


def build_source_metadata(args: argparse.Namespace, snapshot: PreparedSnapshot) -> dict[str, Any]:
    identity = scope_identity(args.municipality_name, args.nis_code)
    return {
        "provider": "Departement Omgeving",
        "source_title": f"Landgebruik - Vlaanderen - toestand {snapshot.year}",
        "catalogue_url": CATALOGUE_URLS[snapshot.year],
        "coverage_id": coverage_id(snapshot.year),
        "authority_level": "authoritative",
        "coverage_scope": args.scope_key,
        **identity,
        "geometry_clipped_to_area": True,
        "attribution": ATTRIBUTION,
        "license_note": "Publieke Vlaamse overheidsdata; raadpleeg de toegangs- en gebruiksvoorwaarden in de bronmetadata.",
        "methodology_version": "3",
        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
        "source_crs": SOURCE_CRS,
        "polygon_crs": OUTPUT_CRS,
        "land_use_class_ids": list(snapshot.theme.class_ids),
        "land_use_class_names": [LAND_USE_CLASSES[class_id] for class_id in snapshot.theme.class_ids],
        "theme": snapshot.theme.reference_layer_name,
        "temporal_series_label": (
            SERIES_LABEL if snapshot.theme.key == "forest" else f"{SERIES_LABEL} - {snapshot.theme.label}"
        ),
        "observation_date_precision": "year",
        "identity_stable": False,
        "semantic_metrics": False,
        "identity_limitation": "Raster-derived polygons can split or merge between source editions; object lineage is not inferred.",
        "selection_aggregation": {
            "method": "intersection_area",
            "metric_key": f"{snapshot.theme.key}_area",
            "label": snapshot.theme.metric_label,
            "unit": "ha",
            "is_estimate": False,
            "warning": "Oppervlakte is exact binnen de officiele 10 m rasterrepresentatie en is niet perceelsnauwkeurig.",
        },
        "comparison_limitation": "Compare editions as 10 m land-use states; source inputs and methodology can evolve between publication years.",
    }


def build_provenance_metadata(args: argparse.Namespace, snapshot: PreparedSnapshot) -> dict[str, Any]:
    return {
        "operator_tool": "provision_official_landuse_timeseries.py",
        "operator_explicit_fetch": True,
        "geometry_clipped_to_area": True,
        "wcs_url": WCS_URL,
        "wcs_version": WCS_VERSION,
        "coverage_id": coverage_id(snapshot.year),
        "catalogue_url": CATALOGUE_URLS[snapshot.year],
        "raw_raster_path": str(snapshot.raster_path),
        "polygon_artifact_path": str(snapshot.vector_path),
        "manifest_path": str(snapshot.manifest_path),
        "partition_boundaries_path": (
            str(getattr(args, "partition_boundaries_path", None))
            if getattr(args, "partition_boundaries_path", None)
            else None
        ),
        "raster_sha256": snapshot.raster_sha256,
        "vector_sha256": snapshot.vector_sha256,
        "source_crs": SOURCE_CRS,
        "output_crs": OUTPUT_CRS,
        "source_resolution_metres": SOURCE_RESOLUTION_METRES,
        "generated_at": utc_now(),
    }


def upload_snapshot(
    session: requests.Session,
    *,
    base_url: str,
    project_id: str,
    area_id: str,
    args: argparse.Namespace,
    snapshot: PreparedSnapshot,
) -> dict[str, Any]:
    observed_at = f"{snapshot.year}-01-01T00:00:00Z"
    with snapshot.vector_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": "department_omgeving_land_use",
                "reference_layer_name": snapshot.theme.reference_layer_name,
                "source_metadata_json": json.dumps(build_source_metadata(args, snapshot), ensure_ascii=False),
                "provenance_metadata_json": json.dumps(build_provenance_metadata(args, snapshot), ensure_ascii=False),
                "area_id": area_id,
                "temporal_series_key": series_key(snapshot.theme, args.scope_key),
                "observed_at": observed_at,
                "valid_from": observed_at,
                "temporal_granularity": "year",
                "source_version": f"{snapshot.year}-v3",
            },
            files={"file": (snapshot.vector_path.name, handle, "application/geo+json")},
            timeout=args.import_timeout,
        )
    return response_data(response)


def main() -> int:
    args = parse_args()
    try:
        years = sorted({int(value.strip()) for value in args.years.split(",") if value.strip()})
    except ValueError:
        print(json.dumps({"status": "error", "message": "Years must be comma-separated integers"}), file=sys.stderr)
        return 2
    requested_themes = {value.strip().lower() for value in args.themes.split(",") if value.strip()}
    definitions = [definition for definition in THEMES if definition.key in requested_themes]
    unsupported_years = [year for year in years if year not in SUPPORTED_YEARS]
    unsupported_themes = requested_themes - {definition.key for definition in THEMES}
    if unsupported_years or unsupported_themes or not years or not definitions:
        print(
            json.dumps(
                {"status": "error", "message": f"Unsupported years={unsupported_years}, themes={sorted(unsupported_themes)}"}
            ),
            file=sys.stderr,
        )
        return 2
    if not args.scope_key.strip() or not args.project_name.strip() or not args.area_name.strip():
        print(json.dumps({"status": "error", "message": "scope-key, project-name and area-name are required"}), file=sys.stderr)
        return 2

    args.output_dir.mkdir(parents=True, exist_ok=True)
    results: list[dict[str, Any]] = []
    try:
        boundary = load_boundary(args.boundary_path)
        boundary_metric = metric_boundary(boundary)
        partition_boundaries = load_partition_boundaries(args.partition_boundaries_path)
        prepared: list[PreparedSnapshot] = []
        with build_session() as source_session:
            verify_coverages(source_session, years, args.request_timeout)
            for year in years:
                for theme in definitions:
                    prepared.append(
                        prepare_snapshot(
                            source_session,
                            args=args,
                            boundary=boundary,
                            boundary_metric=boundary_metric,
                            partition_boundaries=partition_boundaries,
                            year=year,
                            theme=theme,
                            refresh_source=args.force and theme is definitions[0],
                        )
                    )

        if args.fetch_only:
            results = [
                {
                    "year": item.year,
                    "theme": item.theme.key,
                    "status": "prepared",
                    "feature_count": item.feature_count,
                    "raster_path": str(item.raster_path),
                    "vector_path": str(item.vector_path),
                    "manifest_path": str(item.manifest_path),
                }
                for item 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, args)
                for item in prepared:
                    key = series_key(item.theme, args.scope_key)
                    observed_date = f"{item.year}-01-01"
                    dataset = next(
                        (
                            candidate
                            for candidate in existing
                            if candidate.get("temporal_series_key") == key
                            and str(candidate.get("observed_at") or "").startswith(observed_date)
                        ),
                        None,
                    )
                    if dataset:
                        results.append(
                            {
                                "year": item.year,
                                "theme": item.theme.key,
                                "status": "existing",
                                "dataset_id": dataset["id"],
                                "feature_count": dataset.get("feature_count"),
                            }
                        )
                        continue
                    dataset = upload_snapshot(
                        api_session,
                        base_url=base_url,
                        project_id=project_id,
                        area_id=area_id,
                        args=args,
                        snapshot=item,
                    )
                    existing.append(dataset)
                    results.append(
                        {
                            "year": item.year,
                            "theme": item.theme.key,
                            "status": "imported",
                            "dataset_id": dataset["id"],
                            "feature_count": dataset.get("feature_count"),
                        }
                    )
    except (OSError, RuntimeError, requests.RequestException, ValueError, KeyError, 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",
                "scope": args.scope_key,
                "municipality": args.municipality_name,
                "series": [series_key(theme, args.scope_key) for theme in definitions],
                "snapshots": results,
            },
            ensure_ascii=False,
            indent=2,
        )
    )
    return 0


if __name__ == "__main__":
    sys.exit(main())