From 512e15f24b3937ac8d8b781fa703ed34f4968281 Mon Sep 17 00:00:00 2001 From: Codex Date: Thu, 16 Jul 2026 12:47:53 +0200 Subject: [PATCH] Handle DHMV WCS edge grid rounding --- .../app/services/dhmv_acquisition_service.py | 50 ++++++++++++++++--- backend/tests/test_sprint205_dhmv_terrain.py | 40 +++++++++++++++ 2 files changed, 84 insertions(+), 6 deletions(-) diff --git a/backend/app/services/dhmv_acquisition_service.py b/backend/app/services/dhmv_acquisition_service.py index e410e019..a3f2db8f 100644 --- a/backend/app/services/dhmv_acquisition_service.py +++ b/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"], diff --git a/backend/tests/test_sprint205_dhmv_terrain.py b/backend/tests/test_sprint205_dhmv_terrain.py index 6986c35b..94434b9c 100644 --- a/backend/tests/test_sprint205_dhmv_terrain.py +++ b/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)