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