geointel/scripts/provision_spw_terrain_source.py

#!/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