Handle VMM WCS edge grid rounding

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

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"""<?xml version="1.0"?>
     <ExceptionReport xmlns="http://www.opengis.net/ows/1.1">