diff --git a/backend/app/services/dhmv_acquisition_service.py b/backend/app/services/dhmv_acquisition_service.py index b0c4128e..9afe9550 100644 --- a/backend/app/services/dhmv_acquisition_service.py +++ b/backend/app/services/dhmv_acquisition_service.py @@ -46,6 +46,7 @@ class DhmvAcquisitionService: NODATA = -9999.0 ATTRIBUTION = "Bron: Digitaal Vlaanderen, Digitaal Hoogtemodel Vlaanderen II" LICENSE_NOTE = "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen." + WCS_TILE_SIDE_M = 10_000.0 DTM_CATALOG_URL = ( "https://www.vlaanderen.be/datavindplaats/catalogus/" "digitaal-hoogtemodel-vlaanderen-ii-dtm-raster-1-m" @@ -204,6 +205,40 @@ class DhmvAcquisitionService: "height": height, } + @staticmethod + def _wcs_request_url( + settings: Settings, + product: DhmvProduct, + bounds: tuple[float, float, float, float], + resolution_m: float, + ) -> str: + query = [ + ("SERVICE", "WCS"), + ("VERSION", "2.0.1"), + ("REQUEST", "GetCoverage"), + ("COVERAGEID", product.coverage_id), + ("FORMAT", "image/tiff"), + ("SUBSET", f"x({bounds[0]:.3f},{bounds[2]:.3f})"), + ("SUBSET", f"y({bounds[1]:.3f},{bounds[3]:.3f})"), + ("SCALEFACTOR", f"{resolution_m / product.native_resolution_m:g}"), + ] + return f"{settings.dhmv_wcs_url}?{urlencode(query)}" + + @staticmethod + def _tile_bounds(prepared: dict[str, Any]) -> list[tuple[float, float, float, float]]: + min_x, min_y, max_x, max_y = prepared["bbox_epsg31370"] + tiles: list[tuple[float, float, float, float]] = [] + y = min_y + while y < max_y: + tile_max_y = min(y + DhmvAcquisitionService.WCS_TILE_SIDE_M, max_y) + x = min_x + while x < max_x: + tile_max_x = min(x + DhmvAcquisitionService.WCS_TILE_SIDE_M, max_x) + tiles.append((x, y, tile_max_x, tile_max_y)) + x = tile_max_x + y = tile_max_y + return tiles + @staticmethod def _scope_geometry(db, project_id: UUID, area_id: UUID | None, bbox_epsg4326: list[float]): if not db.get(Project, project_id): @@ -291,6 +326,113 @@ class DhmvAcquisitionService: status_code=502, ) + @staticmethod + def _mosaic_geotiffs(coverages: list[bytes]) -> 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: + raise AppError( + code="RASTER_PROCESSING_UNAVAILABLE", + message="Rasterio is required to assemble tiled DHMV coverages", + status_code=503, + ) from exc + + memories = [MemoryFile(content) for content in coverages] + sources = [] + try: + sources = [memory.open() for memory in memories] + crs_values = {str(source.crs) for source in sources} + band_counts = {source.count for source in sources} + resolutions = { + (round(abs(float(source.res[0])), 6), round(abs(float(source.res[1])), 6)) + for source in sources + } + if crs_values != {DhmvAcquisitionService.SOURCE_CRS} or band_counts != {1} or len(resolutions) != 1: + raise AppError( + code="DHMV_TILE_MISMATCH", + message="DHMV coverage tiles do not share one CRS, band layout and resolution", + status_code=502, + ) + mosaic, transform = merge( + sources, + nodata=DhmvAcquisitionService.NODATA, + dtype="float32", + ) + profile = sources[0].profile.copy() + profile.pop("blockxsize", None) + profile.pop("blockysize", None) + profile.update( + driver="GTiff", + width=int(mosaic.shape[2]), + height=int(mosaic.shape[1]), + count=1, + dtype="float32", + crs=DhmvAcquisitionService.SOURCE_CRS, + transform=transform, + nodata=DhmvAcquisitionService.NODATA, + compress="deflate", + predictor=3, + ) + with MemoryFile() as output_memory: + with output_memory.open(**profile) as output: + output.write(mosaic) + return output_memory.read() + except AppError: + raise + except Exception as exc: + raise AppError( + code="DHMV_TILE_MOSAIC_FAILED", + message="DHMV coverage tiles could not be assembled into one georeferenced raster", + details={"reason": str(exc)}, + status_code=502, + ) from exc + finally: + for source in sources: + source.close() + for memory in memories: + memory.close() + + @staticmethod + def _fetch_coverage( + prepared: dict[str, Any], + settings: Settings, + opener: Callable[..., Any] | None = None, + ) -> tuple[bytes, dict[str, Any]]: + product: DhmvProduct = prepared["product"] + request_urls = [ + DhmvAcquisitionService._wcs_request_url( + settings, + product, + bounds, + prepared["resolution_m"], + ) + for bounds in DhmvAcquisitionService._tile_bounds(prepared) + ] + raw_hash = hashlib.sha256() + coverage_hash = hashlib.sha256() + content_types: list[str] = [] + coverages: list[bytes] = [] + for request_url in request_urls: + raw_content, content_type = DhmvAcquisitionService._fetch(request_url, settings, opener) + coverage_content = DhmvAcquisitionService._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_content).to_bytes(8, "big")) + coverage_hash.update(coverage_content) + content_types.append(content_type) + coverages.append(coverage_content) + return DhmvAcquisitionService._mosaic_geotiffs(coverages), { + "tile_count": len(request_urls), + "request_urls": request_urls, + "response_content_types": content_types, + "response_sha256": raw_hash.hexdigest(), + "coverage_sha256": coverage_hash.hexdigest(), + } + @staticmethod def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]: try: @@ -413,8 +555,7 @@ class DhmvAcquisitionService: limitation_message=product.limitation_message, ).model_dump(mode="json") - raw_content, content_type = DhmvAcquisitionService._fetch(prepared["request_url"], resolved_settings, opener) - coverage_content = DhmvAcquisitionService._extract_geotiff(raw_content, content_type) + coverage_content, transfer = DhmvAcquisitionService._fetch_coverage(prepared, resolved_settings, opener) normalized_content, validation = DhmvAcquisitionService._normalize_raster(coverage_content, scope_geometry, prepared) acquired_at = datetime.now(UTC) dataset = DatasetService.import_raster_bytes( @@ -460,13 +601,15 @@ class DhmvAcquisitionService: "coverage_scope": "municipality" if payload.area_id else "bounded_selection", }, provenance_metadata={ - "acquisition": "explicit_bounded_wcs_coverage", + "acquisition": "explicit_bounded_tiled_wcs_coverage", "acquired_at": acquired_at.isoformat(), "request_hash": prepared["request_hash"], "request_url": prepared["request_url"], - "response_content_type": content_type, - "response_sha256": hashlib.sha256(raw_content).hexdigest(), - "coverage_sha256": hashlib.sha256(coverage_content).hexdigest(), + "tile_count": transfer["tile_count"], + "tile_request_urls": transfer["request_urls"], + "response_content_types": transfer["response_content_types"], + "response_sha256": transfer["response_sha256"], + "coverage_sha256": transfer["coverage_sha256"], "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 2ea026e8..4200c985 100644 --- a/backend/tests/test_sprint205_dhmv_terrain.py +++ b/backend/tests/test_sprint205_dhmv_terrain.py @@ -167,6 +167,29 @@ def test_dhmv_request_rejects_unsafe_size_and_resolution() -> None: DhmvAcquireRequest.model_validate(payload.model_dump()) +def test_dhmv_large_scope_is_bounded_into_mosaicable_wcs_tiles() -> None: + prepared = DhmvAcquisitionService._prepared_request( + lambert_bbox_payload(side_m=15_000.0), + Settings(_env_file=None), + ) + + tile_bounds = DhmvAcquisitionService._tile_bounds(prepared) + + assert len(tile_bounds) == 4 + assert all(bounds[2] - bounds[0] <= 10_000.0 for bounds in tile_bounds) + assert all(bounds[3] - bounds[1] <= 10_000.0 for bounds in tile_bounds) + + left = elevation_tiff(left=200_000, top=210_100, width=20, height=20) + right = elevation_tiff(left=200_100, top=210_100, width=20, height=20) + mosaic = DhmvAcquisitionService._mosaic_geotiffs([left, right]) + with MemoryFile(mosaic) as memory, memory.open() as dataset: + assert dataset.crs.to_epsg() == 31370 + assert dataset.res == pytest.approx((5.0, 5.0)) + assert dataset.width == 40 + assert dataset.height == 20 + assert dataset.nodata == -9999.0 + + 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)