"""Validate, crop and import the official SPW bathymetry release. The operator is intentionally bounded and fail-closed. It accepts one pinned official archive, keeps the source checksum in provenance, creates a compact COG for a requested EPSG:4326 scope and persists it through DatasetService's canonical upload API. It never writes directly to PostGIS or dataset storage. """ from __future__ import annotations import argparse from datetime import UTC, datetime from hashlib import sha256 import json import math import mimetypes import os from pathlib import Path, PurePosixPath import tempfile from typing import Any from urllib.error import HTTPError, URLError from urllib.request import Request, urlopen import uuid import zipfile from pyproj import Transformer import rasterio from rasterio.features import geometry_mask from rasterio.mask import mask from rasterio.shutil import copy as raster_copy from shapely.geometry import box, mapping, shape from shapely.ops import transform as shapely_transform SOURCE_URL = ( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "0a544b42-0b30-4c8e-85e7-38149b99eae0/" "BATHY_50CM_ALTITUDE_DNG_GEOTIFF_3812.zip" ) CATALOG_URL = ( "https://geoportail.wallonie.be/catalogue/" "0a544b42-0b30-4c8e-85e7-38149b99eae0.html" ) SOURCE_SHA256 = "04122a1c5cecc7b77be025b580995d024545e82d8310e158b4c112f7fa5f8a6e" SOURCE_MEMBER = "BATHY_50CM_ALTITUDE_DNG.tif" SOURCE_VERSION = "SPW bathymetry 2023-05-23" SURVEY_PERIOD = "2019-2022" SOURCE_CRS = "EPSG:3812" VERTICAL_REFERENCE = "mDNG" NODATA = -9999.0 MAX_ARCHIVE_BYTES = 600 * 1024 * 1024 MAX_UNCOMPRESSED_BYTES = 10 * 1024 * 1024 * 1024 MAX_ARCHIVE_MEMBERS = 32 MAX_COMPRESSION_RATIO = 100.0 DEFAULT_MAX_PIXELS = 30_000_000 DEFAULT_PROJECT_NAME = "Belgium and North Sea Workbench" DEFAULT_API_URL = "http://127.0.0.1:8000" class SpwBathymetryImportError(RuntimeError): pass def sha256_file(path: Path) -> str: digest = sha256() with path.open("rb") as stream: for chunk in iter(lambda: stream.read(8 * 1024 * 1024), b""): digest.update(chunk) return digest.hexdigest() def download_source(path: Path, *, timeout: int) -> None: request = Request( SOURCE_URL, headers={ "Accept": "application/zip", "User-Agent": "GeoIntel/1.0 governed-spw-bathymetry-operator", }, ) path.parent.mkdir(parents=True, exist_ok=True) temporary = path.with_suffix(path.suffix + ".partial") digest = sha256() size = 0 try: with urlopen(request, timeout=timeout) as response, temporary.open("wb") as output: if "zip" not in str(response.headers.get("Content-Type") or "").lower(): raise SpwBathymetryImportError("SPW response is not an official ZIP artifact") for chunk in iter(lambda: response.read(8 * 1024 * 1024), b""): size += len(chunk) if size > MAX_ARCHIVE_BYTES: raise SpwBathymetryImportError("SPW source archive exceeds the configured size limit") output.write(chunk) digest.update(chunk) except (HTTPError, URLError, TimeoutError, OSError) as exc: temporary.unlink(missing_ok=True) raise SpwBathymetryImportError(f"Official SPW source download failed: {exc}") from exc if digest.hexdigest() != SOURCE_SHA256: temporary.unlink(missing_ok=True) raise SpwBathymetryImportError("Downloaded SPW artifact checksum does not match the pinned release") temporary.replace(path) def validate_archive(path: Path) -> zipfile.ZipInfo: if not path.is_file(): raise SpwBathymetryImportError(f"SPW source archive does not exist: {path}") if path.stat().st_size > MAX_ARCHIVE_BYTES: raise SpwBathymetryImportError("SPW source archive exceeds the configured size limit") actual_sha256 = sha256_file(path) if actual_sha256 != SOURCE_SHA256: raise SpwBathymetryImportError( f"SPW source checksum mismatch: expected {SOURCE_SHA256}, received {actual_sha256}" ) try: with zipfile.ZipFile(path) as archive: members = archive.infolist() if len(members) > MAX_ARCHIVE_MEMBERS: raise SpwBathymetryImportError("SPW archive contains too many members") total_size = 0 selected: zipfile.ZipInfo | None = None for member in members: member_path = PurePosixPath(member.filename) if member.is_dir(): continue if member_path.is_absolute() or ".." in member_path.parts: raise SpwBathymetryImportError("SPW archive contains an unsafe member path") total_size += member.file_size compressed = max(1, member.compress_size) if member.file_size / compressed > MAX_COMPRESSION_RATIO: raise SpwBathymetryImportError("SPW archive contains an unsafe compression ratio") if member_path.name == SOURCE_MEMBER: selected = member if total_size > MAX_UNCOMPRESSED_BYTES: raise SpwBathymetryImportError("SPW archive expands beyond the configured size limit") if selected is None: raise SpwBathymetryImportError(f"SPW archive does not contain {SOURCE_MEMBER}") return selected except zipfile.BadZipFile as exc: raise SpwBathymetryImportError("SPW source artifact is not a valid ZIP archive") from exc def parse_bbox(value: str): try: coordinates = [float(item.strip()) for item in value.split(",")] except ValueError as exc: raise argparse.ArgumentTypeError("bbox must contain four numeric EPSG:4326 coordinates") from exc if len(coordinates) != 4: raise argparse.ArgumentTypeError("bbox must contain min_x,min_y,max_x,max_y") min_x, min_y, max_x, max_y = coordinates if ( not all(math.isfinite(item) for item in coordinates) or min_x >= max_x or min_y >= max_y or min_x < -180 or max_x > 180 or min_y < -90 or max_y > 90 ): raise argparse.ArgumentTypeError("bbox is not a valid EPSG:4326 extent") return box(min_x, min_y, max_x, max_y) def unwrap(payload: Any, status_code: int) -> Any: if not isinstance(payload, dict): raise SpwBathymetryImportError(f"GeoIntel returned an invalid JSON envelope (HTTP {status_code})") if status_code >= 400 or payload.get("error"): error = payload.get("error") if isinstance(payload.get("error"), dict) else {} message = error.get("message") or f"GeoIntel request failed with HTTP {status_code}" raise SpwBathymetryImportError(str(message)) if "data" not in payload: raise SpwBathymetryImportError("GeoIntel response does not use the canonical data envelope") return payload["data"] def api_json( base_url: str, path: str, *, timeout: int, method: str = "GET", payload: dict[str, Any] | None = None, ) -> Any: content = json.dumps(payload).encode("utf-8") if payload is not None else None request = Request( f"{base_url.rstrip('/')}{path}", data=content, method=method, headers={ "Accept": "application/json", "Content-Type": "application/json", "User-Agent": "GeoIntel/1.0 governed-spw-bathymetry-operator", }, ) try: with urlopen(request, timeout=timeout) as response: status_code = int(response.status) response_content = response.read() except HTTPError as exc: status_code = exc.code response_content = exc.read() except (URLError, TimeoutError, OSError) as exc: raise SpwBathymetryImportError(f"GeoIntel API is unavailable: {exc}") from exc try: response_payload = json.loads(response_content.decode("utf-8")) except (UnicodeDecodeError, json.JSONDecodeError) as exc: raise SpwBathymetryImportError(f"GeoIntel returned non-JSON HTTP {status_code}") from exc return unwrap(response_payload, status_code) def paged_items(base_url: str, path: str, *, timeout: int) -> list[dict[str, Any]]: items: list[dict[str, Any]] = [] offset = 0 while True: separator = "&" if "?" in path else "?" page = api_json( base_url, f"{path}{separator}limit=200&offset={offset}", timeout=timeout, ) rows = list(page.get("items") or []) items.extend(item for item in rows if isinstance(item, dict)) total = int(page.get("total") or 0) if not rows or len(items) >= total: return items offset += len(rows) def resolve_project( base_url: str, project_name: str, *, project_id: str | None, timeout: int, ) -> dict[str, Any]: if project_id: return api_json(base_url, f"/api/v1/projects/{project_id}", timeout=timeout) projects = paged_items(base_url, "/api/v1/projects?status=active", timeout=timeout) matches = [item for item in projects if str(item.get("name") or "").casefold() == project_name.casefold()] if len(matches) != 1: raise SpwBathymetryImportError( f"Expected one active project named {project_name!r}, found {len(matches)}" ) return matches[0] def resolve_area( base_url: str, project_id: str, *, area_id: str | None, area_name: str | None, timeout: int, ) -> dict[str, Any] | None: if not area_id and not area_name: return None areas = paged_items(base_url, f"/api/v1/projects/{project_id}/areas", timeout=timeout) if area_id: matches = [item for item in areas if str(item.get("id")) == area_id] else: matches = [ item for item in areas if str(area_name).casefold() in str(item.get("name") or "").casefold() ] if len(matches) != 1: raise SpwBathymetryImportError(f"Expected exactly one matching persisted Area, found {len(matches)}") if not isinstance(matches[0].get("geometry"), dict): raise SpwBathymetryImportError("Selected Area has no persisted geometry") return matches[0] def source_path(archive_path: Path, member: zipfile.ZipInfo) -> str: return f"/vsizip/{archive_path.resolve().as_posix()}/{member.filename}" def selection_hash(selection) -> str: payload = json.dumps(mapping(selection), sort_keys=True, separators=(",", ":")) return sha256(payload.encode("utf-8")).hexdigest() def crop_source( archive_path: Path, member: zipfile.ZipInfo, selection_4326, output_path: Path, *, max_pixels: int, ) -> dict[str, Any]: try: with rasterio.open(source_path(archive_path, member)) as source: if source.crs is None or source.crs.to_epsg() != 3812: raise SpwBathymetryImportError("SPW bathymetry source CRS must be EPSG:3812") if source.count != 1 or source.dtypes[0] != "float32": raise SpwBathymetryImportError("SPW bathymetry source must contain one Float32 band") if source.nodata is None or not math.isclose(float(source.nodata), NODATA): raise SpwBathymetryImportError("SPW bathymetry source nodata value must be -9999") if not ( math.isclose(abs(source.res[0]), 0.5, rel_tol=0.01) and math.isclose(abs(source.res[1]), 0.5, rel_tol=0.01) ): raise SpwBathymetryImportError("SPW bathymetry source resolution must be approximately 0.5 m") transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True) selection_metric = shapely_transform(transformer.transform, selection_4326) selection_metric = selection_metric.intersection(box(*source.bounds)) if selection_metric.is_empty or selection_metric.area <= 0: raise SpwBathymetryImportError("Requested scope does not overlap the official SPW bathymetry") min_x, min_y, max_x, max_y = selection_metric.bounds expected_pixels = math.ceil((max_x - min_x) / abs(source.res[0])) * math.ceil( (max_y - min_y) / abs(source.res[1]) ) if expected_pixels > max_pixels: raise SpwBathymetryImportError( f"Requested scope requires {expected_pixels:,} cells; limit is {max_pixels:,}" ) clipped, transform = mask( source, [mapping(selection_metric)], crop=True, filled=True, nodata=NODATA, indexes=[1], ) profile = source.profile.copy() profile.update( driver="GTiff", width=clipped.shape[2], height=clipped.shape[1], count=1, transform=transform, nodata=NODATA, compress="DEFLATE", predictor=3, tiled=True, blockxsize=512, blockysize=512, BIGTIFF="IF_SAFER", ) except SpwBathymetryImportError: raise except Exception as exc: raise SpwBathymetryImportError(f"SPW bathymetry raster could not be read: {exc}") from exc output_path.parent.mkdir(parents=True, exist_ok=True) with tempfile.NamedTemporaryFile(suffix=".tif", dir=output_path.parent, delete=False) as temporary: temporary_path = Path(temporary.name) try: with rasterio.open(temporary_path, "w", **profile) as output: output.write(clipped) raster_copy( temporary_path, output_path, driver="COG", compress="DEFLATE", predictor=3, blocksize=512, overview_resampling="average", ) finally: temporary_path.unlink(missing_ok=True) with rasterio.open(output_path) as output: values = output.read(1, masked=True) valid = geometry_mask( [mapping(selection_metric)], out_shape=values.shape, transform=output.transform, invert=True, ) & ~values.mask valid_values = values.data[valid] if valid_values.size == 0: output_path.unlink(missing_ok=True) raise SpwBathymetryImportError("Selected scope contains no surveyed SPW waterbed cells") if not bool(((valid_values > -500.0) & (valid_values < 1000.0)).all()): output_path.unlink(missing_ok=True) raise SpwBathymetryImportError("SPW bathymetry contains values outside the governed mDNG range") return { "width": output.width, "height": output.height, "resolution_m": max(abs(output.res[0]), abs(output.res[1])), "valid_cell_count": int(valid_values.size), "valid_value_min": float(valid_values.min()), "valid_value_max": float(valid_values.max()), "valid_value_mean": float(valid_values.mean()), "output_sha256": sha256_file(output_path), } def multipart_upload( base_url: str, project_id: str, output_path: Path, fields: dict[str, str], *, timeout: int, ) -> dict[str, Any]: boundary = f"geointel-{uuid.uuid4().hex}" body = bytearray() for name, value in fields.items(): body.extend(f"--{boundary}\r\n".encode()) body.extend(f'Content-Disposition: form-data; name="{name}"\r\n\r\n'.encode()) body.extend(value.encode("utf-8")) body.extend(b"\r\n") body.extend(f"--{boundary}\r\n".encode()) body.extend( ( f'Content-Disposition: form-data; name="file"; filename="{output_path.name}"\r\n' f"Content-Type: {mimetypes.guess_type(output_path.name)[0] or 'image/tiff'}\r\n\r\n" ).encode() ) body.extend(output_path.read_bytes()) body.extend(f"\r\n--{boundary}--\r\n".encode()) request = Request( f"{base_url.rstrip('/')}/api/v1/projects/{project_id}/datasets/upload", data=bytes(body), method="POST", headers={ "Accept": "application/json", "Content-Type": f"multipart/form-data; boundary={boundary}", "User-Agent": "GeoIntel/1.0 governed-spw-bathymetry-operator", }, ) try: with urlopen(request, timeout=timeout) as response: status_code = int(response.status) response_content = response.read() except HTTPError as exc: status_code = exc.code response_content = exc.read() except (URLError, TimeoutError, OSError) as exc: raise SpwBathymetryImportError(f"GeoIntel upload failed: {exc}") from exc try: response_payload = json.loads(response_content.decode("utf-8")) except (UnicodeDecodeError, json.JSONDecodeError) as exc: raise SpwBathymetryImportError(f"GeoIntel upload returned non-JSON HTTP {status_code}") from exc return unwrap(response_payload, status_code) def parse_args() -> argparse.Namespace: parser = argparse.ArgumentParser(description="Import a bounded official SPW bathymetry raster.") 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("--project-id") parser.add_argument("--area") parser.add_argument("--area-id") parser.add_argument("--bbox", required=True, type=parse_bbox, help="min_x,min_y,max_x,max_y in EPSG:4326") parser.add_argument("--raw-zip", required=True, type=Path) parser.add_argument("--download", action="store_true", help="Download the pinned official ZIP when --raw-zip is absent.") parser.add_argument("--output-dir", type=Path, default=Path("storage/operator-evidence/spw-bathymetry/derived")) parser.add_argument("--max-pixels", type=int, default=DEFAULT_MAX_PIXELS) parser.add_argument("--timeout", type=int, default=900) parser.add_argument("--force-refresh", action="store_true") parser.add_argument("--dry-run", action="store_true") return parser.parse_args() def main() -> int: args = parse_args() archive_path = args.raw_zip.resolve() if not archive_path.exists() and args.download: download_source(archive_path, timeout=args.timeout) member = validate_archive(archive_path) project = resolve_project( args.base_url, args.project_name, project_id=args.project_id, timeout=args.timeout, ) project_id = str(project["id"]) area = resolve_area( args.base_url, project_id, area_id=args.area_id, area_name=args.area, timeout=args.timeout, ) selection = args.bbox if area is not None: selection = selection.intersection(shape(area["geometry"])) if selection.is_empty or selection.area <= 0: raise SpwBathymetryImportError("Requested bbox does not overlap the selected persisted Area") scope_hash = selection_hash(selection) output_name = f"spw_bathymetry_{scope_hash[:16]}_3812.tif" output_path = args.output_dir.resolve() / output_name raw_sha256 = sha256_file(archive_path) datasets = paged_items( args.base_url, f"/api/v1/projects/{project_id}/datasets", timeout=args.timeout, ) existing = next( ( item for item in datasets if item.get("status") == "ready" and item.get("source_name") == "spw_bathymetry" and (item.get("source_metadata") or {}).get("source_artifact_sha256") == raw_sha256 and (item.get("source_metadata") or {}).get("selection_hash") == scope_hash ), None, ) if existing and not args.force_refresh: print(json.dumps({"status": "reused", "dataset_id": existing["id"]}, ensure_ascii=False)) return 0 if args.dry_run: print( json.dumps( { "status": "validated", "project_id": project_id, "area_id": area.get("id") if area else None, "bbox_epsg4326": list(selection.bounds), "source_sha256": raw_sha256, "source_member": member.filename, "output_path": str(output_path), }, ensure_ascii=False, ) ) return 0 diagnostics = crop_source( archive_path, member, selection, output_path, max_pixels=args.max_pixels, ) source_metadata = { "product_key": "spw_bathymetry_50cm_mdng", "product_display_name": "SPW waterbodemhoogte 0,5 m", "theme": "bathymetry", "value_semantics": "bed_elevation", "vertical_reference": VERTICAL_REFERENCE, "source_crs": SOURCE_CRS, "native_resolution_m": 0.5, "analysis_resolution_m": diagnostics["resolution_m"], "nodata_value": NODATA, "survey_period": SURVEY_PERIOD, "published_on": "2023-05-23", "bbox_epsg4326": list(selection.bounds), "coverage_zones": ["wallonia"], "coverage_scope": "bounded_selection", "selection_hash": scope_hash, "source_artifact_sha256": raw_sha256, "valid_cell_count": diagnostics["valid_cell_count"], "valid_value_min": diagnostics["valid_value_min"], "valid_value_max": diagnostics["valid_value_max"], "valid_value_mean": diagnostics["valid_value_mean"], "attribution": "Service public de Wallonie (SPW)", "license_note": "CC BY 4.0", "catalog_url": CATALOG_URL, "limitation_message": ( "Waterbodemhoogte in mDNG uit opmetingen 2019-2022; geen actuele waterdiepte of watervolume." ), } provenance_metadata = { "operator": "import_spw_bathymetry.py", "official_source_url": SOURCE_URL, "official_catalog_url": CATALOG_URL, "source_artifact_sha256": raw_sha256, "source_archive_member": member.filename, "derived_artifact_sha256": diagnostics["output_sha256"], "selection_geometry_epsg4326": mapping(selection), "selection_hash": scope_hash, "processed_at": datetime.now(UTC).isoformat(), "water_depth_available": False, "water_volume_available": False, "vertical_datum_conversion_applied": False, } fields = { "dataset_type": "raster", "source": "SPW official operator archive", "dataset_role": "source", "source_name": "spw_bathymetry", "source_metadata_json": json.dumps(source_metadata, separators=(",", ":")), "provenance_metadata_json": json.dumps(provenance_metadata, separators=(",", ":")), "valid_from": "2019-01-01T00:00:00Z", "valid_to": "2022-12-31T23:59:59Z", "temporal_granularity": "period", "source_version": SOURCE_VERSION, } if area is not None: fields["area_id"] = str(area["id"]) created = multipart_upload( args.base_url, project_id, output_path, fields, timeout=args.timeout, ) analysis = api_json( args.base_url, f"/api/v1/projects/{project_id}/datasets/{created['id']}/raster/bathymetry/select", timeout=args.timeout, method="POST", payload={ "bbox": { "min_x": selection.bounds[0], "min_y": selection.bounds[1], "max_x": selection.bounds[2], "max_y": selection.bounds[3], "crs": "EPSG:4326", }, "area_id": area.get("id") if area else None, }, ) print( json.dumps( { "status": "imported", "dataset_id": created["id"], "output_path": str(output_path), "source_sha256": raw_sha256, "derived_sha256": diagnostics["output_sha256"], "summary": analysis["summary"], "coverage_ratio": analysis["coverage_ratio"], }, ensure_ascii=False, ) ) return 0 if __name__ == "__main__": try: raise SystemExit(main()) except SpwBathymetryImportError as exc: raise SystemExit(f"SPW_BATHYMETRY_IMPORT_FAILED: {exc}") from exc