#!/usr/bin/env python3 """Provision the official Wallonia 1 m MNT for bounded GeoIntel analysis.""" from __future__ import annotations import argparse import hashlib import json from pathlib import Path import shutil import sys from urllib.request import Request, urlopen from zipfile import BadZipFile, ZipFile SOURCE_URL = ( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "fe13bc84-e371-46ca-9632-8ad4139f1ee5/RELIEF_WALLONIE_MNT_1M_2021_2022_GEOTIFF_3812.zip" ) TARGET_FILENAME = "spw_mnt_1m_2021_2022_3812.tif" ARCHIVE_FILENAME = "spw_mnt_1m_2021_2022_3812.zip" MIN_ARCHIVE_BYTES = 30_000_000_000 MAX_ARCHIVE_BYTES = 60_000_000_000 MAX_EXTRACTED_BYTES = 80_000_000_000 def sha256_file(path: Path) -> str: digest = hashlib.sha256() with path.open("rb") as handle: while chunk := handle.read(16 * 1024 * 1024): digest.update(chunk) return digest.hexdigest() def write_text_atomic(path: Path, content: str, *, encoding: str) -> None: temporary = path.with_suffix(path.suffix + ".part") temporary.unlink(missing_ok=True) try: temporary.write_text(content, encoding=encoding) temporary.replace(path) except Exception: temporary.unlink(missing_ok=True) raise def download(destination: Path) -> str: temporary = destination.with_suffix(destination.suffix + ".part") temporary.unlink(missing_ok=True) digest = hashlib.sha256() received = 0 request = Request( SOURCE_URL, headers={"User-Agent": "GeoIntel/1.0 SPW-terrain-source-provisioner"}, ) try: with urlopen(request, timeout=300) as response, temporary.open("wb") as output: content_length = int(response.headers.get("Content-Length") or 0) if ( content_length and not MIN_ARCHIVE_BYTES <= content_length <= MAX_ARCHIVE_BYTES ): raise RuntimeError( f"official archive size is outside the governed range: {content_length}" ) while chunk := response.read(16 * 1024 * 1024): received += len(chunk) if received > MAX_ARCHIVE_BYTES: raise RuntimeError( "official archive exceeds the governed 60 GB transfer limit" ) digest.update(chunk) output.write(chunk) if received % (1024 * 1024 * 1024) < len(chunk): print( f" downloaded {received / 1024 / 1024 / 1024:.1f} GiB", flush=True, ) if received < MIN_ARCHIVE_BYTES: raise RuntimeError( f"official archive is unexpectedly small: {received} bytes" ) temporary.replace(destination) return digest.hexdigest() except Exception: temporary.unlink(missing_ok=True) raise def extract_single_geotiff(archive: Path, target: Path) -> None: try: with ZipFile(archive) as bundle: candidates = [ item for item in bundle.infolist() if not item.is_dir() and item.filename.lower().endswith((".tif", ".tiff")) ] if len(candidates) != 1: raise RuntimeError( f"archive must contain exactly one GeoTIFF, found {len(candidates)}" ) member = candidates[0] if member.file_size <= 0 or member.file_size > MAX_EXTRACTED_BYTES: raise RuntimeError( f"GeoTIFF uncompressed size is outside the governed limit: {member.file_size}" ) temporary = target.with_suffix(target.suffix + ".part") temporary.unlink(missing_ok=True) with bundle.open(member) as source, temporary.open("wb") as output: shutil.copyfileobj(source, output, length=16 * 1024 * 1024) temporary.replace(target) except BadZipFile as exc: raise RuntimeError("official SPW MNT archive is not a valid ZIP file") from exc def validate_raster(path: Path) -> dict: try: import numpy as np import rasterio from rasterio.windows import Window except ImportError as exc: raise RuntimeError( "rasterio and numpy are required to validate the SPW MNT source" ) from exc with rasterio.open(path) as source: if source.crs is None or source.crs.to_epsg() != 3812: raise RuntimeError(f"SPW MNT must use EPSG:3812, found {source.crs}") if source.count != 1: raise RuntimeError(f"SPW MNT must have one band, found {source.count}") if not all(abs(abs(float(value)) - 1.0) <= 0.05 for value in source.res): raise RuntimeError(f"SPW MNT must retain 1 m cells, found {source.res}") sample_windows = [] sample_size = 512 for x_fraction, y_fraction in ( (0.1, 0.1), (0.5, 0.5), (0.9, 0.9), (0.1, 0.9), (0.9, 0.1), ): col = max( 0, min( source.width - sample_size, round(source.width * x_fraction - sample_size / 2), ), ) row = max( 0, min( source.height - sample_size, round(source.height * y_fraction - sample_size / 2), ), ) sample_windows.append( Window( col, row, min(sample_size, source.width), min(sample_size, source.height), ) ) samples = [ source.read(1, window=window, masked=True).compressed().astype("float64") for window in sample_windows ] values = np.concatenate([sample for sample in samples if sample.size]) values = values[np.isfinite(values)] if not values.size: raise RuntimeError( "SPW MNT validation samples contain no finite elevation values" ) if float(values.min()) < -100.0 or float(values.max()) > 1000.0: raise RuntimeError( f"SPW MNT samples contain implausible values: {values.min()}..{values.max()}" ) return { "path": str(path), "crs": str(source.crs), "width": int(source.width), "height": int(source.height), "resolution": [float(value) for value in source.res], "bounds": [float(value) for value in source.bounds], "nodata": None if source.nodata is None else float(source.nodata), "dtype": source.dtypes[0], "sample_min_m": float(values.min()), "sample_max_m": float(values.max()), } def provision(destination: Path, force: bool, keep_archive: bool) -> dict: destination.mkdir(parents=True, exist_ok=True) target = destination / TARGET_FILENAME archive = destination / ARCHIVE_FILENAME archive_digest = None if target.is_file() and not force: print(f"SPW MNT: validating existing source {target}") else: if archive.is_file(): archive_size = archive.stat().st_size if not MIN_ARCHIVE_BYTES <= archive_size <= MAX_ARCHIVE_BYTES: raise RuntimeError( f"existing archive size is outside the governed range: {archive_size}" ) print(f"SPW MNT: using existing archive {archive}") archive_digest = sha256_file(archive) else: print("SPW MNT: downloading official archive") archive_digest = download(archive) print(f"SPW MNT: archive sha256 {archive_digest}") extract_single_geotiff(archive, target) if not keep_archive: archive.unlink(missing_ok=True) validation = validate_raster(target) source_digest = sha256_file(target) write_text_atomic( target.with_suffix(".sha256"), f"{source_digest} {target.name}\n", encoding="ascii", ) validation.update( { "source_sha256": source_digest, "archive_sha256": archive_digest, "download_url": SOURCE_URL, "catalog_url": "https://geoportail.wallonie.be/catalogue/fe13bc84-e371-46ca-9632-8ad4139f1ee5.html", } ) report_path = destination / "provisioning-report.json" write_text_atomic( report_path, json.dumps(validation, indent=2) + "\n", encoding="utf-8", ) print(f"SPW MNT: ready ({target.stat().st_size / 1024 / 1024 / 1024:.1f} GiB)") print(f"Provisioning report: {report_path}") return validation def main() -> int: parser = argparse.ArgumentParser( description="Provision the official Wallonia 1 m MNT 2021-2022 GeoTIFF." ) parser.add_argument( "--destination", type=Path, default=Path("storage/source-cache/spw-terrain") ) parser.add_argument("--force", action="store_true") parser.add_argument("--keep-archive", action="store_true") args = parser.parse_args() provision(args.destination.resolve(), args.force, args.keep_archive) return 0 if __name__ == "__main__": try: raise SystemExit(main()) except Exception as exc: print(f"SPW_TERRAIN_PROVISIONING_FAILED: {exc}", file=sys.stderr) raise SystemExit(1) from exc