Complete regional raster exploration

backend/app/services/flood_hazard_analysis_service.py

@@ -16,11 +16,13 @@ from app.core.errors import AppError
 from app.models import Area, Dataset
 from app.schemas.flood_hazard import (
     FloodHazardMetric,
+    FloodHazardPartitionSelectionRequest,
     FloodHazardSelectionRequest,
     FloodHazardSelectionResponse,
     FloodHazardSelectionSummary,
 )
 from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService
+from app.services.raster_partition_analysis_service import RasterPartitionAnalysisService
 
 
 class FloodHazardAnalysisService:
@@ -170,6 +172,8 @@ class FloodHazardAnalysisService:
         primary = metrics[0]
         response = FloodHazardSelectionResponse(
             dataset_id=dataset.id,
+            dataset_ids=[dataset.id],
+            partition_count=1,
             product_key=product.key,
             mechanism=product.mechanism,
             climate_context=product.climate_context,
@@ -195,6 +199,129 @@ class FloodHazardAnalysisService:
         )
         return response.model_dump(mode="json")
 
+    @staticmethod
+    def analyze_partitions(
+        db,
+        project_id: UUID,
+        payload: FloodHazardPartitionSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        product = FloodHazardAcquisitionService._products().get(payload.product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="FLOOD_HAZARD_PRODUCT_NOT_SUPPORTED",
+                message="Select a governed VMM fluvial or pluvial flood-depth scenario",
+                details={"product_key": payload.product_key},
+                status_code=422,
+            )
+        selection_4326 = FloodHazardAnalysisService._selection_geometry(db, project_id, payload)
+        partition = RasterPartitionAnalysisService.select(
+            db,
+            project_id,
+            source_name=FloodHazardAcquisitionService.PROVIDER,
+            product_key=product.key,
+            selection_geometry_4326=selection_4326,
+            nodata=FloodHazardAcquisitionService.NODATA,
+            max_pixels=resolved_settings.flood_hazard_max_pixels,
+        )
+        try:
+            import numpy as np
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Numpy is required for partitioned flood-hazard analysis",
+                status_code=503,
+            ) from exc
+
+        raw = partition.values
+        valid = (
+            partition.selected_cells
+            & np.isfinite(raw)
+            & (raw != FloodHazardAcquisitionService.NODATA)
+            & (raw > 0.0)
+        )
+        values = raw[valid]
+        selected_cell_count = int(partition.selected_cells.sum())
+        inundated_cell_count = int(values.size)
+        cell_area_m2 = partition.resolution_x * partition.resolution_y
+
+        def metric(key: str, label: str, value: float, unit: str, method: str) -> FloodHazardMetric:
+            return FloodHazardMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        inundated_area_ha = inundated_cell_count * cell_area_m2 / 10_000.0
+        metrics = [
+            metric(
+                "modelled_inundated_area_ha",
+                "Gemodelleerd overstroomd oppervlak",
+                inundated_area_ha,
+                "ha",
+                "positive_depth_cells_times_cell_area",
+            ),
+            metric(
+                "modelled_inundated_share_pct",
+                "Aandeel selectie met gemodelleerde diepte",
+                inundated_cell_count / max(1, selected_cell_count) * 100.0,
+                "%",
+                "positive_depth_cells_divided_by_selected_cells",
+            ),
+        ]
+        if inundated_cell_count:
+            metrics.extend(
+                [
+                    metric("modelled_depth_mean_m", "Gemiddelde gemodelleerde maximumdiepte", values.mean(), "m", "mean_positive_depth_cells"),
+                    metric("modelled_depth_p90_m", "90e percentiel gemodelleerde maximumdiepte", np.percentile(values, 90), "m", "percentile_90_positive_depth_cells"),
+                    metric("modelled_depth_max_m", "Hoogste gemodelleerde maximumdiepte", values.max(), "m", "maximum_positive_depth_cells"),
+                    metric(
+                        "modelled_max_depth_area_integral_m3",
+                        "Diepte-oppervlakte-integraal (geen gelijktijdig volume)",
+                        values.sum() * cell_area_m2,
+                        "m3",
+                        "sum_local_max_depth_times_cell_area",
+                    ),
+                ]
+            )
+        primary = metrics[0]
+        first_dataset = partition.datasets[0]
+        response = FloodHazardSelectionResponse(
+            dataset_id=first_dataset.id,
+            dataset_ids=[dataset.id for dataset in partition.datasets],
+            partition_count=len(partition.datasets),
+            product_key=product.key,
+            mechanism=product.mechanism,
+            climate_context=product.climate_context,
+            probability_class=product.probability_class,
+            return_period_years=product.return_period_years,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            selected_cell_count=selected_cell_count,
+            inundated_cell_count=inundated_cell_count,
+            inundated_fraction=round(inundated_cell_count / max(1, selected_cell_count), 6),
+            resolution_m=round(max(partition.resolution_x, partition.resolution_y), 4),
+            summary=FloodHazardSelectionSummary(
+                metric_label=primary.metric_label,
+                metric_value=primary.metric_value,
+                metric_unit=primary.metric_unit,
+                aggregation_method=primary.aggregation_method,
+                primary_metric_key=primary.metric_key,
+                metrics=metrics,
+            ),
+            unsupported_metrics=FloodHazardAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=(
+                f"{FloodHazardAnalysisService.LIMITATION} De selectie werd exact berekend over "
+                f"{len(partition.datasets)} persistente gemeentelijke rasterpartities."
+            ),
+            generated_at=datetime.now(UTC).isoformat(),
+        )
+        return response.model_dump(mode="json")
+
     @staticmethod
     def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
         dataset = FloodHazardAnalysisService._load_dataset(db, project_id, dataset_id)

backend/app/services/raster_partition_analysis_service.py

@@ -0,0 +1,200 @@
+from __future__ import annotations
+
+import math
+from contextlib import ExitStack
+from dataclasses import dataclass
+from pathlib import Path
+from typing import Any
+from uuid import UUID
+
+from pyproj import Transformer
+from shapely.geometry import mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.errors import AppError
+from app.models import Dataset
+
+
+@dataclass(frozen=True)
+class RasterPartitionSelection:
+    datasets: list[Dataset]
+    values: Any
+    selected_cells: Any
+    resolution_x: float
+    resolution_y: float
+
+
+class RasterPartitionAnalysisService:
+    MAX_PARTITIONS = 64
+
+    @staticmethod
+    def _bbox_intersects(dataset: Dataset, bbox: tuple[float, float, float, float]) -> bool:
+        source_bbox = (dataset.source_metadata or {}).get("bbox_epsg4326")
+        if not isinstance(source_bbox, list) or len(source_bbox) != 4:
+            return True
+        try:
+            min_x, min_y, max_x, max_y = (float(value) for value in source_bbox)
+        except (TypeError, ValueError):
+            return True
+        return not (
+            max_x <= bbox[0]
+            or min_x >= bbox[2]
+            or max_y <= bbox[1]
+            or min_y >= bbox[3]
+        )
+
+    @staticmethod
+    def _candidate_datasets(
+        db,
+        project_id: UUID,
+        *,
+        source_name: str,
+        product_key: str,
+        bbox: tuple[float, float, float, float],
+    ) -> list[Dataset]:
+        rows = (
+            db.query(Dataset)
+            .filter(
+                Dataset.project_id == project_id,
+                Dataset.source_name == source_name,
+                Dataset.dataset_type == "raster",
+                Dataset.status == "ready",
+            )
+            .all()
+        )
+        candidates = [
+            dataset
+            for dataset in rows
+            if str((dataset.source_metadata or {}).get("product_key") or "") == product_key
+            and dataset.storage_path
+            and Path(dataset.storage_path).is_file()
+            and RasterPartitionAnalysisService._bbox_intersects(dataset, bbox)
+        ]
+        candidates.sort(key=lambda dataset: (str(dataset.area_id or ""), str(dataset.id)))
+        if not candidates:
+            raise AppError(
+                code="RASTER_PARTITIONS_NOT_FOUND",
+                message="No persisted raster partitions cover this selection",
+                details={"source_name": source_name, "product_key": product_key},
+                status_code=404,
+            )
+        if len(candidates) > RasterPartitionAnalysisService.MAX_PARTITIONS:
+            raise AppError(
+                code="RASTER_PARTITION_LIMIT_EXCEEDED",
+                message="The selection intersects too many raster partitions",
+                details={
+                    "partition_count": len(candidates),
+                    "max_partitions": RasterPartitionAnalysisService.MAX_PARTITIONS,
+                },
+                status_code=422,
+            )
+        return candidates
+
+    @staticmethod
+    def select(
+        db,
+        project_id: UUID,
+        *,
+        source_name: str,
+        product_key: str,
+        selection_geometry_4326,
+        nodata: float,
+        max_pixels: int,
+    ) -> RasterPartitionSelection:
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.features import geometry_mask
+            from rasterio.merge import merge
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio and numpy are required for partitioned raster analysis",
+                status_code=503,
+            ) from exc
+
+        bbox = tuple(float(value) for value in selection_geometry_4326.bounds)
+        datasets = RasterPartitionAnalysisService._candidate_datasets(
+            db,
+            project_id,
+            source_name=source_name,
+            product_key=product_key,
+            bbox=bbox,
+        )
+        transformer = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True)
+        selection_metric = shapely_transform(transformer.transform, selection_geometry_4326)
+        min_x, min_y, max_x, max_y = selection_metric.bounds
+
+        try:
+            with ExitStack() as stack:
+                sources = [stack.enter_context(rasterio.open(dataset.storage_path)) for dataset in datasets]
+                invalid_sources = [
+                    index
+                    for index, source in enumerate(sources)
+                    if source.crs is None or source.crs.to_epsg() != 31370 or source.count != 1
+                ]
+                if invalid_sources:
+                    raise AppError(
+                        code="RASTER_PARTITION_MISMATCH",
+                        message="Raster partitions do not share the governed CRS and band layout",
+                        details={"invalid_partition_indexes": invalid_sources},
+                        status_code=409,
+                    )
+                target_resolution = max(abs(float(sources[0].res[0])), abs(float(sources[0].res[1])))
+                invalid_resolutions = [
+                    {
+                        "partition_index": index,
+                        "resolution": [abs(float(source.res[0])), abs(float(source.res[1]))],
+                    }
+                    for index, source in enumerate(sources)
+                    if not all(
+                        math.isclose(abs(float(value)), target_resolution, rel_tol=0.001, abs_tol=0.01)
+                        for value in source.res
+                    )
+                ]
+                if invalid_resolutions:
+                    raise AppError(
+                        code="RASTER_PARTITION_MISMATCH",
+                        message="Raster partitions do not share one analysis resolution",
+                        details={"invalid_resolutions": invalid_resolutions},
+                        status_code=409,
+                    )
+                width = max(1, math.ceil((max_x - min_x) / target_resolution))
+                height = max(1, math.ceil((max_y - min_y) / target_resolution))
+                if width * height > max_pixels:
+                    raise AppError(
+                        code="RASTER_PARTITION_SELECTION_TOO_LARGE",
+                        message="Select a smaller rectangle for regional raster analysis",
+                        details={"pixel_count": width * height, "max_pixels": max_pixels},
+                        status_code=422,
+                    )
+                mosaic, transform = merge(
+                    sources,
+                    bounds=(min_x, min_y, max_x, max_y),
+                    res=(target_resolution, target_resolution),
+                    nodata=nodata,
+                    dtype="float32",
+                )
+                values = np.asarray(mosaic[0], dtype="float64")
+                selected_cells = geometry_mask(
+                    [mapping(selection_metric)],
+                    out_shape=values.shape,
+                    transform=transform,
+                    invert=True,
+                )
+                return RasterPartitionSelection(
+                    datasets=datasets,
+                    values=values,
+                    selected_cells=selected_cells,
+                    resolution_x=target_resolution,
+                    resolution_y=target_resolution,
+                )
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="RASTER_PARTITION_ANALYSIS_FAILED",
+                message="Persisted raster partitions could not be assembled for this selection",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc

backend/app/services/terrain_analysis_service.py

@@ -14,8 +14,15 @@ from shapely.ops import transform as shapely_transform
 from app.core.config import Settings, get_settings
 from app.core.errors import AppError
 from app.models import Area, Dataset
-from app.schemas.dhmv import TerrainMetric, TerrainSelectionRequest, TerrainSelectionResponse, TerrainSelectionSummary
+from app.schemas.dhmv import (
+    TerrainMetric,
+    TerrainPartitionSelectionRequest,
+    TerrainSelectionRequest,
+    TerrainSelectionResponse,
+    TerrainSelectionSummary,
+)
 from app.services.dhmv_acquisition_service import DhmvAcquisitionService
+from app.services.raster_partition_analysis_service import RasterPartitionAnalysisService
 
 
 class TerrainAnalysisService:
@@ -177,6 +184,8 @@ class TerrainAnalysisService:
         selected_cell_count = int(selected_cells.sum())
         response = TerrainSelectionResponse(
             dataset_id=dataset.id,
+            dataset_ids=[dataset.id],
+            partition_count=1,
             product_key=product_key,
             surface_model=surface_model,
             selection_bbox=payload.bbox,
@@ -200,6 +209,139 @@ class TerrainAnalysisService:
         )
         return response.model_dump(mode="json")
 
+    @staticmethod
+    def analyze_partitions(
+        db,
+        project_id: UUID,
+        payload: TerrainPartitionSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        product = DhmvAcquisitionService._products().get(payload.product_key.strip().lower())
+        if product is None:
+            raise AppError(
+                code="DHMV_PRODUCT_NOT_SUPPORTED",
+                message="Select a governed DHMV terrain or surface product",
+                details={"product_key": payload.product_key},
+                status_code=422,
+            )
+        selection_4326 = TerrainAnalysisService._selection_geometry(db, project_id, payload)
+        partition = RasterPartitionAnalysisService.select(
+            db,
+            project_id,
+            source_name=DhmvAcquisitionService.PROVIDER,
+            product_key=product.key,
+            selection_geometry_4326=selection_4326,
+            nodata=DhmvAcquisitionService.NODATA,
+            max_pixels=resolved_settings.dhmv_max_pixels,
+        )
+        surface_models = {
+            str((dataset.source_metadata or {}).get("surface_model") or "")
+            for dataset in partition.datasets
+        }
+        if surface_models != {product.surface_model}:
+            raise AppError(
+                code="INVALID_TERRAIN_METADATA",
+                message="DHMV partition provenance is incomplete",
+                details={"surface_models": sorted(surface_models)},
+                status_code=409,
+            )
+
+        try:
+            import numpy as np
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Numpy is required for partitioned terrain analysis",
+                status_code=503,
+            ) from exc
+
+        raw = partition.values
+        invalid = ~np.isfinite(raw) | (raw == DhmvAcquisitionService.NODATA)
+        valid_mask = partition.selected_cells & ~invalid
+        values = raw[valid_mask]
+        if values.size == 0:
+            raise AppError(
+                code="TERRAIN_NO_VALID_DATA",
+                message="No valid DHMV height cells occur in this selection",
+                status_code=422,
+            )
+        slope_values = np.asarray([], dtype="float64")
+        if raw.shape[0] >= 2 and raw.shape[1] >= 2:
+            surface = np.where(valid_mask, raw, np.nan)
+            gradient_y, gradient_x = np.gradient(
+                surface,
+                partition.resolution_y,
+                partition.resolution_x,
+            )
+            slope = np.degrees(np.arctan(np.hypot(gradient_x, gradient_y)))
+            slope_values = slope[np.isfinite(slope) & valid_mask]
+
+        def metric(key: str, label: str, value: float, unit: str, method: str) -> TerrainMetric:
+            return TerrainMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        prefix = "terrain" if product.surface_model == "terrain" else "surface"
+        elevation_label = (
+            "Gemiddelde maaiveldhoogte"
+            if product.surface_model == "terrain"
+            else "Gemiddelde oppervlaktehoogte"
+        )
+        metrics = [
+            metric(f"{prefix}_elevation_mean_m", elevation_label, values.mean(), "m TAW", "mean_valid_cells"),
+            metric(f"{prefix}_elevation_min_m", "Laagste hoogte", values.min(), "m TAW", "minimum_valid_cells"),
+            metric(f"{prefix}_elevation_max_m", "Hoogste hoogte", values.max(), "m TAW", "maximum_valid_cells"),
+            metric(f"{prefix}_elevation_p10_m", "10e percentiel hoogte", np.percentile(values, 10), "m TAW", "percentile_10_valid_cells"),
+            metric(f"{prefix}_elevation_p90_m", "90e percentiel hoogte", np.percentile(values, 90), "m TAW", "percentile_90_valid_cells"),
+            metric("relief_m", "Reliëfverschil", values.max() - values.min(), "m", "maximum_minus_minimum"),
+        ]
+        if slope_values.size:
+            metrics.extend(
+                [
+                    metric("slope_mean_deg", "Gemiddelde helling", slope_values.mean(), "°", "mean_finite_gradient"),
+                    metric("slope_p90_deg", "90e percentiel helling", np.percentile(slope_values, 90), "°", "percentile_90_finite_gradient"),
+                    metric("slope_max_deg", "Steilste helling", slope_values.max(), "°", "maximum_finite_gradient"),
+                ]
+            )
+        primary = metrics[0]
+        selected_cell_count = int(partition.selected_cells.sum())
+        first_dataset = partition.datasets[0]
+        response = TerrainSelectionResponse(
+            dataset_id=first_dataset.id,
+            dataset_ids=[dataset.id for dataset in partition.datasets],
+            partition_count=len(partition.datasets),
+            product_key=product.key,
+            surface_model=product.surface_model,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            sample_count=int(values.size),
+            slope_sample_count=int(slope_values.size),
+            coverage_ratio=round(float(values.size / max(1, selected_cell_count)), 6),
+            resolution_m=round(max(partition.resolution_x, partition.resolution_y), 4),
+            vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
+            summary=TerrainSelectionSummary(
+                metric_label=primary.metric_label,
+                metric_value=primary.metric_value,
+                metric_unit=primary.metric_unit,
+                aggregation_method=primary.aggregation_method,
+                primary_metric_key=primary.metric_key,
+                metrics=metrics,
+            ),
+            unsupported_metrics=TerrainAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=(
+                f"{TerrainAnalysisService.LIMITATION} De selectie werd exact berekend over "
+                f"{len(partition.datasets)} persistente gemeentelijke rasterpartities."
+            ),
+            generated_at=datetime.now(UTC).isoformat(),
+        )
+        return response.model_dump(mode="json")
+
     @staticmethod
     def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes:
         dataset = TerrainAnalysisService._load_dataset(db, project_id, dataset_id)