Tile large DHMV coverage requests

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