from __future__ import annotations from datetime import UTC, datetime import hashlib import json import math from pathlib import Path from typing import Any from uuid import UUID from geoalchemy2.shape import to_shape from pyproj import Transformer from shapely.geometry import box, mapping from shapely.ops import transform as shapely_transform from app.core.config import Settings, get_settings from app.core.errors import AppError from app.models import Area, Dataset, Project from app.schemas.spw_terrain import ( SpwTerrainAcquireRequest, SpwTerrainAcquisitionResult, SpwTerrainProductRead, ) from app.services.dataset_service import DatasetService class SpwTerrainService: PROVIDER = "spw_terrain" PRODUCT_KEY = "spw_mnt_1m_2021_2022" DISPLAY_NAME = "SPW terreinmodel (MNT) 2021-2022" SOURCE_FILENAME = "spw_mnt_1m_2021_2022_3812.tif" SOURCE_SHA256_FILENAME = "spw_mnt_1m_2021_2022_3812.sha256" SOURCE_CRS = "EPSG:3812" SOURCE_RESOLUTION_M = 1.0 SURFACE_MODEL = "terrain" VERTICAL_REFERENCE = "DNG / Deuxieme Nivellement General (EPSG:5710)" VERTICAL_UNIT_LABEL = "m DNG" ACQUISITION_PERIOD = "2021-02-19/2022-03-05" CATALOG_URL = "https://geoportail.wallonie.be/catalogue/fe13bc84-e371-46ca-9632-8ad4139f1ee5.html" DOWNLOAD_URL = ( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "fe13bc84-e371-46ca-9632-8ad4139f1ee5/RELIEF_WALLONIE_MNT_1M_2021_2022_GEOTIFF_3812.zip" ) ATTRIBUTION = ( "Service public de Wallonie (SPW) - Relief de la Wallonie MNT 2021-2022" ) LICENSE_NOTE = "CC BY 4.0; cite SPW and identify modifications." NODATA = -9999.0 LIMITATION = ( "GeoIntel leest uitsluitend een begrensd venster uit het checksum-gevalideerde officiele 1 m MNT en " "bewaart een analyse-afgeleide op de gekozen resolutie. Het MNT beschrijft maaiveldhoogte in DNG, niet " "oppervlaktehoogte, afstroming, waterdiepte of watervolume. Kleine bronzones zijn door SPW geinterpoleerd." ) @staticmethod def _source_path(settings: Settings) -> Path: return Path(settings.spw_terrain_source_dir) / SpwTerrainService.SOURCE_FILENAME @staticmethod def _source_sha256(settings: Settings) -> str | None: checksum_path = ( Path(settings.spw_terrain_source_dir) / SpwTerrainService.SOURCE_SHA256_FILENAME ) if not checksum_path.is_file(): return None parts = checksum_path.read_text(encoding="ascii").strip().split() digest = parts[0].lower() if parts else "" if len(digest) != 64 or any(character not in "0123456789abcdef" for character in digest): return None return digest @staticmethod def list_products(*, settings: Settings | None = None) -> list[dict[str, Any]]: resolved = settings or get_settings() configured = ( resolved.spw_terrain_enabled and SpwTerrainService._source_path(resolved).is_file() and SpwTerrainService._source_sha256(resolved) is not None ) product = SpwTerrainProductRead( key=SpwTerrainService.PRODUCT_KEY, display_name=SpwTerrainService.DISPLAY_NAME, surface_model=SpwTerrainService.SURFACE_MODEL, source_filename=SpwTerrainService.SOURCE_FILENAME, native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M, analysis_resolution_m=resolved.spw_terrain_analysis_resolution_m, source_crs=SpwTerrainService.SOURCE_CRS, vertical_reference=SpwTerrainService.VERTICAL_REFERENCE, acquisition_period=SpwTerrainService.ACQUISITION_PERIOD, catalog_url=SpwTerrainService.CATALOG_URL, attribution=SpwTerrainService.ATTRIBUTION, license_note=SpwTerrainService.LICENSE_NOTE, limitation_message=SpwTerrainService.LIMITATION, coverage_zones=["wallonia"], configured=configured, status="configured" if configured else "source_not_provisioned", ) return [product.model_dump()] @staticmethod def _scope_geometry(db, project_id: UUID, payload: SpwTerrainAcquireRequest): if not db.get(Project, project_id): raise AppError( code="PROJECT_NOT_FOUND", message="Project not found", status_code=404 ) if payload.product_key.strip().lower() != SpwTerrainService.PRODUCT_KEY: raise AppError( code="SPW_TERRAIN_PRODUCT_NOT_SUPPORTED", message="Select the governed SPW MNT 2021-2022 product", details={"product_key": payload.product_key}, status_code=422, ) if payload.bbox.crs.upper() != "EPSG:4326": raise AppError( code="INVALID_BBOX_CRS", message="SPW terrain acquisition requires EPSG:4326", status_code=400, ) values = [ payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y, ] if ( not all(math.isfinite(value) for value in values) or values[0] >= values[2] or values[1] >= values[3] ): raise AppError( code="INVALID_BBOX", message="SPW terrain selection must be a finite non-empty rectangle", status_code=400, ) selection = box(*values) if payload.area_id is None: return selection, values area = db.get(Area, payload.area_id) if area is None or area.project_id != project_id: raise AppError( code="AREA_NOT_FOUND", message="Area not found", status_code=404 ) selection = selection.intersection(to_shape(area.geometry)) if selection.is_empty or selection.area <= 0: raise AppError( code="SPW_TERRAIN_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422, ) return selection, values @staticmethod def _read_source_window( source_path: Path, scope_4326, resolution: float, settings: Settings ) -> tuple[bytes, dict[str, Any]]: try: import numpy as np import rasterio from rasterio.enums import Resampling from rasterio.features import geometry_mask from rasterio.io import MemoryFile from rasterio.transform import from_bounds from rasterio.windows import from_bounds as window_from_bounds except ImportError as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for SPW terrain", status_code=503, ) from exc scope_metric = shapely_transform( Transformer.from_crs( "EPSG:4326", SpwTerrainService.SOURCE_CRS, always_xy=True ).transform, scope_4326, ) try: with rasterio.open(source_path) as source: if ( source.crs is None or source.crs.to_epsg() != 3812 or source.count != 1 ): raise AppError( code="SPW_TERRAIN_SOURCE_INVALID", message="SPW MNT must be a one-band EPSG:3812 raster", status_code=409, ) if not all( math.isclose(abs(float(value)), 1.0, abs_tol=0.05) for value in source.res ): raise AppError( code="SPW_TERRAIN_SOURCE_INVALID", message="SPW MNT must retain the official 1 m resolution", status_code=409, ) clipped_geometry = scope_metric.intersection(box(*source.bounds)) if clipped_geometry.is_empty or clipped_geometry.area <= 0: raise AppError( code="SPW_TERRAIN_SELECTION_OUTSIDE_COVERAGE", message="Selection does not overlap SPW MNT coverage", status_code=422, ) min_x, min_y, max_x, max_y = clipped_geometry.bounds bounds = ( math.floor(min_x / resolution) * resolution, math.floor(min_y / resolution) * resolution, math.ceil(max_x / resolution) * resolution, math.ceil(max_y / resolution) * resolution, ) width_m, height_m = bounds[2] - bounds[0], bounds[3] - bounds[1] if ( width_m > settings.spw_terrain_max_side_m or height_m > settings.spw_terrain_max_side_m ): raise AppError( code="SPW_TERRAIN_SELECTION_TOO_LARGE", message="SPW terrain selection exceeds the configured side limit", status_code=422, ) width, height = ( max(1, round(width_m / resolution)), max(1, round(height_m / resolution)), ) if width * height > settings.spw_terrain_max_pixels: raise AppError( code="SPW_TERRAIN_SELECTION_TOO_LARGE", message="SPW terrain selection exceeds the configured cell limit", details={ "pixel_count": width * height, "max_pixels": settings.spw_terrain_max_pixels, }, status_code=422, ) window = window_from_bounds(*bounds, transform=source.transform) band = source.read( 1, window=window, out_shape=(height, width), masked=True, resampling=Resampling.bilinear, ) output_transform = from_bounds(*bounds, width, height) outside_scope = geometry_mask( [mapping(clipped_geometry)], out_shape=(height, width), transform=output_transform, invert=False, ) values = np.asarray(np.ma.getdata(band), dtype="float32") invalid = ( np.ma.getmaskarray(band) | outside_scope | ~np.isfinite(values) ) if source.nodata is not None: invalid |= np.isclose( values.astype("float64"), float(source.nodata) ) values[invalid] = SpwTerrainService.NODATA valid = values[~invalid] if valid.size == 0: raise AppError( code="SPW_TERRAIN_NO_VALID_DATA", message="SPW MNT contains no valid cells in this selection", status_code=422, ) if float(valid.min()) < -100.0 or float(valid.max()) > 1000.0: raise AppError( code="SPW_TERRAIN_SOURCE_INVALID_VALUES", message="SPW MNT contains implausible elevations for Wallonia", details={ "minimum": float(valid.min()), "maximum": float(valid.max()), }, status_code=409, ) profile = { "driver": "GTiff", "width": width, "height": height, "count": 1, "dtype": "float32", "crs": SpwTerrainService.SOURCE_CRS, "transform": output_transform, "nodata": SpwTerrainService.NODATA, "compress": "deflate", "predictor": 3, } with MemoryFile() as memory: with memory.open(**profile) as output: output.write(values, 1) content = memory.read() return content, { "width": width, "height": height, "valid_pixel_count": int(valid.size), "bbox_epsg3812": list(bounds), "source_width": int(source.width), "source_height": int(source.height), "source_nodata": None if source.nodata is None else float(source.nodata), "source_resolution_m": 1.0, "analysis_resolution_m": resolution, "elevation_min_m": float(valid.min()), "elevation_max_m": float(valid.max()), } except AppError: raise except Exception as exc: raise AppError( code="SPW_TERRAIN_SOURCE_READ_FAILED", message="The provisioned SPW MNT could not be read", details={"reason": str(exc)}, status_code=500, ) from exc @staticmethod def _cached_dataset(db, project_id: UUID, filename: str) -> Dataset | None: candidate = ( db.query(Dataset) .filter( Dataset.project_id == project_id, Dataset.name == filename, Dataset.source_name == SpwTerrainService.PROVIDER, Dataset.status == "ready", ) .order_by(Dataset.imported_at.desc()) .first() ) return ( candidate if candidate and candidate.storage_path and Path(candidate.storage_path).is_file() else None ) @staticmethod def acquire( db, project_id: UUID, payload: SpwTerrainAcquireRequest, *, settings: Settings | None = None, ) -> dict[str, Any]: resolved = settings or get_settings() if not resolved.spw_terrain_enabled: raise AppError( code="SPW_TERRAIN_NOT_CONFIGURED", message="SPW terrain bounded analysis is disabled", status_code=503, ) source_path = SpwTerrainService._source_path(resolved) source_sha256 = SpwTerrainService._source_sha256(resolved) if not source_path.is_file() or source_sha256 is None: raise AppError( code="SPW_TERRAIN_SOURCE_NOT_PROVISIONED", message="The official SPW MNT source and checksum have not been provisioned on this runtime", details={ "expected_path": str(source_path), "expected_checksum_path": str( source_path.with_name(SpwTerrainService.SOURCE_SHA256_FILENAME) ), "operator_command": "python scripts/provision_spw_terrain_source.py", }, status_code=503, ) scope, bbox_4326 = SpwTerrainService._scope_geometry(db, project_id, payload) resolution = float( payload.resolution_m or resolved.spw_terrain_analysis_resolution_m ) identity = { "product_key": SpwTerrainService.PRODUCT_KEY, "bbox_epsg4326": [round(float(value), 8) for value in bbox_4326], "area_id": str(payload.area_id) if payload.area_id else None, "analysis_resolution_m": resolution, } request_hash = hashlib.sha256( json.dumps(identity, sort_keys=True).encode() ).hexdigest() filename = f"spw_mnt_2021_2022_{request_hash[:12]}_3812.tif" if not payload.force_refresh: cached = SpwTerrainService._cached_dataset(db, project_id, filename) if cached is not None: metadata = cached.source_metadata or {} return SpwTerrainAcquisitionResult( output_dataset_id=cached.id, reused=True, provider=SpwTerrainService.PROVIDER, product_key=SpwTerrainService.PRODUCT_KEY, display_name=SpwTerrainService.DISPLAY_NAME, surface_model=SpwTerrainService.SURFACE_MODEL, native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M, resolution_m=resolution, width=int((cached.metadata_json or {}).get("width", 0)), height=int((cached.metadata_json or {}).get("height", 0)), valid_pixel_count=int(metadata.get("valid_pixel_count", 0)), nodata_value=SpwTerrainService.NODATA, bbox_epsg4326=bbox_4326, bbox_epsg3812=list(metadata.get("bbox_epsg3812") or []), vertical_reference=SpwTerrainService.VERTICAL_REFERENCE, acquisition_period=SpwTerrainService.ACQUISITION_PERIOD, attribution=SpwTerrainService.ATTRIBUTION, limitation_message=SpwTerrainService.LIMITATION, ).model_dump(mode="json") content, validation = SpwTerrainService._read_source_window( source_path, scope, resolution, resolved ) acquired_at = datetime.now(UTC) dataset = DatasetService.import_raster_bytes( db, project_id=project_id, area_id=payload.area_id, filename=filename, content=content, source="SPW Relief de la Wallonie MNT 2021-2022 operator-provisioned GeoTIFF", source_name=SpwTerrainService.PROVIDER, temporal_series_key=f"spw:terrain:mnt:{request_hash[:24]}", observed_at=datetime(2022, 3, 5, 23, 59, 59, tzinfo=UTC), valid_from=datetime(2021, 2, 19, tzinfo=UTC), valid_to=datetime(2022, 3, 5, 23, 59, 59, tzinfo=UTC), temporal_granularity="period", source_version="RELIEF_WALLONIE_MNT_1M_2021_2022", source_metadata={ "provider": SpwTerrainService.PROVIDER, "product_key": SpwTerrainService.PRODUCT_KEY, "product_display_name": SpwTerrainService.DISPLAY_NAME, "surface_model": SpwTerrainService.SURFACE_MODEL, "source_crs": SpwTerrainService.SOURCE_CRS, "source_resolution_m": SpwTerrainService.SOURCE_RESOLUTION_M, "analysis_resolution_m": validation["analysis_resolution_m"], "valid_pixel_count": validation["valid_pixel_count"], "bbox_epsg4326": bbox_4326, "bbox_epsg3812": validation["bbox_epsg3812"], "coverage_zones": ["wallonia"], "vertical_reference": SpwTerrainService.VERTICAL_REFERENCE, "vertical_unit": "m", "vertical_unit_label": SpwTerrainService.VERTICAL_UNIT_LABEL, "acquisition_period": SpwTerrainService.ACQUISITION_PERIOD, "catalog_url": SpwTerrainService.CATALOG_URL, "download_url": SpwTerrainService.DOWNLOAD_URL, "attribution": SpwTerrainService.ATTRIBUTION, "license_note": SpwTerrainService.LICENSE_NOTE, "limitation_message": SpwTerrainService.LIMITATION, }, provenance_metadata={ "acquisition": "operator_provisioned_official_archive_bounded_window", "acquired_at": acquired_at.isoformat(), "request_hash": request_hash, "source_filename": SpwTerrainService.SOURCE_FILENAME, "source_sha256": source_sha256, "derived_sha256": hashlib.sha256(content).hexdigest(), "resampling": "bilinear", "validation": validation, }, ) return SpwTerrainAcquisitionResult( output_dataset_id=dataset.id, reused=False, provider=SpwTerrainService.PROVIDER, product_key=SpwTerrainService.PRODUCT_KEY, display_name=SpwTerrainService.DISPLAY_NAME, surface_model=SpwTerrainService.SURFACE_MODEL, native_resolution_m=SpwTerrainService.SOURCE_RESOLUTION_M, resolution_m=validation["analysis_resolution_m"], width=validation["width"], height=validation["height"], valid_pixel_count=validation["valid_pixel_count"], nodata_value=SpwTerrainService.NODATA, bbox_epsg4326=bbox_4326, bbox_epsg3812=validation["bbox_epsg3812"], vertical_reference=SpwTerrainService.VERTICAL_REFERENCE, acquisition_period=SpwTerrainService.ACQUISITION_PERIOD, attribution=SpwTerrainService.ATTRIBUTION, limitation_message=SpwTerrainService.LIMITATION, ).model_dump(mode="json")