diff --git a/backend/app/services/flood_hazard_acquisition_service.py b/backend/app/services/flood_hazard_acquisition_service.py index 13846b5e..1da8b1df 100644 --- a/backend/app/services/flood_hazard_acquisition_service.py +++ b/backend/app/services/flood_hazard_acquisition_service.py @@ -57,6 +57,10 @@ class FloodHazardAcquisitionService: WCS_REQUEST_INTERVAL_SECONDS = 1.0 WCS_RETRY_DELAY_SECONDS = 3.0 WCS_TRANSIENT_STATUS_CODES = frozenset({400, 429, 502, 503, 504}) + # The VMM WCS rounds the grid size of partial edge tiles to an integer + # number of cells. Keep that provider artefact bounded and auditable. + WCS_EDGE_RESOLUTION_REL_TOLERANCE = 0.05 + WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M = 0.25 SERVICE_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/publieke-inspire-coverage-service-van-ogrk" LIMITATION = ( "Gemodelleerde maximale overstromingsdiepte voor een vast kans- en klimaatscenario. " @@ -319,7 +323,11 @@ class FloodHazardAcquisitionService: ) @staticmethod - def _mosaic_geotiffs(coverages: list[bytes], expected_resolution_m: float) -> bytes: + def _mosaic_geotiffs( + coverages: list[bytes], + expected_resolution_m: float, + diagnostics: dict[str, Any] | None = None, + ) -> bytes: if len(coverages) == 1: return coverages[0] try: @@ -331,11 +339,63 @@ class FloodHazardAcquisitionService: sources = [] try: sources = [memory.open() for memory in memories] - for source in sources: - if source.crs is None or source.crs.to_epsg() != 31370 or source.count != 1: - raise AppError(code="FLOOD_HAZARD_TILE_MISMATCH", message="VMM coverage tiles do not share the governed CRS and band layout", status_code=502) - if not all(math.isclose(abs(float(value)), expected_resolution_m, rel_tol=0.02, abs_tol=0.05) for value in source.res): - raise AppError(code="FLOOD_HAZARD_TILE_MISMATCH", message="VMM coverage tile resolution differs from the governed request", status_code=502) + invalid_crs = [ + index + for index, source in enumerate(sources) + if source.crs is None or source.crs.to_epsg() != 31370 + ] + invalid_bands = [index for index, source in enumerate(sources) if source.count != 1] + tile_resolutions = [ + [abs(float(source.res[0])), abs(float(source.res[1]))] + for source in sources + ] + invalid_resolution = [ + { + "tile_index": index, + "resolution": tile_resolutions[index], + } + for index, source in enumerate(sources) + if not all( + math.isclose( + abs(float(value)), + expected_resolution_m, + rel_tol=FloodHazardAcquisitionService.WCS_EDGE_RESOLUTION_REL_TOLERANCE, + abs_tol=FloodHazardAcquisitionService.WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M, + ) + for value in source.res + ) + ] + if invalid_crs or invalid_bands or invalid_resolution: + raise AppError( + code="FLOOD_HAZARD_TILE_MISMATCH", + message="VMM coverage tiles do not match the governed CRS, band layout and resolution", + details={ + "invalid_crs_tile_indexes": invalid_crs, + "invalid_band_tile_indexes": invalid_bands, + "invalid_resolution_tiles": invalid_resolution, + "expected_resolution_m": expected_resolution_m, + }, + status_code=502, + ) + harmonized_tile_indexes = [ + index + for index, resolution in enumerate(tile_resolutions) + if not all( + math.isclose(value, expected_resolution_m, rel_tol=0.02, abs_tol=0.05) + for value in resolution + ) + ] + if diagnostics is not None: + diagnostics.update( + { + "source_tile_resolutions_m": tile_resolutions, + "target_resolution_m": expected_resolution_m, + "harmonized_tile_indexes": harmonized_tile_indexes, + "harmonization_method": ( + "rasterio_merge_target_resolution" if harmonized_tile_indexes else None + ), + } + ) mosaic, transform = merge(sources, res=(expected_resolution_m, expected_resolution_m), nodata=0.0, dtype="float32") profile = sources[0].profile.copy() profile.pop("blockxsize", None) @@ -378,16 +438,25 @@ class FloodHazardAcquisitionService: time.sleep(FloodHazardAcquisitionService.WCS_RETRY_DELAY_SECONDS) raw_content, content_type = FloodHazardAcquisitionService._fetch(request_url, settings, opener) coverage = FloodHazardAcquisitionService._extract_geotiff(raw_content, content_type) - raw_hash.update(len(raw_content).to_bytes(8, "big")); raw_hash.update(raw_content) - coverage_hash.update(len(coverage).to_bytes(8, "big")); coverage_hash.update(coverage) + raw_hash.update(len(raw_content).to_bytes(8, "big")) + raw_hash.update(raw_content) + coverage_hash.update(len(coverage).to_bytes(8, "big")) + coverage_hash.update(coverage) content_types.append(content_type) coverages.append(coverage) - return FloodHazardAcquisitionService._mosaic_geotiffs(coverages, prepared["resolution_m"]), { + mosaic_diagnostics: dict[str, Any] = {} + mosaic = FloodHazardAcquisitionService._mosaic_geotiffs( + coverages, + prepared["resolution_m"], + diagnostics=mosaic_diagnostics, + ) + return mosaic, { "tile_count": len(request_urls), "request_urls": request_urls, "response_content_types": content_types, "response_sha256": raw_hash.hexdigest(), "coverage_sha256": coverage_hash.hexdigest(), + "grid_harmonization": mosaic_diagnostics, } @staticmethod @@ -417,7 +486,8 @@ class FloodHazardAcquisitionService: normalized_m[positive] = raw_cm[positive] / 100.0 valid_values = normalized_m[positive].astype("float64") profile = source.profile.copy() - profile.pop("blockxsize", None); profile.pop("blockysize", None) + profile.pop("blockxsize", None) + profile.pop("blockysize", None) profile.update(driver="GTiff", width=normalized_m.shape[1], height=normalized_m.shape[0], count=1, dtype="float32", crs=FloodHazardAcquisitionService.SOURCE_CRS, transform=transform, nodata=FloodHazardAcquisitionService.NODATA, compress="deflate", predictor=3) with MemoryFile() as output_memory: with output_memory.open(**profile) as output: @@ -536,6 +606,7 @@ class FloodHazardAcquisitionService: "response_content_types": transfer["response_content_types"], "response_sha256": transfer["response_sha256"], "coverage_sha256": transfer["coverage_sha256"], + "grid_harmonization": transfer["grid_harmonization"], "normalized_sha256": hashlib.sha256(normalized).hexdigest(), "bbox_epsg4326": prepared["bbox_epsg4326"], "bbox_epsg31370": prepared["bbox_epsg31370"], diff --git a/backend/tests/test_sprint208_vmm_flood_hazard.py b/backend/tests/test_sprint208_vmm_flood_hazard.py index d8543475..e9025d37 100644 --- a/backend/tests/test_sprint208_vmm_flood_hazard.py +++ b/backend/tests/test_sprint208_vmm_flood_hazard.py @@ -5,7 +5,6 @@ from uuid import uuid4 import numpy as np import pytest -import rasterio from fastapi.testclient import TestClient from pyproj import Transformer from rasterio.io import MemoryFile @@ -104,6 +103,23 @@ def depth_tiff(*, normalized_metres: bool = False) -> bytes: return memory.read() +def edge_depth_tiff(*, left: float, top: float, x_resolution: float, y_resolution: float = 5.0) -> bytes: + values = np.full((20, 20), 100.0, dtype="float32") + with MemoryFile() as memory: + with memory.open( + driver="GTiff", + width=20, + height=20, + count=1, + dtype="float32", + crs="EPSG:31370", + transform=from_origin(left, top, x_resolution, y_resolution), + nodata=0.0, + ) as output: + output.write(values, 1) + return memory.read() + + def test_flood_hazard_registry_is_complete_and_semantically_honest() -> None: products = FloodHazardAcquisitionService.list_products() @@ -151,6 +167,39 @@ def test_flood_hazard_tiles_stay_below_the_observed_vmm_coverage_limit() -> None ) +def test_flood_hazard_mosaic_harmonizes_only_bounded_wcs_edge_grid_rounding() -> None: + regular = edge_depth_tiff(left=200_000, top=210_100, x_resolution=5.0) + rounded_edge = edge_depth_tiff( + left=200_100, + top=210_100, + x_resolution=4.76555, + y_resolution=5.0008, + ) + diagnostics: dict[str, object] = {} + + mosaic = FloodHazardAcquisitionService._mosaic_geotiffs( + [regular, rounded_edge], + expected_resolution_m=5.0, + diagnostics=diagnostics, + ) + + with MemoryFile(mosaic) as memory, memory.open() as dataset: + assert dataset.res == pytest.approx((5.0, 5.0)) + assert diagnostics["harmonized_tile_indexes"] == [1] + assert diagnostics["harmonization_method"] == "rasterio_merge_target_resolution" + + unsafe_edge = edge_depth_tiff(left=200_100, top=210_100, x_resolution=4.5) + with pytest.raises(AppError) as exc_info: + FloodHazardAcquisitionService._mosaic_geotiffs( + [regular, unsafe_edge], + expected_resolution_m=5.0, + ) + assert exc_info.value.code == "FLOOD_HAZARD_TILE_MISMATCH" + assert exc_info.value.details["invalid_resolution_tiles"] == [ + {"tile_index": 1, "resolution": [4.5, 5.0]} + ] + + def test_flood_hazard_xml_provider_error_is_exposed_without_losing_the_canonical_error() -> None: response = b"""