"""Provision the official DOV digital soil map for the municipality of Mol. The operator follows every bounded WFS page, retains checksummed source responses, clips soil polygons to the persisted Mol Area in EPSG:31370 and imports the result through GeoIntel's canonical dataset upload route. It does not write directly to vector_features and it does not treat the historical 1949-1971 field survey as a current drainage observation. """ from __future__ import annotations import argparse import hashlib import json import os import sys from collections import defaultdict from datetime import datetime, timezone from pathlib import Path from typing import Any, Iterable import requests from pyproj import Transformer from requests.adapters import HTTPAdapter from shapely.geometry import 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 WFS_URL = "https://www.dov.vlaanderen.be/geoserver/wfs" CATALOG_URL = ( "https://www.vlaanderen.be/datavindplaats/catalogus/" "digitale-bodemkaart-van-het-vlaams-gewest-bodemtypes" ) TYPE_NAME = "bodemkaart:bodemtypes" SOURCE_NAME = "dov_soil_map" SOURCE_VERSION = "Digitale uitgave juni 2017" SURVEY_PERIOD = "1949-1971" OBSERVED_AT = "1971-12-31T23:59:59Z" VALID_FROM = "1949-01-01T00:00:00Z" VALID_TO = OBSERVED_AT ATTRIBUTION = "Databank Ondergrond Vlaanderen - Digitale bodemkaart: bodemtypes" DEFAULT_API_URL = "http://127.0.0.1:8000" DEFAULT_OUTPUT_DIR = "/app/storage/operator-evidence/dov-soil-map/mol" DEFAULT_PROJECT_NAME = "Mol Municipality Workbench" DEFAULT_AREA_FRAGMENT = "Gemeente Mol" DATASET_FILENAME = "dov_soil_map_mol.geojson" MANIFEST_FILENAME = "dov_soil_map_mol.manifest.json" SCHEMA_VERSION = 1 TO_LAMBERT72 = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True) TO_WGS84 = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description="Provision the official DOV soil map for Mol.") 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-fragment", default=DEFAULT_AREA_FRAGMENT) parser.add_argument( "--output-dir", type=Path, default=Path(os.environ.get("GEOINTEL_SOIL_MAP_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)), ) parser.add_argument("--page-limit", type=int, default=500) parser.add_argument("--max-features", type=int, default=20_000) 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 utc_now() -> str: return datetime.now(timezone.utc).isoformat() def sha256_bytes(value: bytes) -> str: return hashlib.sha256(value).hexdigest() 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_bytes_atomic(path: Path, value: bytes) -> None: path.parent.mkdir(parents=True, exist_ok=True) temporary = path.with_suffix(path.suffix + ".tmp") temporary.write_bytes(value) temporary.replace(path) def write_json_atomic(path: Path, value: Any, *, pretty: bool = False) -> None: encoded = json.dumps( value, ensure_ascii=False, indent=2 if pretty else None, separators=None if pretty else (",", ":"), ).encode("utf-8") write_bytes_atomic(path, encoded) def source_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-DOV-Soil-Mol-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 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 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 = list(page.get("items") or []) page_total = int(page.get("total") or 0) if total is None: total = page_total elif page_total != total: raise RuntimeError("GeoIntel pagination total changed while locating the Mol workspace") 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 pagination returned {len(items)} of {total} items") return items def polygonal_geometry(geometry): if geometry is None or geometry.is_empty: return None if not geometry.is_valid: geometry = make_valid(geometry) polygons: list[Polygon] = [] def collect(candidate) -> None: if candidate is None or candidate.is_empty: return if isinstance(candidate, Polygon): polygons.append(candidate) elif isinstance(candidate, MultiPolygon): polygons.extend(part for part in candidate.geoms if not part.is_empty) elif hasattr(candidate, "geoms"): for part in candidate.geoms: collect(part) collect(geometry) if not polygons: return None result = unary_union(polygons) if not result.is_valid: result = make_valid(result) return result if not result.is_empty and result.is_valid else None def locate_workspace( session: requests.Session, base_url: str, project_name: str, area_fragment: str, timeout: int, ) -> tuple[str, str, Any, list[dict[str, Any]]]: projects = paginated_items(session, f"{base_url}/api/v1/projects", timeout=timeout) project = next((item for item in projects 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 = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", timeout=timeout) area = next( (item for item in areas if area_fragment.casefold() in str(item.get("name") or "").casefold()), None, ) if not area or not area.get("geometry"): raise RuntimeError(f"Persisted Mol Area containing {area_fragment!r} is missing") boundary_wgs84 = polygonal_geometry(shape(area["geometry"])) if boundary_wgs84 is None: raise RuntimeError("Persisted Mol Area is not valid polygonal geometry") datasets = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/datasets", timeout=timeout) return project_id, str(area["id"]), boundary_wgs84, datasets def iter_wfs_pages( session: requests.Session, bbox_lambert72: tuple[float, float, float, float], *, page_limit: int, timeout: int, ) -> Iterable[tuple[dict[str, Any], str, bytes]]: start_index = 0 expected_total: int | None = None while True: params = { "service": "WFS", "version": "2.0.0", "request": "GetFeature", "typeNames": TYPE_NAME, "srsName": "EPSG:4326", "bbox": ",".join(f"{value:.3f}" for value in bbox_lambert72) + ",EPSG:31370", "count": str(page_limit), "startIndex": str(start_index), "sortBy": "gid", "outputFormat": "application/json", } response = session.get(WFS_URL, params=params, timeout=timeout) response.raise_for_status() payload = response.json() if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection": raise RuntimeError("DOV WFS returned an invalid FeatureCollection") features = list(payload.get("features") or []) matched = int(payload.get("numberMatched") or payload.get("totalFeatures") or 0) if expected_total is None: expected_total = matched elif matched != expected_total: raise RuntimeError("DOV WFS numberMatched changed during pagination") yield payload, response.url, response.content returned = int(payload.get("numberReturned") or len(features)) if returned != len(features): raise RuntimeError("DOV WFS numberReturned does not match its feature payload") start_index += returned if returned == 0 or start_index >= expected_total: if start_index != expected_total: raise RuntimeError(f"DOV WFS returned {start_index} of {expected_total} matched features") break def normalize_feature( feature: dict[str, Any], boundary_lambert72, *, municipality: str = "Mol", nis_code: str = "13025", coverage_scope: str = "municipality", feature_id_suffix: str | None = None, ) -> tuple[dict[str, Any] | None, bool]: geometry_payload = feature.get("geometry") if not geometry_payload: return None, False source_wgs84 = polygonal_geometry(shape(geometry_payload)) if source_wgs84 is None: return None, False source_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, source_wgs84)) if source_lambert72 is None or not source_lambert72.intersects(boundary_lambert72): return None, False was_clipped = not source_lambert72.within(boundary_lambert72) clipped_lambert72 = polygonal_geometry(source_lambert72.intersection(boundary_lambert72)) if clipped_lambert72 is None or clipped_lambert72.area <= 0: return None, was_clipped clipped_wgs84 = polygonal_geometry(transform_geometry(TO_WGS84.transform, clipped_lambert72)) if clipped_wgs84 is None: return None, was_clipped raw = dict(feature.get("properties") or {}) gid = raw.get("gid") map_polygon_id = raw.get("id_kaartvlak") source_id = str(feature.get("id") or f"{TYPE_NAME}:{gid or map_polygon_id}") persisted_id = f"{source_id}:{feature_id_suffix}" if feature_id_suffix else source_id properties = { "source_name": SOURCE_NAME, "source_collection": TYPE_NAME, "source_feature_id": persisted_id, "source_gid": gid, "source_map_polygon_id": map_polygon_id, "reference_layer_name": "soil", "theme": "soil", "authority_level": "authoritative_historical_baseline", "coverage_scope": coverage_scope, "municipality": municipality, "nis_code": nis_code, "source_version": SOURCE_VERSION, "survey_period": SURVEY_PERIOD, "soil_type_code": raw.get("Bodemtype"), "unified_soil_type_code": raw.get("Unibodemtype"), "soil_series_code": raw.get("Bodemserie"), "soil_series_description": raw.get("Beknopte_omschrijving_bodemserie"), "soil_generalized_legend": raw.get("Gegeneraliseerde_legende"), "soil_texture_class_code": raw.get("Textuurklasse_code"), "soil_texture_class": raw.get("Textuurklasse"), "soil_drainage_class_code": raw.get("Drainageklasse_code"), "soil_drainage_class": raw.get("Drainageklasse"), "soil_profile_group_code": raw.get("Profielontwikkelingsgroep_code"), "soil_profile_group": raw.get("Profielontwikkelingsgroep"), "soil_substrate_code": raw.get("Substraat_code"), "soil_substrate": raw.get("Substraat_Vlaanderen") or raw.get("Substraat_legende"), "soil_region": raw.get("Streek"), "classification_type": raw.get("Type_classificatie"), "soil_map_title": raw.get("Eenduidige_legende_titel"), "clipped_area_ha": round(float(clipped_lambert72.area) / 10_000.0, 8), "attribution": ATTRIBUTION, "historical_drainage_limitation": ( "Drainage class derives from field data collected between 1949 and 1971 and may differ today." ), } return { "type": "Feature", "id": persisted_id, "geometry": mapping(clipped_wgs84), "properties": properties, }, was_clipped def prepare_artifact( session: requests.Session, boundary_wgs84, output_dir: Path, *, page_limit: int, max_features: int, timeout: int, ) -> tuple[Path, Path, dict[str, Any]]: output_dir.mkdir(parents=True, exist_ok=True) raw_dir = output_dir / "raw" raw_dir.mkdir(parents=True, exist_ok=True) boundary_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, boundary_wgs84)) if boundary_lambert72 is None: raise RuntimeError("Mol boundary could not be transformed to EPSG:31370") retained: list[dict[str, Any]] = [] raw_pages: list[dict[str, Any]] = [] source_urls: list[str] = [] seen_ids: set[str] = set() raw_feature_count = 0 duplicate_count = 0 rejected_count = 0 clipped_count = 0 area_by_legend: dict[str, float] = defaultdict(float) area_by_texture: dict[str, float] = defaultdict(float) area_by_drainage: dict[str, float] = defaultdict(float) for page_number, (payload, source_url, raw_bytes) in enumerate( iter_wfs_pages( session, boundary_lambert72.bounds, page_limit=page_limit, timeout=timeout, ), start=1, ): page_path = raw_dir / f"dov_soil_map_page_{page_number:05d}.json" write_bytes_atomic(page_path, raw_bytes) features = list(payload.get("features") or []) raw_feature_count += len(features) if raw_feature_count > max_features: raise RuntimeError( f"DOV WFS exceeded the {max_features} feature safety limit; refusing a truncated import" ) raw_pages.append( { "path": str(page_path.relative_to(output_dir)), "sha256": sha256_bytes(raw_bytes), "size_bytes": len(raw_bytes), "feature_count": len(features), "source_url": source_url, } ) source_urls.append(source_url) for feature in features: raw = dict(feature.get("properties") or {}) source_id = str(feature.get("id") or f"{TYPE_NAME}:{raw.get('gid')}") if source_id in seen_ids: duplicate_count += 1 continue seen_ids.add(source_id) normalized, was_clipped = normalize_feature(feature, boundary_lambert72) if normalized is None: rejected_count += 1 continue if was_clipped: clipped_count += 1 retained.append(normalized) properties = normalized["properties"] area = float(properties["clipped_area_ha"]) area_by_legend[str(properties.get("soil_generalized_legend") or "Onbekend")] += area area_by_texture[str(properties.get("soil_texture_class") or "Onbekend")] += area area_by_drainage[str(properties.get("soil_drainage_class") or "Onbekend")] += area if not retained: raise RuntimeError("DOV WFS returned no valid soil polygons inside the persisted Mol Area") generated_at = utc_now() artifact = { "type": "FeatureCollection", "name": "Digitale bodemkaart - Gemeente Mol", "features": retained, "source": ATTRIBUTION, "source_version": SOURCE_VERSION, "survey_period": SURVEY_PERIOD, "catalog_url": CATALOG_URL, "generated_at": generated_at, } artifact_path = output_dir / DATASET_FILENAME write_json_atomic(artifact_path, artifact) manifest = { "schema_version": SCHEMA_VERSION, "source_version": SOURCE_VERSION, "survey_period": SURVEY_PERIOD, "source_type_name": TYPE_NAME, "wfs_url": WFS_URL, "catalog_url": CATALOG_URL, "attribution": ATTRIBUTION, "generated_at": generated_at, "crs_source_service": "EPSG:31370", "crs_response_and_persisted": "EPSG:4326", "crs_clip_and_area_measurement": "EPSG:31370", "boundary_sha256": sha256_bytes(json.dumps(mapping(boundary_wgs84), sort_keys=True).encode("utf-8")), "boundary_bbox_wgs84": list(boundary_wgs84.bounds), "boundary_bbox_epsg31370": list(boundary_lambert72.bounds), "page_limit": page_limit, "page_count": len(raw_pages), "raw_source_feature_count": raw_feature_count, "feature_count": len(retained), "duplicate_count": duplicate_count, "rejected_or_outside_count": rejected_count, "clipped_feature_count": clipped_count, "reference_truncated": False, "raw_pages": raw_pages, "source_urls": source_urls, "area_by_generalized_legend_ha": {key: round(value, 6) for key, value in sorted(area_by_legend.items())}, "area_by_texture_ha": {key: round(value, 6) for key, value in sorted(area_by_texture.items())}, "area_by_drainage_ha": {key: round(value, 6) for key, value in sorted(area_by_drainage.items())}, "artifact_path": str(artifact_path), "artifact_sha256": sha256_file(artifact_path), "artifact_size_bytes": artifact_path.stat().st_size, "limitations": [ "The map is based on field data collected between 1949 and 1971.", "Current drainage, land use and local soil disturbance may differ from the mapped class.", "The 1:20,000 source is contextual evidence and not a parcel-scale soil investigation.", ], } manifest_path = output_dir / MANIFEST_FILENAME write_json_atomic(manifest_path, manifest, pretty=True) return artifact_path, manifest_path, manifest def selection_metrics() -> list[dict[str, Any]]: return [ { "metric_key": "soil_dry_sand_area", "method": "intersection_area", "label": "Gekarteerd als droog zand", "unit": "ha", "geometry_dimension": 2, "filter_property": "soil_generalized_legend", "filter_values": ["Droog zand", "Zeer droog zand"], }, { "metric_key": "soil_moist_sand_area", "method": "intersection_area", "label": "Gekarteerd als vochtig zand", "unit": "ha", "geometry_dimension": 2, "filter_property": "soil_generalized_legend", "filter_values": ["Vochtig zand"], }, { "metric_key": "soil_wet_sand_area", "method": "intersection_area", "label": "Gekarteerd als nat zand", "unit": "ha", "geometry_dimension": 2, "filter_property": "soil_generalized_legend", "filter_values": ["Nat zand", "Zeer nat zand"], }, { "metric_key": "soil_anthropogenic_area", "method": "intersection_area", "label": "Antropogene bodemklasse", "unit": "ha", "geometry_dimension": 2, "filter_property": "soil_generalized_legend", "filter_values": ["Antropogeen"], }, ] def upload_artifact( session: requests.Session, *, base_url: str, project_id: str, area_id: str, artifact_path: Path, manifest_path: Path, manifest: dict[str, Any], timeout: int, ) -> dict[str, Any]: limitation = ( "Historische bodemkartering op schaal 1:20.000 op basis van veldwerk 1949-1971; " "de huidige drainage en lokale bodemtoestand kunnen afwijken." ) source_metadata = { "provider": SOURCE_NAME, "theme": "soil", "layer_type": "soil", "source_collection": TYPE_NAME, "source_crs": "EPSG:31370", "persisted_crs": "EPSG:4326", "authority_level": "authoritative_historical_baseline", "coverage_scope": "municipality", "municipality": "Mol", "nis_code": "13025", "feature_count": manifest["feature_count"], "geometry_clipped_to_area": True, "semantic_metrics": False, "survey_period": SURVEY_PERIOD, "source_scale": "1:20,000", "attribution": ATTRIBUTION, "catalog_url": CATALOG_URL, "license_note": "DOV standard attribution and public GDI reuse conditions apply.", "limitation_message": limitation, "selection_aggregation": { "metric_key": "soil_mapped_area", "method": "intersection_area", "label": "Bodemkaartoppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": limitation, }, "selection_metrics": selection_metrics(), } provenance_metadata = { "operator_tool": "provision_mol_soil_map.py", "operator_explicit_fetch": True, "geometry_clipped_to_area": True, "source_type_name": TYPE_NAME, "wfs_url": WFS_URL, "catalog_url": CATALOG_URL, "manifest_path": str(manifest_path), "artifact_sha256": manifest["artifact_sha256"], "raw_page_checksums": {page["path"]: page["sha256"] for page in manifest["raw_pages"]}, "source_urls": manifest["source_urls"], "reference_truncated": False, "generated_at": manifest["generated_at"], "limitations": manifest["limitations"], } with artifact_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": SOURCE_NAME, "reference_layer_name": "soil", "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": "dov:digital-soil-map:mol", "observed_at": OBSERVED_AT, "valid_from": VALID_FROM, "valid_to": VALID_TO, "temporal_granularity": "period", "source_version": SOURCE_VERSION, }, files={"file": (artifact_path.name, handle, "application/geo+json")}, timeout=timeout, ) return response_data(response) def main() -> int: args = parse_args() if args.page_limit < 1 or args.page_limit > 2000 or args.max_features < args.page_limit: print(json.dumps({"status": "error", "message": "Invalid page or feature safety limits"}), file=sys.stderr) return 2 base_url = args.base_url.rstrip("/") api_session = requests.Session() try: project_id, area_id, boundary, datasets = locate_workspace( api_session, base_url, args.project_name, args.area_fragment, args.request_timeout, ) existing = next( ( item for item in datasets if item.get("source_name") == SOURCE_NAME and str(item.get("area_id") or "") == area_id and item.get("status") == "ready" ), None, ) if existing and not args.force: result = { "status": "reused", "project_id": project_id, "area_id": area_id, "dataset_id": existing["id"], "feature_count": existing.get("feature_count"), } else: artifact_path, manifest_path, manifest = prepare_artifact( source_session(), boundary, args.output_dir, page_limit=args.page_limit, max_features=args.max_features, timeout=args.request_timeout, ) if args.fetch_only: result = { "status": "prepared", "project_id": project_id, "area_id": area_id, "artifact_path": str(artifact_path), "feature_count": manifest["feature_count"], } else: dataset = upload_artifact( api_session, base_url=base_url, project_id=project_id, area_id=area_id, artifact_path=artifact_path, manifest_path=manifest_path, manifest=manifest, timeout=args.import_timeout, ) result = { "status": "imported", "project_id": project_id, "area_id": area_id, "dataset_id": dataset["id"], "feature_count": dataset.get("feature_count") or manifest["feature_count"], "artifact_path": str(artifact_path), } print(json.dumps(result, ensure_ascii=False, indent=2)) return 0 except (OSError, RuntimeError, requests.RequestException, ValueError) as exc: print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr) return 1 if __name__ == "__main__": raise SystemExit(main())