360 lines
15 KiB
Python
360 lines
15 KiB
Python
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.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.features import geometry_mask
|
|
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,
|
|
)
|
|
clipped, clipped_transform = mask(
|
|
source,
|
|
[mapping(analysis_geometry)],
|
|
crop=True,
|
|
filled=False,
|
|
indexes=[1],
|
|
)
|
|
band = np.ma.asarray(clipped[0], dtype="float64")
|
|
raw = band.filled(np.nan)
|
|
selected_cells = geometry_mask(
|
|
[mapping(analysis_geometry)],
|
|
out_shape=band.shape,
|
|
transform=clipped_transform,
|
|
invert=True,
|
|
)
|
|
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),
|
|
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
|