Files
geointel/scripts/import_spw_bathymetry.py
T
Codex 70c5e34ecf
GeoIntel release gates / Compile, test, contracts and builds (push) Canceled after 0s
GeoIntel release gates / Python and npm vulnerability policy (push) Canceled after 0s
GeoIntel release gates / GIS image, SBOM and container scan (push) Canceled after 0s
Add governed SPW bathymetry analysis
2026-07-19 04:07:01 +02:00

630 lines
24 KiB
Python

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