Files
geointel/scripts/provision_mol_historical_landuse.py
T
Codex c0943fe7d4
GeoIntel CI / docs-smoke (push) Canceled after 0s
GeoIntel CI / contract-smoke (push) Canceled after 0s
feat: add temporal Mol explorer
2026-07-14 15:19:42 +02:00

553 lines
24 KiB
Python

"""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"}),
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 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:
response = session.get(
WFS_URL,
params={
"service": "WFS",
"version": "2.0.0",
"request": "GetFeature",
"typeNames": collection,
"outputFormat": "application/json",
"srsName": "EPSG:4326",
"FILTER": filter_xml,
"count": page_size,
"startIndex": start_index,
},
timeout=timeout,
)
response.raise_for_status()
page = response.json().get("features") or []
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:
response = session.get(
WFS_URL,
params={
"service": "WFS",
"version": "2.0.0",
"request": "GetFeature",
"typeNames": collection,
"outputFormat": "application/json",
"srsName": "EPSG:4326",
"FILTER": filter_xml,
"count": page_size,
"startIndex": start_index,
},
timeout=timeout,
)
response.raise_for_status()
page = response.json().get("features") or []
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())