geointel/scripts/provision_agricultural_parcel_history.py

"""Provision definitive Flemish agricultural-use parcels for an explicit scope.

The operator downloads only the official annual ALZ archives, retains those
archives as checksummed evidence, validates the GeoPackage schema and CRS,
clips parcel geometry to the persisted GeoIntel Area in EPSG:31370 and uploads
one canonical vector Dataset per year. It never writes to vector_features
directly and never uses the provisional current-campaign snapshot.
"""

from __future__ import annotations

import argparse
from dataclasses import dataclass
import hashlib
import json
import math
import os
import re
import sys
import tempfile
import zipfile
from datetime import date, datetime, timezone
from pathlib import Path
from typing import Any, Iterable
from urllib.parse import urlsplit

import requests
from pyproj import Transformer
from requests.adapters import HTTPAdapter
from shapely.geometry import GeometryCollection, MultiPolygon, Polygon, mapping, shape
from shapely.ops import transform as transform_geometry
from shapely.ops import unary_union
from shapely.validation import make_valid
from urllib3.util.retry import Retry


SCRIPT_DIR = Path(__file__).resolve().parent
if str(SCRIPT_DIR) not in sys.path:
    sys.path.insert(0, str(SCRIPT_DIR))

from geographic_scopes import GEOGRAPHIC_SCOPES, GeographicScope  # noqa: E402


DEFAULT_API_URL = "http://127.0.0.1:8000"
DEFAULT_SCOPE_KEY = "kempen-transport-region"
DEFAULT_OUTPUT_ROOT = Path("/app/storage/operator-evidence/agricultural-use-parcels")
CATALOG_URL = "https://landbouwcijfers.vlaanderen.be/open-geodata-landbouwgebruikspercelen"
DATA_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/open-geodata-landbouwgebruikspercelen"
ATTRIBUTION = "Agentschap Landbouw en Zeevisserij - Landbouwcijfers"
SOURCE_NAME = "agentschap_landbouw_zeevisserij_agricultural_parcels"
REFERENCE_LAYER_NAME = "agriculture"
SOURCE_CRS = "EPSG:31370"
OUTPUT_CRS = "EPSG:4326"
SCHEMA_VERSION = 1
MAX_ARCHIVE_BYTES = 250 * 1024 * 1024
MAX_EXTRACTED_BYTES = 2 * 1024 * 1024 * 1024
MAX_ARCHIVE_MEMBERS = 32
ARCHIVE_HOST = "www.landbouwvlaanderen.be"
ARCHIVE_PATH_PATTERN = re.compile(r"^/bestanden/gis/agpa_(20[0-9]{2})_(20[0-9]{2}-[0-9]{2}-[0-9]{2})_public\.zip$")

ARCHIVE_URLS = {
    2008: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2008_2022-03-23_public.zip",
    2009: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2009_2022-03-23_public.zip",
    2010: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2010_2022-03-23_public.zip",
    2011: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2011_2022-03-23_public.zip",
    2012: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2012_2022-03-23_public.zip",
    2013: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2013_2022-03-23_public.zip",
    2014: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2014_2022-03-23_public.zip",
    2015: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2015_2022-03-23_public.zip",
    2016: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2016_2022-03-23_public.zip",
    2017: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2017_2022-03-23_public.zip",
    2018: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2018_2022-03-23_public.zip",
    2019: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2019_2020-03-20_public.zip",
    2020: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2020_2021-03-19_public.zip",
    2021: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2021_2022-03-15_public.zip",
    2022: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2022_2023-06-26_public.zip",
    2023: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2023_2024-03-28_public.zip",
    2024: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2024_2025-03-27_public.zip",
    2025: "https://www.landbouwvlaanderen.be/bestanden/gis/agpa_2025_2026-05-13_public.zip",
}
SUPPORTED_YEARS = tuple(ARCHIVE_URLS)


@dataclass(frozen=True)
class AgriculturalReleaseConfig:
    year: int
    archive_url: str

    @property
    def definitive_version(self) -> str:
        return f"{self.year}-v3"

STABLE_REQUIRED_FIELDS = {
    "agpakey",
    "parcelnumber",
    "area_ha",
    "reference_id",
    "springcrop_code",
    "springcrop_title",
    "maincrop_code",
    "maincrop_title",
    "maincropgroup_title",
    "productionmethod_code",
    "productionmethod_title",
    "municipality_code",
    "municipality_title",
}

GROUP_ALIASES = {
    "grasland": "grassland",
    "mais": "maize",
    "granen/zaden/peulvruchten": "grains_seeds_legumes",
    "granen, zaden en peulvruchten": "grains_seeds_legumes",
    "aardappelen": "potatoes",
    "groenten/kruiden/sierplanten": "horticulture",
    "groenten, kruiden en sierplanten": "horticulture",
    "suikerbieten": "sugar_beets",
    "voedergewassen": "fodder_crops",
    "fruit en noten": "fruit_nuts",
    "overige gewassen": "other_crops",
    "vlas en hennep": "flax_hemp",
    "houtachtige gewassen": "woody_crops",
    "landbouwinfrastructuur": "agricultural_infrastructure",
    "water": "water",
}

METRIC_GROUPS = (
    ("grassland_area", "Grasland", ("grassland",)),
    ("maize_area", "Mais", ("maize",)),
    ("grains_seeds_legumes_area", "Granen, zaden en peulvruchten", ("grains_seeds_legumes",)),
    ("potatoes_area", "Aardappelen", ("potatoes",)),
    ("horticulture_area", "Groenten, kruiden en sierplanten", ("horticulture",)),
    ("sugar_beets_area", "Suikerbieten", ("sugar_beets",)),
    ("fodder_crops_area", "Voedergewassen", ("fodder_crops",)),
    ("fruit_nuts_area", "Fruit en noten", ("fruit_nuts",)),
    ("other_crops_area", "Overige gewassen", ("other_crops", "flax_hemp", "woody_crops")),
    ("agricultural_infrastructure_area", "Landbouwinfrastructuur", ("agricultural_infrastructure",)),
    ("agricultural_water_area", "Water binnen de aangifte", ("water",)),
)

TO_LAMBERT72 = Transformer.from_crs(OUTPUT_CRS, SOURCE_CRS, always_xy=True)
TO_WGS84 = Transformer.from_crs(SOURCE_CRS, OUTPUT_CRS, always_xy=True)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description="Provision definitive ALZ agricultural-use parcel history.")
    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL))
    parser.add_argument("--scope", choices=sorted(GEOGRAPHIC_SCOPES), default=DEFAULT_SCOPE_KEY)
    parser.add_argument("--years", default=",".join(str(year) for year in SUPPORTED_YEARS))
    parser.add_argument(
        "--archive-url",
        help="Exact official archive URL for one explicitly confirmed future definitive edition.",
    )
    parser.add_argument("--output-root", type=Path, default=Path(os.environ.get("GEOINTEL_AGRICULTURE_OUTPUT_ROOT", DEFAULT_OUTPUT_ROOT)))
    parser.add_argument("--request-timeout", type=int, default=900)
    parser.add_argument("--import-timeout", type=int, default=3600)
    parser.add_argument("--max-features", type=int, default=250_000)
    parser.add_argument("--max-archive-mb", type=int, default=250)
    parser.add_argument("--force", action="store_true", help="Redownload and rebuild retained evidence artifacts.")
    parser.add_argument("--fetch-only", action="store_true", help="Prepare evidence without importing Datasets.")
    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: Any, *, pretty: bool = False) -> None:
    path.parent.mkdir(parents=True, exist_ok=True)
    temporary = path.with_suffix(path.suffix + ".partial")
    temporary.write_text(
        json.dumps(payload, ensure_ascii=False, indent=2 if pretty else None, sort_keys=pretty, separators=None if pretty else (",", ":")),
        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-ALZ-Agricultural-Parcels-Operator/1.0"})
    adapter = HTTPAdapter(max_retries=retry)
    session.mount("https://", adapter)
    session.mount("http://", adapter)
    return session


def response_data(response: requests.Response) -> Any:
    try:
        payload = response.json()
    except requests.JSONDecodeError as exc:
        raise RuntimeError(f"GeoIntel API returned non-JSON content ({response.status_code})") from exc
    if not response.ok:
        if isinstance(payload, dict):
            message = payload.get("message") or payload.get("error") or response.text
        else:
            message = response.text
        raise RuntimeError(f"GeoIntel API request failed ({response.status_code}): {message}")
    if isinstance(payload, dict) and "data" in payload:
        return payload["data"]
    return payload


def api_items(session: requests.Session, url: str, timeout: int) -> list[dict[str, Any]]:
    items: list[dict[str, Any]] = []
    offset = 0
    expected_total: int | None = None
    while expected_total is None or offset < expected_total:
        response = session.get(url, params={"limit": 200, "offset": offset}, timeout=timeout)
        data = response_data(response)
        if isinstance(data, dict):
            page = data.get("items") or data.get("results") or []
            page_total = int(data.get("total") if data.get("total") is not None else len(page))
        else:
            page = data
            page_total = len(page) if isinstance(page, list) else 0
        if not isinstance(page, list):
            raise RuntimeError(f"Expected a list response from {url}")
        if expected_total is None:
            expected_total = page_total
        elif page_total != expected_total:
            raise RuntimeError("GeoIntel pagination total changed while reading the agricultural workspace")
        items.extend(item for item in page if isinstance(item, dict))
        if not page:
            break
        offset += len(page)
    if expected_total is not None and len(items) != expected_total:
        raise RuntimeError(f"GeoIntel list response for {url} returned {len(items)} of {expected_total} items")
    return items


def locate_workspace(
    session: requests.Session,
    base_url: str,
    scope: GeographicScope,
    timeout: int,
) -> tuple[str, str, Any, list[dict[str, Any]]]:
    projects = api_items(session, f"{base_url}/api/v1/projects", timeout)
    project = next((item for item in projects if item.get("name") == scope.project_name), None)
    if project is None:
        raise RuntimeError(f"Project {scope.project_name!r} is missing; provision the geographic scope first")
    project_id = str(project["id"])
    areas = api_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", timeout)
    area = next((item for item in areas if item.get("name") == scope.area_name), None)
    if area is None:
        raise RuntimeError(f"Area {scope.area_name!r} is missing from project {scope.project_name!r}")
    boundary = polygonal_geometry(shape(area.get("geometry")))
    if boundary is None:
        raise RuntimeError(f"Area {scope.area_name!r} has no valid polygon geometry")
    datasets = api_items(session, f"{base_url}/api/v1/projects/{project_id}/datasets", timeout)
    return project_id, str(area["id"]), boundary, datasets


def parse_years(raw: str) -> list[int]:
    try:
        years = sorted({int(value.strip()) for value in raw.split(",") if value.strip()})
    except ValueError as exc:
        raise ValueError("Years must be a comma-separated list of integers") from exc
    unsupported = [year for year in years if year not in ARCHIVE_URLS]
    if not years or unsupported:
        raise ValueError(f"Supported definitive years are {SUPPORTED_YEARS}; unsupported: {unsupported}")
    return years


def validate_archive_url(url: str, *, expected_year: int) -> str:
    parsed = urlsplit(url)
    match = ARCHIVE_PATH_PATTERN.fullmatch(parsed.path)
    if (
        parsed.scheme != "https"
        or parsed.hostname != ARCHIVE_HOST
        or parsed.port not in {None, 443}
        or parsed.username
        or parsed.password
        or parsed.query
        or parsed.fragment
        or not match
        or int(match.group(1)) != expected_year
    ):
        raise ValueError("Agricultural release archive is outside the official ALZ URL contract")
    try:
        datetime.strptime(match.group(2), "%Y-%m-%d")
    except ValueError as exc:
        raise ValueError("Agricultural release archive contains an invalid publication date") from exc
    return url


def resolve_release_config(year: int, *, archive_url: str | None = None) -> AgriculturalReleaseConfig:
    known_url = ARCHIVE_URLS.get(year)
    if known_url is not None:
        if archive_url is not None and archive_url != known_url:
            raise ValueError(f"The official retained archive identity for {year} may not be overridden")
        return AgriculturalReleaseConfig(year=year, archive_url=known_url)
    if year <= max(SUPPORTED_YEARS) or not archive_url:
        raise ValueError(
            f"One future definitive edition after {max(SUPPORTED_YEARS)} may be supplied with --archive-url"
        )
    return AgriculturalReleaseConfig(year=year, archive_url=validate_archive_url(archive_url, expected_year=year))


def resolve_release_configs(raw_years: str, *, archive_url: str | None = None) -> list[AgriculturalReleaseConfig]:
    if archive_url is None:
        return [resolve_release_config(year) for year in parse_years(raw_years)]
    try:
        years = sorted({int(value.strip()) for value in raw_years.split(",") if value.strip()})
    except ValueError as exc:
        raise ValueError("Years must be a comma-separated list of integers") from exc
    if len(years) != 1:
        raise ValueError("--archive-url requires exactly one explicitly selected definitive year")
    return [resolve_release_config(years[0], archive_url=archive_url)]


def polygonal_geometry(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 normalized_group_title(value: Any) -> str | None:
    if value is None:
        return None
    normalized = str(value).strip().casefold().replace("ï", "i")
    return GROUP_ALIASES.get(normalized, normalized.replace(" ", "_")) if normalized else None


def json_value(value: Any) -> Any:
    if value is None:
        return None
    if isinstance(value, (datetime, date)):
        return value.isoformat()
    if isinstance(value, float) and not math.isfinite(value):
        return None
    if isinstance(value, (str, int, float, bool)):
        return value
    if hasattr(value, "item"):
        try:
            return json_value(value.item())
        except (TypeError, ValueError):
            pass
    return str(value)


def download_archive(
    session: requests.Session,
    url: str,
    destination: Path,
    *,
    timeout: int,
    max_bytes: int,
    force: bool,
) -> dict[str, Any]:
    destination.parent.mkdir(parents=True, exist_ok=True)
    if destination.is_file() and not force:
        validate_archive(destination)
        return {"status": "reused", "sha256": sha256_file(destination), "size_bytes": destination.stat().st_size}
    temporary = destination.with_suffix(destination.suffix + ".partial")
    temporary.unlink(missing_ok=True)
    response = session.get(url, timeout=timeout, stream=True)
    response.raise_for_status()
    final_url = str(getattr(response, "url", "") or url)
    expected_year = int(ARCHIVE_PATH_PATTERN.fullmatch(urlsplit(url).path).group(1))
    if validate_archive_url(final_url, expected_year=expected_year) != url:
        raise RuntimeError("Official archive download redirected to a different release identity")
    content_length = int(response.headers.get("content-length") or 0)
    if content_length > max_bytes:
        raise RuntimeError(f"Official archive exceeds the configured {max_bytes // (1024 * 1024)} MiB limit")
    size = 0
    try:
        with temporary.open("wb") as handle:
            for chunk in response.iter_content(chunk_size=1024 * 1024):
                if not chunk:
                    continue
                size += len(chunk)
                if size > max_bytes:
                    raise RuntimeError(f"Official archive exceeded the configured {max_bytes // (1024 * 1024)} MiB limit while streaming")
                handle.write(chunk)
        validate_archive(temporary)
        temporary.replace(destination)
    except Exception:
        temporary.unlink(missing_ok=True)
        raise
    return {"status": "downloaded", "sha256": sha256_file(destination), "size_bytes": destination.stat().st_size}


def archive_geopackage_member(path: Path) -> str:
    with zipfile.ZipFile(path) as archive:
        entries = archive.infolist()
        if len(entries) > MAX_ARCHIVE_MEMBERS:
            raise RuntimeError(f"Official archive contains more than {MAX_ARCHIVE_MEMBERS} members")
        if sum(item.file_size for item in entries if not item.is_dir()) > MAX_EXTRACTED_BYTES:
            raise RuntimeError("Official archive exceeds the extracted-size safety limit")
        members = [item.filename for item in entries if not item.is_dir() and item.filename.lower().endswith(".gpkg")]
    if len(members) != 1:
        raise RuntimeError(f"Official archive must contain exactly one GeoPackage; found {len(members)}")
    member = members[0]
    if Path(member).name != member or ".." in Path(member).parts:
        raise RuntimeError("Official archive contains an unsafe GeoPackage path")
    return member


def validate_archive(path: Path) -> str:
    try:
        return archive_geopackage_member(path)
    except zipfile.BadZipFile as exc:
        raise RuntimeError(f"Official source archive {path.name} is not a valid ZIP") from exc


def extract_geopackage(archive_path: Path, destination_dir: Path) -> Path:
    member = validate_archive(archive_path)
    destination = destination_dir / Path(member).name
    with zipfile.ZipFile(archive_path) as archive, archive.open(member) as source, destination.open("wb") as target:
        extracted_bytes = 0
        for chunk in iter(lambda: source.read(1024 * 1024), b""):
            extracted_bytes += len(chunk)
            if extracted_bytes > MAX_EXTRACTED_BYTES:
                raise RuntimeError("Official GeoPackage exceeded the extracted-size safety limit while streaming")
            target.write(chunk)
    if destination.stat().st_size == 0:
        raise RuntimeError("Extracted official GeoPackage is empty")
    return destination


def load_pyogrio():
    try:
        import pyogrio  # type: ignore[import-not-found]
    except ImportError as exc:
        raise RuntimeError("Agricultural parcel preparation requires the optional GeoIntel GIS dependencies (pyogrio/geopandas)") from exc
    return pyogrio


def inspect_geopackage(pyogrio, path: Path) -> tuple[str, set[str], int]:
    layers = pyogrio.list_layers(path)
    polygon_layers = [str(row[0]) for row in layers if "Polygon" in str(row[1])]
    if len(polygon_layers) != 1:
        raise RuntimeError(f"Expected one polygon layer in {path.name}; found {polygon_layers}")
    layer = polygon_layers[0]
    info = pyogrio.read_info(path, layer=layer)
    crs = str(info.get("crs") or "")
    if "31370" not in crs:
        raise RuntimeError(f"Expected EPSG:31370 agricultural parcels; received {crs or 'no CRS'}")
    info_fields = info.get("fields")
    fields = {str(value).lower() for value in info_fields} if info_fields is not None else set()
    missing = sorted(STABLE_REQUIRED_FIELDS - fields)
    if missing:
        raise RuntimeError(f"Agricultural parcel schema is missing stable fields: {', '.join(missing)}")
    return layer, fields, int(info.get("features") or 0)


def build_crop_code_list(records: Iterable[dict[str, Any]], *, year: int) -> dict[str, Any]:
    crops: dict[tuple[str, str, str], dict[str, Any]] = {}
    groups: set[str] = set()
    title_by_code: dict[str, set[str]] = {}
    for record in records:
        code = str(json_value(record.get("maincrop_code")) or "").strip()
        title = str(json_value(record.get("maincrop_title")) or "").strip()
        group_title = str(json_value(record.get("maincropgroup_title")) or "").strip()
        if not code and not title and not group_title:
            continue
        crops[(code, title, group_title)] = {"code": code or None, "title": title or None, "group_title": group_title or None}
        if group_title:
            groups.add(group_title)
        if code and title:
            title_by_code.setdefault(code, set()).add(title)
    conflicts = {code: sorted(titles) for code, titles in title_by_code.items() if len(titles) > 1}
    return {
        "year": year,
        "source": ATTRIBUTION,
        "generated_from": "full official annual GeoPackage attributes",
        "crop_entries": sorted(crops.values(), key=lambda item: (str(item["code"]), str(item["title"]), str(item["group_title"]))),
        "main_crop_groups": sorted(groups),
        "code_title_conflicts": conflicts,
        "historical_comparison_rule": "Use maincropgroup_title for comparable grouped hectares; detailed crop code/title remains source-faithful per year.",
    }


def dataframe_records(frame, columns: Iterable[str]) -> Iterable[dict[str, Any]]:
    for row in frame.itertuples(index=False, name=None):
        yield {column: json_value(value) for column, value in zip(columns, row)}


def normalize_frame(frame, *, year: int, boundary_lambert72, max_features: int) -> tuple[list[dict[str, Any]], dict[str, Any]]:
    features: list[dict[str, Any]] = []
    seen_ids: set[str] = set()
    clipped_count = 0
    source_area_ha = 0.0
    exact_area_ha = 0.0
    for _, row in frame.iterrows():
        source_geometry = polygonal_geometry(row.geometry)
        if source_geometry is None or not source_geometry.intersects(boundary_lambert72):
            continue
        clipped = polygonal_geometry(source_geometry.intersection(boundary_lambert72))
        if clipped is None or clipped.area <= 0:
            continue
        properties = {str(column).lower(): json_value(row[column]) for column in frame.columns if str(column).lower() != "geometry"}
        agpa_key = str(properties.get("agpakey") or "").strip()
        if not agpa_key:
            raise RuntimeError(f"Agricultural parcel {year} contains a feature without agpakey")
        source_feature_id = f"alz:{year}:{agpa_key}"
        if source_feature_id in seen_ids:
            raise RuntimeError(f"Agricultural parcel {year} contains duplicate agpakey {agpa_key}")
        seen_ids.add(source_feature_id)
        clipped_wgs84 = polygonal_geometry(transform_geometry(TO_WGS84.transform, clipped))
        if clipped_wgs84 is None:
            raise RuntimeError(f"Agricultural parcel {agpa_key} could not be transformed to EPSG:4326")
        source_area = float(source_geometry.area) / 10_000.0
        exact_area = float(clipped.area) / 10_000.0
        source_area_ha += source_area
        exact_area_ha += exact_area
        was_clipped = not source_geometry.within(boundary_lambert72)
        clipped_count += int(was_clipped)
        properties.update(
            {
                "source_feature_id": source_feature_id,
                "source_year": year,
                "source_agpa_key": agpa_key,
                "main_crop_group_key": normalized_group_title(properties.get("maincropgroup_title")),
                "source_geometry_area_ha": round(source_area, 8),
                "clipped_area_ha": round(exact_area, 8),
                "geometry_was_clipped": was_clipped,
                "historical_parcel_identity_stable": False,
            }
        )
        features.append({"type": "Feature", "id": source_feature_id, "geometry": mapping(clipped_wgs84), "properties": properties})
        if len(features) > max_features:
            raise RuntimeError(f"Clipped agricultural parcel count exceeds the configured limit of {max_features}")
    if not features:
        raise RuntimeError(f"Official agricultural parcel archive {year} has no features inside the persisted scope")
    return features, {
        "feature_count": len(features),
        "clipped_feature_count": clipped_count,
        "source_geometry_area_ha": round(source_area_ha, 4),
        "clipped_area_ha": round(exact_area_ha, 4),
    }


def artifact_paths(
    output_root: Path,
    scope_key: str,
    year: int,
    *,
    archive_url: str | None = None,
) -> dict[str, Path]:
    release = resolve_release_config(year, archive_url=archive_url)
    directory = output_root / scope_key / str(year)
    return {
        "directory": directory,
        "archive": directory / Path(urlsplit(release.archive_url).path).name,
        "artifact": directory / f"agricultural_use_parcels_{year}_{scope_key}.geojson",
        "codelist": directory / f"agricultural_use_parcels_{year}_crop_codes.json",
        "manifest": directory / f"agricultural_use_parcels_{year}_{scope_key}.manifest.json",
    }


def reusable_artifact(
    paths: dict[str, Path],
    *,
    year: int,
    scope_key: str,
    archive_url: str | None = None,
) -> dict[str, Any] | None:
    if not all(paths[key].is_file() for key in ("archive", "artifact", "codelist", "manifest")):
        return None
    manifest = json.loads(paths["manifest"].read_text(encoding="utf-8"))
    if manifest.get("schema_version") != SCHEMA_VERSION or manifest.get("year") != year or manifest.get("scope_key") != scope_key:
        return None
    release = resolve_release_config(year, archive_url=archive_url)
    if manifest.get("source_url") != release.archive_url:
        return None
    if manifest.get("source_archive_sha256") != sha256_file(paths["archive"]):
        return None
    if manifest.get("artifact_sha256") != sha256_file(paths["artifact"]):
        return None
    if manifest.get("crop_code_list_sha256") != sha256_file(paths["codelist"]):
        return None
    validate_archive(paths["archive"])
    return manifest


def prepare_year(
    session: requests.Session,
    *,
    year: int,
    scope: GeographicScope,
    boundary_wgs84,
    output_root: Path,
    request_timeout: int,
    max_archive_bytes: int,
    max_features: int,
    force: bool,
    archive_url: str | None = None,
) -> tuple[dict[str, Path], dict[str, Any]]:
    release = resolve_release_config(year, archive_url=archive_url)
    paths = artifact_paths(output_root, scope.key, year, archive_url=release.archive_url)
    paths["directory"].mkdir(parents=True, exist_ok=True)
    if not force:
        reused = reusable_artifact(paths, year=year, scope_key=scope.key, archive_url=release.archive_url)
        if reused is not None:
            return paths, reused
    download = download_archive(
        session,
        release.archive_url,
        paths["archive"],
        timeout=request_timeout,
        max_bytes=max_archive_bytes,
        force=force,
    )
    boundary_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, boundary_wgs84))
    if boundary_lambert72 is None:
        raise RuntimeError(f"Scope {scope.key} could not be transformed to EPSG:31370")
    pyogrio = load_pyogrio()
    with tempfile.TemporaryDirectory(prefix=f"agpa-{year}-", dir=paths["directory"]) as temporary:
        gpkg_path = extract_geopackage(paths["archive"], Path(temporary))
        layer, fields, source_feature_count = inspect_geopackage(pyogrio, gpkg_path)
        code_columns = [field for field in ("maincrop_code", "maincrop_title", "maincropgroup_title") if field in fields]
        code_frame = pyogrio.read_dataframe(gpkg_path, layer=layer, columns=code_columns, read_geometry=False)
        code_list = build_crop_code_list(dataframe_records(code_frame, code_columns), year=year)
        write_json_atomic(paths["codelist"], code_list, pretty=True)
        frame = pyogrio.read_dataframe(gpkg_path, layer=layer, bbox=boundary_lambert72.bounds)
        features, summary = normalize_frame(frame, year=year, boundary_lambert72=boundary_lambert72, max_features=max_features)
    artifact = {"type": "FeatureCollection", "name": paths["artifact"].stem, "features": features}
    write_json_atomic(paths["artifact"], artifact)
    manifest = {
        "schema_version": SCHEMA_VERSION,
        "generated_at": utc_now(),
        "year": year,
        "scope_key": scope.key,
        "scope_name": scope.display_name,
        "scope_type": scope.scope_type,
        "member_nis_codes": list(scope.nis_codes),
        "source_url": release.archive_url,
        "catalog_url": CATALOG_URL,
        "data_catalog_url": DATA_CATALOG_URL,
        "attribution": ATTRIBUTION,
        "source_crs": SOURCE_CRS,
        "output_crs": OUTPUT_CRS,
        "source_archive_path": str(paths["archive"]),
        "source_archive_sha256": download["sha256"],
        "source_archive_size_bytes": download["size_bytes"],
        "source_feature_count": source_feature_count,
        "source_fields": sorted(fields),
        "crop_code_list_path": str(paths["codelist"]),
        "crop_code_list_sha256": sha256_file(paths["codelist"]),
        "artifact_path": str(paths["artifact"]),
        "artifact_sha256": sha256_file(paths["artifact"]),
        **summary,
        "limitations": [
            "Parcel identities are not stable across campaign years; evolution compares grouped area totals, not parcel lineage.",
            "The layer describes declared agricultural use at the annual reference deadline and can include water, hedges, buildings and infrastructure.",
            "Detailed crop codes and titles remain source-faithful per year; only official main-crop groups are used for comparable historical metrics.",
            "The current-campaign provisional snapshot is intentionally excluded.",
        ],
    }
    write_json_atomic(paths["manifest"], manifest, pretty=True)
    return paths, manifest


def selection_metrics() -> list[dict[str, Any]]:
    warning = "Historische evolutie vergelijkt officiele hoofdteeltgroepen; individuele perceelidentiteiten zijn niet stabiel tussen campagnejaren."
    return [
        {
            "metric_key": metric_key,
            "method": "intersection_area",
            "label": label,
            "unit": "ha",
            "geometry_dimension": 2,
            "filter_property": "main_crop_group_key",
            "filter_values": list(group_keys),
            "warning": warning,
        }
        for metric_key, label, group_keys in METRIC_GROUPS
    ]


def series_key(scope: GeographicScope) -> str:
    return f"alz:agricultural-use-parcels:{scope.key}"


def upload_artifact(
    session: requests.Session,
    *,
    base_url: str,
    project_id: str,
    area_id: str,
    scope: GeographicScope,
    year: int,
    paths: dict[str, Path],
    manifest: dict[str, Any],
    timeout: int,
) -> dict[str, Any]:
    source_metadata = {
        "provider": "Agentschap Landbouw en Zeevisserij",
        "theme": "agriculture",
        "layer_name": f"Landbouwgebruikspercelen {year}",
        "authority_level": "authoritative",
        "coverage_scope": scope.scope_type,
        "scope_key": scope.key,
        "scope_name": scope.display_name,
        "member_nis_codes": list(scope.nis_codes),
        "feature_count": manifest["feature_count"],
        "geometry_clipped_to_area": True,
        "identity_stable": False,
        "semantic_metrics": False,
        "attribution": ATTRIBUTION,
        "catalog_url": CATALOG_URL,
        "selection_aggregation": {
            "metric_key": "declared_agricultural_use_area",
            "method": "intersection_area",
            "label": "Aangegeven gebruiksoppervlakte",
            "unit": "ha",
            "geometry_dimension": 2,
            "warning": "De aangifte bevat naast teelten ook onder meer water, hagen, gebouwen en landbouwinfrastructuur; dit is geen eigendoms- of kadastrale oppervlakte.",
        },
        "selection_metrics": selection_metrics(),
    }
    provenance_metadata = {
        "operator_tool": "provision_agricultural_parcel_history.py",
        "operator_explicit_fetch": True,
        "geometry_clipped_to_area": True,
        "source_archive_url": manifest["source_url"],
        "source_archive_path": str(paths["archive"]),
        "source_archive_sha256": manifest["source_archive_sha256"],
        "crop_code_list_path": str(paths["codelist"]),
        "crop_code_list_sha256": manifest["crop_code_list_sha256"],
        "manifest_path": str(paths["manifest"]),
        "artifact_sha256": manifest["artifact_sha256"],
        "catalog_url": CATALOG_URL,
        "limitations": manifest["limitations"],
    }
    with paths["artifact"].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": SOURCE_NAME,
                "reference_layer_name": REFERENCE_LAYER_NAME,
                "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(scope),
                "observed_at": f"{year}-01-01T00:00:00Z",
                "temporal_granularity": "year",
                "source_version": f"{year}-definitive",
            },
            files={"file": (paths["artifact"].name, handle, "application/geo+json")},
            timeout=timeout,
        )
    return response_data(response)


def existing_dataset_for_year(datasets: list[dict[str, Any]], *, area_id: str, year: int) -> dict[str, Any] | None:
    return next(
        (
            dataset
            for dataset in datasets
            if dataset.get("source_name") == SOURCE_NAME
            and dataset.get("source_version") == f"{year}-definitive"
            and str(dataset.get("area_id") or "") == area_id
        ),
        None,
    )


def main() -> int:
    args = parse_args()
    try:
        releases = resolve_release_configs(args.years, archive_url=args.archive_url)
        if args.max_features < 1 or args.max_archive_mb < 1:
            raise ValueError("Feature and archive safety limits must be positive")
        scope = GEOGRAPHIC_SCOPES[args.scope]
        base_url = args.base_url.rstrip("/")
        with requests.Session() as api_session:
            project_id, area_id, boundary, datasets = locate_workspace(api_session, base_url, scope, args.import_timeout)
            results: list[dict[str, Any]] = []
            with build_session() as official_session:
                for release in releases:
                    year = release.year
                    paths, manifest = prepare_year(
                        official_session,
                        year=year,
                        scope=scope,
                        boundary_wgs84=boundary,
                        output_root=args.output_root,
                        request_timeout=args.request_timeout,
                        max_archive_bytes=min(args.max_archive_mb * 1024 * 1024, MAX_ARCHIVE_BYTES),
                        max_features=args.max_features,
                        force=args.force,
                        archive_url=release.archive_url,
                    )
                    existing = existing_dataset_for_year(datasets, area_id=area_id, year=year)
                    if existing is not None:
                        persisted_checksum = str(existing.get("checksum_sha256") or "")
                        if persisted_checksum and persisted_checksum != manifest["artifact_sha256"]:
                            raise RuntimeError(f"A different {year} agricultural parcel artifact is already persisted; refusing silent replacement")
                        result = {"year": year, "status": "existing", "dataset_id": existing["id"], "feature_count": existing.get("feature_count")}
                    elif args.fetch_only:
                        result = {"year": year, "status": "prepared", "artifact_path": str(paths["artifact"]), "feature_count": manifest["feature_count"]}
                    else:
                        dataset = upload_artifact(
                            api_session,
                            base_url=base_url,
                            project_id=project_id,
                            area_id=area_id,
                            scope=scope,
                            year=year,
                            paths=paths,
                            manifest=manifest,
                            timeout=args.import_timeout,
                        )
                        result = {"year": year, "status": "imported", "dataset_id": dataset["id"], "feature_count": dataset.get("feature_count")}
                    results.append({**result, "area_ha": manifest["clipped_area_ha"], "manifest_path": str(paths["manifest"])})
    except (KeyError, OSError, RuntimeError, ValueError, zipfile.BadZipFile, requests.RequestException) as exc:
        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
        return 1
    print(
        json.dumps(
            {
                "status": "ok",
                "scope": scope.key,
                "project": scope.project_name,
                "area": scope.area_name,
                "series_key": series_key(scope),
                "years": results,
                "historical_identity_stable": False,
            },
            ensure_ascii=False,
            indent=2,
        )
    )
    return 0


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