Initial public release

backend/app/services/bathymetry_raster_analysis_service.py

@@ -0,0 +1,366 @@
+from __future__ import annotations
+
+import io
+import math
+from datetime import UTC, datetime
+from pathlib import Path
+from uuid import UUID
+
+from geoalchemy2.shape import to_shape
+from pyproj import Transformer
+from shapely.geometry import box, mapping
+from shapely.ops import transform as shapely_transform
+
+from app.core.config import Settings, get_settings
+from app.core.errors import AppError
+from app.services.raster_cell_selection import select_cells
+from app.models import Area, Dataset
+from app.schemas.bathymetry import (
+    BathymetryRasterMetric,
+    BathymetryRasterSelectionRequest,
+    BathymetryRasterSelectionResponse,
+    BathymetryRasterSelectionSummary,
+)
+
+
+class BathymetryRasterAnalysisService:
+    SOURCE_NAME = "spw_bathymetry"
+    PRODUCT_KEY = "spw_bathymetry_50cm_mdng"
+    UNSUPPORTED_METRICS = [
+        "current_water_depth_m",
+        "water_volume_m3",
+        "vertical_datum_conversion",
+    ]
+    LIMITATION = (
+        "De rasterwaarden zijn waterbodemhoogtes in mDNG uit een samengestelde SPW-opmeting "
+        "(2019-2022). Zonder een gelijktijdig waterpeil zijn actuele waterdiepte en watervolume "
+        "niet berekenbaar. mDNG wordt niet stilzwijgend naar TAW, LAT of een ander verticaal datum omgezet."
+    )
+
+    @staticmethod
+    def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> Dataset:
+        dataset = db.get(Dataset, dataset_id)
+        if not dataset or dataset.project_id != project_id:
+            raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
+        if dataset.dataset_type != "raster" or dataset.source_name != BathymetryRasterAnalysisService.SOURCE_NAME:
+            raise AppError(
+                code="INVALID_BATHYMETRY_RASTER_DATASET",
+                message="Bathymetry analysis requires a governed SPW bathymetry raster dataset",
+                status_code=400,
+            )
+        if dataset.status != "ready" or not dataset.storage_path or not Path(dataset.storage_path).is_file():
+            raise AppError(
+                code="DATASET_FILE_MISSING",
+                message="Persisted bathymetry raster file is unavailable",
+                status_code=404,
+            )
+        return dataset
+
+    @staticmethod
+    def _metadata(dataset: Dataset) -> dict:
+        metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
+        if (
+            metadata.get("product_key") != BathymetryRasterAnalysisService.PRODUCT_KEY
+            or metadata.get("theme") != "bathymetry"
+            or metadata.get("value_semantics") != "bed_elevation"
+            or metadata.get("vertical_reference") != "mDNG"
+            or metadata.get("source_crs") != "EPSG:3812"
+        ):
+            raise AppError(
+                code="INVALID_BATHYMETRY_RASTER_METADATA",
+                message="Bathymetry raster provenance or value semantics are incomplete",
+                status_code=409,
+            )
+        return metadata
+
+    @staticmethod
+    def _selection_geometry(db, project_id: UUID, payload: BathymetryRasterSelectionRequest):
+        selection = box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)
+        if payload.area_id is None:
+            return selection
+        area = db.get(Area, payload.area_id)
+        if not area:
+            raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
+        if area.project_id != project_id:
+            raise AppError(
+                code="INVALID_DATASET_SCOPE",
+                message="Area does not belong to this project",
+                status_code=400,
+            )
+        selection = selection.intersection(to_shape(area.geometry))
+        if selection.is_empty or selection.area <= 0:
+            raise AppError(
+                code="BATHYMETRY_SELECTION_OUTSIDE_AREA",
+                message="Selection does not overlap the selected work area",
+                status_code=422,
+            )
+        return selection
+
+    @staticmethod
+    def analyze(
+        db,
+        project_id: UUID,
+        dataset_id: UUID,
+        payload: BathymetryRasterSelectionRequest,
+        *,
+        settings: Settings | None = None,
+    ) -> dict:
+        resolved_settings = settings or get_settings()
+        dataset = BathymetryRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        source_metadata = BathymetryRasterAnalysisService._metadata(dataset)
+        selection_4326 = BathymetryRasterAnalysisService._selection_geometry(db, project_id, payload)
+        try:
+            import numpy as np
+            import rasterio
+            from rasterio.mask import mask
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio and numpy are required for bathymetry analysis",
+                status_code=503,
+            ) from exc
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                if source.crs is None or source.crs.to_epsg() != 3812:
+                    raise AppError(
+                        code="INVALID_DATASET_CRS",
+                        message="SPW bathymetry raster CRS must be EPSG:3812",
+                        status_code=409,
+                    )
+                if source.count != 1:
+                    raise AppError(
+                        code="INVALID_BATHYMETRY_RASTER_BANDS",
+                        message="SPW bathymetry requires one bed-elevation band",
+                        status_code=409,
+                    )
+                transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True)
+                selection_metric = shapely_transform(transformer.transform, selection_4326)
+                analysis_geometry = selection_metric.intersection(box(*source.bounds))
+                if analysis_geometry.is_empty or analysis_geometry.area <= 0:
+                    raise AppError(
+                        code="BATHYMETRY_SELECTION_OUTSIDE_DATASET",
+                        message="Selection does not overlap the persisted bathymetry raster",
+                        status_code=422,
+                    )
+                min_x, min_y, max_x, max_y = analysis_geometry.bounds
+                expected_cells = math.ceil((max_x - min_x) / abs(source.res[0])) * math.ceil(
+                    (max_y - min_y) / abs(source.res[1])
+                )
+                if expected_cells > resolved_settings.bathymetry_raster_max_pixels:
+                    raise AppError(
+                        code="BATHYMETRY_SELECTION_TOO_LARGE",
+                        message="Bathymetry analysis exceeds the configured raster cell limit",
+                        details={
+                            "pixel_count": expected_cells,
+                            "max_pixels": resolved_settings.bathymetry_raster_max_pixels,
+                        },
+                        status_code=422,
+                    )
+                # ``all_touched`` keeps the values of cells the selection only
+                # clips, so a selection finer than one cell still has data to
+                # read. Which of those cells actually count is decided by
+                # ``select_cells`` below, so the normal result is unchanged.
+                clipped, clipped_transform = mask(
+                    source,
+                    [mapping(analysis_geometry)],
+                    crop=True,
+                    filled=False,
+                    indexes=[1],
+                    all_touched=True,
+                )
+                band = np.ma.asarray(clipped[0], dtype="float64")
+                raw = band.filled(np.nan)
+                cell_selection = select_cells(
+                    analysis_geometry,
+                    out_shape=band.shape,
+                    transform=clipped_transform,
+                    cell_area_m2=abs(float(source.res[0])) * abs(float(source.res[1])),
+                )
+                selected_cells = cell_selection.mask
+                valid_cells = selected_cells & ~np.ma.getmaskarray(band) & np.isfinite(raw)
+                if source.nodata is not None:
+                    valid_cells &= ~np.isclose(raw, float(source.nodata))
+                values = raw[valid_cells]
+                if values.size == 0:
+                    raise AppError(
+                        code="BATHYMETRY_NO_VALID_DATA",
+                        message="No surveyed waterbed cells occur in this selection",
+                        status_code=422,
+                    )
+                resolution_x = abs(float(source.res[0]))
+                resolution_y = abs(float(source.res[1]))
+                cell_area_m2 = resolution_x * resolution_y
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="BATHYMETRY_ANALYSIS_FAILED",
+                message="The persisted bathymetry raster could not be analysed",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc
+
+        def metric(key: str, label: str, value: float, unit: str, method: str) -> BathymetryRasterMetric:
+            return BathymetryRasterMetric(
+                metric_key=key,
+                metric_label=label,
+                metric_value=round(float(value), 4),
+                metric_unit=unit,
+                aggregation_method=method,
+            )
+
+        selected_cell_count = int(selected_cells.sum())
+        valid_cell_count = int(values.size)
+        vertical_unit = str(source_metadata["vertical_reference"])
+        coverage_ratio = valid_cell_count / max(1, selected_cell_count)
+        metrics = [
+            metric(
+                "bed_elevation_mean_m",
+                "Gemiddelde waterbodemhoogte",
+                values.mean(),
+                f"m {vertical_unit}",
+                "mean_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_min_m",
+                "Laagste waterbodemhoogte",
+                values.min(),
+                f"m {vertical_unit}",
+                "minimum_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_max_m",
+                "Hoogste waterbodemhoogte",
+                values.max(),
+                f"m {vertical_unit}",
+                "maximum_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_p10_m",
+                "10e percentiel waterbodemhoogte",
+                np.percentile(values, 10),
+                f"m {vertical_unit}",
+                "percentile_10_valid_source_cells",
+            ),
+            metric(
+                "bed_elevation_p90_m",
+                "90e percentiel waterbodemhoogte",
+                np.percentile(values, 90),
+                f"m {vertical_unit}",
+                "percentile_90_valid_source_cells",
+            ),
+            metric(
+                "surveyed_bed_surface_ha",
+                "Oppervlakte met gemeten waterbodem",
+                valid_cell_count * cell_area_m2 / 10_000.0,
+                "ha",
+                "valid_source_cells_times_cell_area",
+            ),
+            metric(
+                "bathymetry_coverage_pct",
+                "Dekking waterbodemmeting",
+                coverage_ratio * 100.0,
+                "%",
+                "valid_source_cells_divided_by_selected_cells",
+            ),
+        ]
+        primary = metrics[0]
+        response = BathymetryRasterSelectionResponse(
+            dataset_id=dataset.id,
+            product_key=BathymetryRasterAnalysisService.PRODUCT_KEY,
+            selection_bbox=payload.bbox,
+            selection_area_id=payload.area_id,
+            selected_cell_count=selected_cell_count,
+            valid_cell_count=valid_cell_count,
+            coverage_ratio=round(coverage_ratio, 6),
+            cell_selection_warning=cell_selection.warning,
+            resolution_m=round(max(resolution_x, resolution_y), 4),
+            vertical_reference=vertical_unit,
+            survey_period=str(source_metadata.get("survey_period") or "2019-2022"),
+            summary=BathymetryRasterSelectionSummary(
+                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=BathymetryRasterAnalysisService.UNSUPPORTED_METRICS,
+            limitation_message=BathymetryRasterAnalysisService.LIMITATION,
+            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 = BathymetryRasterAnalysisService._load_dataset(db, project_id, dataset_id)
+        BathymetryRasterAnalysisService._metadata(dataset)
+        try:
+            import numpy as np
+            import rasterio
+            from PIL import Image
+            from rasterio.enums import Resampling
+        except ImportError as exc:
+            raise AppError(
+                code="RASTER_PROCESSING_UNAVAILABLE",
+                message="Rasterio, numpy and Pillow are required for bathymetry rendering",
+                status_code=503,
+            ) from exc
+
+        try:
+            with rasterio.open(dataset.storage_path) as source:
+                scale = min(1.0, max_dimension / max(source.width, source.height))
+                width = max(1, round(source.width * scale))
+                height = max(1, round(source.height * scale))
+                data = source.read(
+                    1,
+                    out_shape=(height, width),
+                    masked=True,
+                    resampling=Resampling.bilinear,
+                )
+                values = np.asarray(data.filled(np.nan), dtype="float64")
+                valid = np.isfinite(values) & ~np.ma.getmaskarray(data)
+                if source.nodata is not None:
+                    valid &= ~np.isclose(values, float(source.nodata))
+                if not valid.any():
+                    raise AppError(
+                        code="BATHYMETRY_NO_VALID_DATA",
+                        message="Bathymetry raster contains no renderable cells",
+                        status_code=422,
+                    )
+                low, high = np.percentile(values[valid], [2, 98])
+                if high <= low:
+                    high = low + 1.0
+                normalized = np.clip((values - low) / (high - low), 0.0, 1.0)
+                normalized = np.where(valid, normalized, 0.0)
+                stops = np.asarray([0.0, 0.35, 0.7, 1.0])
+                colors = np.asarray(
+                    [
+                        [8, 47, 73],
+                        [15, 118, 140],
+                        [103, 190, 170],
+                        [236, 224, 163],
+                    ],
+                    dtype="float64",
+                )
+                rgba = np.zeros((height, width, 4), dtype="uint8")
+                for channel in range(3):
+                    rgba[:, :, channel] = np.interp(
+                        normalized,
+                        stops,
+                        colors[:, channel],
+                    ).astype("uint8")
+                rgba[:, :, 3] = np.where(valid, 220, 0).astype("uint8")
+                output = io.BytesIO()
+                Image.fromarray(rgba).save(output, format="PNG", optimize=True)
+                return output.getvalue()
+        except AppError:
+            raise
+        except Exception as exc:
+            raise AppError(
+                code="BATHYMETRY_PREVIEW_FAILED",
+                message="The persisted bathymetry raster could not be rendered",
+                details={"reason": str(exc)},
+                status_code=500,
+            ) from exc