"""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())