geointel/scripts/provision_mol_historical_landuse.py

"""Provision official historical land-use theme snapshots for Mol.

The command reads the Digitaal Vlaanderen historical land-use WFS for 1778,
1873 and 1969, clips features to the official Mol boundary, separates the
supported map themes and imports every snapshot through the normal dataset API.
It is an explicit, idempotent operator command and never runs at startup.
"""

from __future__ import annotations

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

import requests
from requests.adapters import HTTPAdapter
from shapely.geometry import mapping, shape
from shapely.validation import make_valid
from urllib3.util.retry import Retry


MUNICIPALITY_NAME = "Mol"
MUNICIPALITY_NIS_CODE = "13025"
PROJECT_NAME = "Mol Municipality Workbench"
DEFAULT_API_URL = "http://127.0.0.1:8000"
DEFAULT_OUTPUT_DIR = Path("/app/storage/operator-data/mol-historical-landuse")
DEFAULT_BOUNDARY_PATH = Path("/app/storage/operator-data/mol-municipality/mol_municipality_boundary.geojson")
WFS_URL = "https://geo.api.vlaanderen.be/HistLandgebruik/wfs"
ATTRIBUTION = "Bron: Historisch landgebruik Vlaanderen, Digitaal Vlaanderen"
COLLECTIONS = {1778: "HistLandgebruik:Lgbrk1778", 1873: "HistLandgebruik:Lgbrk1873", 1969: "HistLandgebruik:Lgbrk1969"}


@dataclass(frozen=True)
class ThemeDefinition:
    key: str
    label: str
    matches: Callable[[str], bool]
    filter_value: str
    filter_mode: str


THEMES = (
    ThemeDefinition("buildings", "Historische bebouwing", lambda value: value.startswith("bebouwing"), "bebouwing*", "like"),
    ThemeDefinition("forest", "Historisch bos", lambda value: value.startswith("bos-") or value == "bos", "bos*", "like"),
    ThemeDefinition("water", "Historisch water", lambda value: value == "water", "water", "equal"),
    ThemeDefinition("roads", "Historische wegen", lambda value: value.startswith("weg"), "weg*", "like"),
)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description="Provision official historical land-use snapshots for Mol.")
    parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL))
    parser.add_argument("--project-name", default=PROJECT_NAME)
    parser.add_argument("--years", default="1778,1873,1969")
    parser.add_argument("--themes", default="buildings,forest,water,roads")
    parser.add_argument("--output-dir", type=Path, default=Path(os.environ.get("MOL_HISTORICAL_LANDUSE_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)))
    parser.add_argument("--boundary-path", type=Path, default=Path(os.environ.get("MOL_BOUNDARY_PATH", DEFAULT_BOUNDARY_PATH)))
    parser.add_argument("--page-size", type=int, default=200)
    parser.add_argument("--max-features", type=int, default=100000)
    parser.add_argument("--simplify-tolerance-degrees", type=float, default=0.00001)
    parser.add_argument("--request-timeout", type=int, default=180)
    parser.add_argument("--import-timeout", type=int, default=1800)
    parser.add_argument("--force", action="store_true")
    parser.add_argument("--fetch-only", action="store_true")
    return parser.parse_args()


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", "POST"}),
        raise_on_status=True,
    )
    session = requests.Session()
    session.headers.update({"User-Agent": "GeoIntel-Mol-Historical-Landuse-Operator/1.0"})
    adapter = HTTPAdapter(max_retries=retry)
    session.mount("https://", adapter)
    session.mount("http://", adapter)
    return session


def load_boundary(path: Path):
    if not path.exists():
        raise RuntimeError(f"Mol boundary is missing at {path}; run provision_mol_municipality_workspace.py first")
    payload = json.loads(path.read_text(encoding="utf-8"))
    features = payload.get("features") or []
    if len(features) != 1:
        raise RuntimeError("Mol boundary artifact must contain exactly one feature")
    boundary = shape(features[0]["geometry"])
    if not boundary.is_valid:
        boundary = make_valid(boundary)
    if boundary.is_empty or not boundary.is_valid:
        raise RuntimeError("Mol boundary artifact is invalid")
    return boundary


def wfs_filter_xml(definition: ThemeDefinition, bounds: tuple[float, float, float, float]) -> str:
    min_x, min_y, max_x, max_y = bounds
    comparison = (
        f"<fes:PropertyIsLike wildCard='*' singleChar='?' escapeChar='!'>"
        f"<fes:ValueReference>KLASSE</fes:ValueReference><fes:Literal>{definition.filter_value}</fes:Literal>"
        f"</fes:PropertyIsLike>"
        if definition.filter_mode == "like"
        else (
            "<fes:PropertyIsEqualTo><fes:ValueReference>KLASSE</fes:ValueReference>"
            f"<fes:Literal>{definition.filter_value}</fes:Literal></fes:PropertyIsEqualTo>"
        )
    )
    return (
        "<fes:Filter xmlns:fes='http://www.opengis.net/fes/2.0' xmlns:gml='http://www.opengis.net/gml/3.2'>"
        "<fes:And><fes:BBOX><fes:ValueReference>SHAPE</fes:ValueReference>"
        "<gml:Envelope srsName='EPSG:4326'>"
        f"<gml:lowerCorner>{min_x:.8f} {min_y:.8f}</gml:lowerCorner>"
        f"<gml:upperCorner>{max_x:.8f} {max_y:.8f}</gml:upperCorner>"
        "</gml:Envelope></fes:BBOX>"
        f"{comparison}</fes:And></fes:Filter>"
    )


def wfs_request_xml(
    *,
    collection: str,
    filter_xml: str,
    page_size: int,
    start_index: int,
) -> str:
    return (
        "<?xml version='1.0' encoding='UTF-8'?>"
        "<wfs:GetFeature service='WFS' version='2.0.0' outputFormat='application/json' "
        f"count='{page_size}' startIndex='{start_index}' "
        "xmlns:wfs='http://www.opengis.net/wfs/2.0' "
        "xmlns:fes='http://www.opengis.net/fes/2.0' "
        "xmlns:gml='http://www.opengis.net/gml/3.2' "
        "xmlns:HistLandgebruik='https://geo.api.vlaanderen.be/HistLandgebruik'>"
        f"<wfs:Query typeNames='{collection}' srsName='EPSG:4326'>{filter_xml}</wfs:Query>"
        "</wfs:GetFeature>"
    )


def fetch_wfs_page(
    session: requests.Session,
    *,
    collection: str,
    filter_xml: str,
    page_size: int,
    start_index: int,
    timeout: int,
) -> 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()
    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 features


def fetch_year(
    session: requests.Session,
    year: int,
    boundary,
    definitions: list[ThemeDefinition],
    *,
    page_size: int,
    max_features: int,
    simplify_tolerance_degrees: float,
    timeout: int,
):
    collection = COLLECTIONS[year]
    features_by_theme: dict[str, list[dict[str, Any]]] = {definition.key: [] for definition in definitions}
    for definition in definitions:
        filter_xml = wfs_filter_xml(definition, boundary.bounds)
        start_index = 0
        seen: set[str] = set()
        while True:
            page = fetch_wfs_page(
                session,
                collection=collection,
                filter_xml=filter_xml,
                page_size=page_size,
                start_index=start_index,
                timeout=timeout,
            )
            for raw_feature in page:
                feature_id = str(raw_feature.get("id") or "")
                if not feature_id or feature_id in seen:
                    continue
                seen.add(feature_id)
                properties = raw_feature.get("properties") or {}
                landuse_class = str(properties.get("KLASSE") or "").strip().lower()
                if not definition.matches(landuse_class):
                    continue
                geometry = shape(raw_feature["geometry"])
                if not geometry.is_valid:
                    geometry = make_valid(geometry)
                geometry = geometry.intersection(boundary)
                if geometry.is_empty:
                    continue
                if not geometry.is_valid:
                    geometry = make_valid(geometry)
                if geometry.is_empty or not geometry.is_valid:
                    continue
                if simplify_tolerance_degrees > 0:
                    geometry = geometry.simplify(simplify_tolerance_degrees, preserve_topology=True)
                theme_features = features_by_theme[definition.key]
                if len(theme_features) >= max_features:
                    raise RuntimeError(
                        f"Historical land use {definition.key} {year} exceeds the {max_features} feature safety limit"
                    )
                theme_features.append({**raw_feature, "id": feature_id, "geometry": mapping(geometry)})
            if len(page) < page_size:
                break
            start_index += len(page)
    for definition in definitions:
        if not features_by_theme[definition.key]:
            raise RuntimeError(f"Historical land-use WFS returned no {definition.key} features for Mol in {year}")
    return features_by_theme


def write_theme_snapshot(
    session: requests.Session,
    *,
    year: int,
    definition: ThemeDefinition,
    boundary,
    path: Path,
    page_size: int,
    max_features: int,
    simplify_tolerance_degrees: float,
    timeout: int,
) -> int:
    collection = COLLECTIONS[year]
    filter_xml = wfs_filter_xml(definition, boundary.bounds)
    start_index = 0
    seen: set[str] = set()
    feature_count = 0
    first_feature = True
    try:
        with path.open("w", encoding="utf-8") as output:
            output.write(
                json.dumps(
                    {
                        "type": "FeatureCollection",
                        "name": f"{definition.label} - Mol {year}",
                        "municipality": MUNICIPALITY_NAME,
                        "nis_code": MUNICIPALITY_NIS_CODE,
                        "observation_year": year,
                        "attribution": ATTRIBUTION,
                    },
                    ensure_ascii=False,
                    separators=(",", ":"),
                )[:-1]
            )
            output.write(',"features":[')
            while True:
                page = fetch_wfs_page(
                    session,
                    collection=collection,
                    filter_xml=filter_xml,
                    page_size=page_size,
                    start_index=start_index,
                    timeout=timeout,
                )
                for raw_feature in page:
                    feature_id = str(raw_feature.get("id") or "")
                    if not feature_id or feature_id in seen:
                        continue
                    seen.add(feature_id)
                    properties = dict(raw_feature.get("properties") or {})
                    landuse_class = str(properties.get("KLASSE") or "").strip().lower()
                    if not definition.matches(landuse_class):
                        continue
                    geometry = shape(raw_feature["geometry"])
                    if not geometry.is_valid:
                        geometry = make_valid(geometry)
                    geometry = geometry.intersection(boundary)
                    if geometry.is_empty:
                        continue
                    if not geometry.is_valid:
                        geometry = make_valid(geometry)
                    if geometry.is_empty or not geometry.is_valid:
                        continue
                    if simplify_tolerance_degrees > 0:
                        geometry = geometry.simplify(simplify_tolerance_degrees, preserve_topology=True)
                    properties.update(
                        {
                            "source_name": "historical_landuse",
                            "source_feature_id": feature_id,
                            "reference_layer_name": definition.key,
                            "authority_level": "authoritative",
                            "municipality": MUNICIPALITY_NAME,
                            "nis_code": MUNICIPALITY_NIS_CODE,
                            "observation_year": year,
                            "historical_landuse_class": landuse_class,
                            "attribution": ATTRIBUTION,
                        }
                    )
                    prepared = {
                        "type": "Feature",
                        "id": feature_id,
                        "geometry": mapping(geometry),
                        "properties": properties,
                    }
                    if not first_feature:
                        output.write(",")
                    output.write(json.dumps(prepared, ensure_ascii=False, separators=(",", ":")))
                    first_feature = False
                    feature_count += 1
                    if feature_count > max_features:
                        raise RuntimeError(
                            f"Historical land use {definition.key} {year} exceeds the {max_features} feature safety limit"
                        )
                if len(page) < page_size:
                    break
                start_index += len(page)
            output.write("]}")
    except Exception:
        path.unlink(missing_ok=True)
        path.with_suffix(".manifest.json").unlink(missing_ok=True)
        raise
    if feature_count == 0:
        path.unlink(missing_ok=True)
        raise RuntimeError(f"Historical land-use WFS returned no {definition.key} features for Mol in {year}")
    path.with_suffix(".manifest.json").write_text(
        json.dumps(
            {
                "year": year,
                "theme": definition.key,
                "feature_count": feature_count,
                "collection": collection,
                "simplify_tolerance_degrees": simplify_tolerance_degrees,
                "generated_at": datetime.now(timezone.utc).isoformat(),
            },
            ensure_ascii=False,
            indent=2,
        ),
        encoding="utf-8",
    )
    return feature_count


def build_theme_snapshot(year: int, definition: ThemeDefinition, source_features: list[dict[str, Any]]) -> dict[str, Any]:
    features: list[dict[str, Any]] = []
    for source_feature in source_features:
        properties = dict(source_feature.get("properties") or {})
        landuse_class = str(properties.get("KLASSE") or "").strip().lower()
        if not definition.matches(landuse_class):
            continue
        feature_id = str(source_feature["id"])
        properties.update(
            {
                "source_name": "historical_landuse",
                "source_feature_id": feature_id,
                "reference_layer_name": definition.key,
                "authority_level": "authoritative",
                "municipality": MUNICIPALITY_NAME,
                "nis_code": MUNICIPALITY_NIS_CODE,
                "observation_year": year,
                "historical_landuse_class": landuse_class,
                "attribution": ATTRIBUTION,
            }
        )
        features.append(
            {"type": "Feature", "id": feature_id, "geometry": source_feature["geometry"], "properties": properties}
        )
    if not features:
        raise RuntimeError(f"No {definition.key} features were classified for Mol in {year}")
    return {
        "type": "FeatureCollection",
        "name": f"{definition.label} - Mol {year}",
        "features": features,
        "municipality": MUNICIPALITY_NAME,
        "nis_code": MUNICIPALITY_NIS_CODE,
        "observation_year": year,
        "attribution": ATTRIBUTION,
    }


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 locate_workspace(session: requests.Session, base_url: str, project_name: str, 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") == project_name), None)
    if not project:
        raise RuntimeError(f"Project {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 "gemeente mol" in str(item.get("name", "")).lower()), None)
    if not area:
        raise RuntimeError("Official Mol area 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,
    year: int,
    definition: ThemeDefinition,
    path: Path,
    simplify_tolerance_degrees: float,
    timeout: int,
) -> dict[str, Any]:
    observed_at = f"{year}-01-01T00:00:00Z"
    series_key = f"digitaal-vlaanderen:historical-landuse:{definition.key}:mol"
    source_metadata = {
        "provider": "Digitaal Vlaanderen",
        "collection": COLLECTIONS[year],
        "authority_level": "authoritative",
        "coverage_scope": "municipality",
        "municipality": MUNICIPALITY_NAME,
        "nis_code": MUNICIPALITY_NIS_CODE,
        "attribution": ATTRIBUTION,
        "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 en karteermethodes verschillen per bronjaar; interpreteer trends binnen die methodologische context.",
        },
    }
    provenance_metadata = {
        "operator_tool": "provision_mol_historical_landuse.py",
        "operator_explicit_fetch": True,
        "wfs_url": WFS_URL,
        "collection": COLLECTIONS[year],
        "geometry_simplification_tolerance_degrees": simplify_tolerance_degrees,
        "generated_at": datetime.now(timezone.utc).isoformat(),
    }
    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 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 COLLECTIONS]
    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

    args.output_dir.mkdir(parents=True, exist_ok=True)
    results: list[dict[str, Any]] = []
    try:
        boundary = load_boundary(args.boundary_path)
        prepared: list[tuple[int, ThemeDefinition, Path, int]] = []
        with build_session() as source_session:
            for year in years:
                target_paths = {
                    definition.key: args.output_dir / f"mol_historical_{definition.key}_{year}.geojson"
                    for definition in definitions
                }
                for definition in definitions:
                    path = target_paths[definition.key]
                    manifest_path = path.with_suffix(".manifest.json")
                    if args.force or not path.exists() or not manifest_path.exists():
                        count = write_theme_snapshot(
                            source_session,
                            year=year,
                            definition=definition,
                            boundary=boundary,
                            path=path,
                            page_size=args.page_size,
                            max_features=args.max_features,
                            simplify_tolerance_degrees=args.simplify_tolerance_degrees,
                            timeout=args.request_timeout,
                        )
                    else:
                        count = int(json.loads(manifest_path.read_text(encoding="utf-8"))["feature_count"])
                    prepared.append((year, definition, path, count))

        if args.fetch_only:
            results = [
                {"year": year, "theme": definition.key, "path": str(path), "feature_count": count, "status": "prepared"}
                for year, definition, path, count 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.project_name, args.import_timeout)
                for year, definition, path, count in prepared:
                    series_key = f"digitaal-vlaanderen:historical-landuse:{definition.key}:mol"
                    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:
                        results.append({"year": year, "theme": definition.key, "dataset_id": dataset["id"], "feature_count": dataset.get("feature_count"), "status": "existing"})
                        continue
                    dataset = upload_snapshot(
                        api_session,
                        base_url,
                        project_id,
                        area_id,
                        year,
                        definition,
                        path,
                        args.simplify_tolerance_degrees,
                        args.import_timeout,
                    )
                    results.append({"year": year, "theme": definition.key, "dataset_id": dataset["id"], "feature_count": dataset.get("feature_count"), "status": "imported"})
    except (OSError, RuntimeError, requests.RequestException, ValueError, KeyError) as exc:
        print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr)
        return 1

    print(json.dumps({"status": "ok", "municipality": MUNICIPALITY_NAME, "snapshots": results}, ensure_ascii=False, indent=2))
    return 0


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