Handle DHMV WCS edge grid rounding

backend/app/services/dhmv_acquisition_service.py

@@ -51,6 +51,8 @@ class DhmvAcquisitionService:
     WCS_REQUEST_INTERVAL_SECONDS = 2.0
     WCS_RETRY_DELAY_SECONDS = 4.0
     WCS_TRANSIENT_STATUS_CODES = frozenset({400, 429, 502, 503, 504})
+    WCS_EDGE_RESOLUTION_REL_TOLERANCE = 0.05
+    WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M = 0.25
     DTM_CATALOG_URL = (
         "https://www.vlaanderen.be/datavindplaats/catalogus/"
         "digitaal-hoogtemodel-vlaanderen-ii-dtm-raster-1-m"
@@ -349,11 +351,14 @@ class DhmvAcquisitionService:
         )
 
     @staticmethod
-    def _mosaic_geotiffs(coverages: list[bytes], expected_resolution_m: float | None = None) -> bytes:
+    def _mosaic_geotiffs(
+        coverages: list[bytes],
+        expected_resolution_m: float | None = None,
+        diagnostics: dict[str, Any] | None = None,
+    ) -> bytes:
         if len(coverages) == 1:
             return coverages[0]
         try:
-            import rasterio
             from rasterio.io import MemoryFile
             from rasterio.merge import merge
         except ImportError as exc:
@@ -370,14 +375,23 @@ class DhmvAcquisitionService:
             target_resolution = float(expected_resolution_m or abs(float(sources[0].res[0])))
             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": [abs(float(source.res[0])), abs(float(source.res[1]))],
+                    "resolution": tile_resolutions[index],
                 }
                 for index, source in enumerate(sources)
                 if not all(
-                    math.isclose(abs(float(value)), target_resolution, rel_tol=0.02, abs_tol=0.05)
+                    math.isclose(
+                        abs(float(value)),
+                        target_resolution,
+                        rel_tol=DhmvAcquisitionService.WCS_EDGE_RESOLUTION_REL_TOLERANCE,
+                        abs_tol=DhmvAcquisitionService.WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M,
+                    )
                     for value in source.res
                 )
             ]
@@ -393,6 +407,23 @@ class DhmvAcquisitionService:
                     },
                     status_code=502,
                 )
+            harmonized_tile_indexes = [
+                index
+                for index, resolution in enumerate(tile_resolutions)
+                if not all(
+                    math.isclose(value, target_resolution, 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": target_resolution,
+                        "harmonized_tile_indexes": harmonized_tile_indexes,
+                        "harmonization_method": "rasterio_merge_target_resolution" if harmonized_tile_indexes else None,
+                    }
+                )
             mosaic, transform = merge(
                 sources,
                 res=(target_resolution, target_resolution),
@@ -471,19 +502,25 @@ class DhmvAcquisitionService:
             coverage_hash.update(coverage_content)
             content_types.append(content_type)
             coverages.append(coverage_content)
-        return DhmvAcquisitionService._mosaic_geotiffs(coverages, prepared["resolution_m"]), {
+        mosaic_diagnostics: dict[str, Any] = {}
+        mosaic = DhmvAcquisitionService._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
     def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]:
         try:
             import numpy as np
-            import rasterio
             from rasterio.io import MemoryFile
             from rasterio.mask import mask
         except ImportError as exc:
@@ -656,6 +693,7 @@ class DhmvAcquisitionService:
                 "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_content).hexdigest(),
                 "bbox_epsg4326": prepared["bbox_epsg4326"],
                 "bbox_epsg31370": prepared["bbox_epsg31370"],

backend/tests/test_sprint205_dhmv_terrain.py

@@ -115,6 +115,24 @@ def elevation_tiff(*, left: float, top: float, width: int, height: int, resoluti
         return memory.read()
 
 
+def edge_elevation_tiff(*, left: float, top: float, x_resolution: float, y_resolution: float = 5.0) -> bytes:
+    rows, columns = np.indices((20, 20))
+    values = (20.0 + columns * 0.5 + rows).astype("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=-9999.0,
+        ) as output:
+            output.write(values, 1)
+        return memory.read()
+
+
 def multipart_tiff(content: bytes) -> tuple[bytes, str]:
     boundary = "wcs-test"
     payload = (
@@ -190,6 +208,28 @@ def test_dhmv_large_scope_is_bounded_into_mosaicable_wcs_tiles() -> None:
         assert dataset.nodata == -9999.0
 
 
+def test_dhmv_mosaic_harmonizes_only_bounded_wcs_edge_grid_rounding() -> None:
+    regular = elevation_tiff(left=200_000, top=210_100, width=20, height=20)
+    rounded_edge = edge_elevation_tiff(left=200_100, top=210_100, x_resolution=4.76555, y_resolution=5.0008)
+    diagnostics: dict[str, object] = {}
+
+    mosaic = DhmvAcquisitionService._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_elevation_tiff(left=200_100, top=210_100, x_resolution=4.5)
+    with pytest.raises(AppError) as exc_info:
+        DhmvAcquisitionService._mosaic_geotiffs([regular, unsafe_edge], expected_resolution_m=5.0)
+    assert exc_info.value.code == "DHMV_TILE_MISMATCH"
+
+
 def test_dhmv_multipart_geotiff_is_extracted_and_invalid_response_fails_closed() -> None:
     tiff = elevation_tiff(left=200_000, top=210_100, width=20, height=20)
     multipart, content_type = multipart_tiff(tiff)