JensandClaude Opus 5
dd87a62e8f
Four ways a selection produced a confident number about a different area than the operator drew: Flood hazard divided the inundated cells by every cell in the drawn rectangle, including cells the VMM raster does not model at all. A selection reaching past the modelled extent therefore reported a diluted risk share, turning missing data into an implied absence of risk. Terrain, bathymetry and thematic raster already divided by valid cells; flood hazard was the outlier. It now reports the three populations separately, states model coverage next to the drawn area, and returns a null fraction rather than a zero when nothing was modelled. geometry_mask selects a cell when its centre falls inside the geometry, so a rectangle smaller than one cell — or one landing between four centres — selected nothing and the analysis returned zeros indistinguishable on screen from "we looked and there is nothing here". On a 100 m population raster a 40 m rectangle over a city block reported no inhabitants. Selection now falls back to the touched cells and says that it did, since the answer then covers more ground than was requested. rasterio.mask applies the same centre rule when cropping, so that call is widened too; the cells that count are still decided by the centre rule wherever it selects anything. The object count treated any feature touching the selection as whole, while intersection_area clipped it — two headline numbers on one panel describing different populations. The count stays whole-feature, which is what "objecten" means to an operator, but now reports how many the edge cuts and is marked an estimate when it does. The area_weighted_sum branch reuses that same count instead of issuing its own near-identical query. Partitioned selection de-duplicated the count on source_feature_id but returned the raw rows, so a building on a municipal boundary was counted once and drawn twice. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
368 lines
16 KiB
Python
368 lines
16 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.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.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,
|
|
)
|
|
# ``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
|