diff --git a/backend/app/schemas/bathymetry.py b/backend/app/schemas/bathymetry.py index d2bbceda..5ff3b115 100644 --- a/backend/app/schemas/bathymetry.py +++ b/backend/app/schemas/bathymetry.py @@ -160,6 +160,10 @@ class BathymetryRasterSelectionResponse(BaseModel): selected_cell_count: int = Field(ge=1) valid_cell_count: int = Field(ge=1) coverage_ratio: float = Field(ge=0, le=1) + # Set when the drawn selection is smaller than one source cell and the + # analysis was widened to the cells it touches, so the value covers more + # ground than was requested. + cell_selection_warning: str | None = None resolution_m: float = Field(gt=0) vertical_reference: str survey_period: str diff --git a/backend/app/schemas/dhmv.py b/backend/app/schemas/dhmv.py index 74cc1a39..c00da18f 100644 --- a/backend/app/schemas/dhmv.py +++ b/backend/app/schemas/dhmv.py @@ -90,6 +90,10 @@ class TerrainSelectionResponse(BaseModel): sample_count: int slope_sample_count: int coverage_ratio: float + # Set when the drawn selection is smaller than one source cell and the + # analysis was widened to the cells it touches, so the value covers more + # ground than was requested. + cell_selection_warning: str | None = None resolution_m: float vertical_reference: str summary: TerrainSelectionSummary diff --git a/backend/app/schemas/flood_hazard.py b/backend/app/schemas/flood_hazard.py index 452e404c..8bc9bae2 100644 --- a/backend/app/schemas/flood_hazard.py +++ b/backend/app/schemas/flood_hazard.py @@ -93,9 +93,17 @@ class FloodHazardSelectionResponse(BaseModel): return_period_years: int selection_bbox: VectorSelectionBBox selection_area_id: UUID | None = None + # Three populations kept apart: cells drawn, cells the model covers, and + # cells with a positive modelled depth. ``inundated_fraction`` is a share + # of the modelled cells, and is null when nothing was modelled — absence + # of a model is not evidence of zero risk. selected_cell_count: int + valid_cell_count: int = 0 + no_data_cell_count: int = 0 + data_coverage_ratio: float = 1.0 inundated_cell_count: int - inundated_fraction: float + inundated_fraction: float | None = None + coverage_warning: str | None = None resolution_m: float summary: FloodHazardSelectionSummary unsupported_metrics: list[str] diff --git a/backend/app/schemas/operations.py b/backend/app/schemas/operations.py index fa4e0fe1..f3e70532 100644 --- a/backend/app/schemas/operations.py +++ b/backend/app/schemas/operations.py @@ -236,7 +236,13 @@ class VectorSelectionSummary(BaseModel): metric_unit: str aggregation_method: str primary_metric_key: str | None = None + # ``feature_count`` counts whole features that touch the selection, while + # area and length metrics clip to it. These fields say how far the two + # populations diverge, so the numbers on one panel can be read together. feature_count: int + fully_covered_feature_count: int | None = None + partially_covered_feature_count: int | None = None + selection_edge_warning: str | None = None is_estimate: bool = False warning: str | None = None metrics: list[VectorSelectionMetric] = Field(default_factory=list) diff --git a/backend/app/schemas/thematic_raster.py b/backend/app/schemas/thematic_raster.py index 405dfdb5..ce51cbe5 100644 --- a/backend/app/schemas/thematic_raster.py +++ b/backend/app/schemas/thematic_raster.py @@ -113,6 +113,10 @@ class ThematicRasterSelectionResponse(BaseModel): selected_cell_count: int valid_cell_count: int coverage_ratio: float + # Set when the drawn selection is smaller than one source cell and the + # analysis was widened to the cells it touches, so the value covers more + # ground than was requested. + cell_selection_warning: str | None = None resolution_m: float observation_year: int summary: ThematicRasterSelectionSummary diff --git a/backend/app/services/bathymetry_raster_analysis_service.py b/backend/app/services/bathymetry_raster_analysis_service.py index 28d7b0ca..bb7e1612 100644 --- a/backend/app/services/bathymetry_raster_analysis_service.py +++ b/backend/app/services/bathymetry_raster_analysis_service.py @@ -13,6 +13,7 @@ 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, @@ -157,21 +158,27 @@ class BathymetryRasterAnalysisService: }, 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) - selected_cells = geometry_mask( - [mapping(analysis_geometry)], + cell_selection = select_cells( + analysis_geometry, out_shape=band.shape, transform=clipped_transform, - invert=True, + 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)) @@ -268,6 +275,7 @@ class BathymetryRasterAnalysisService: 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"), diff --git a/backend/app/services/flood_hazard_analysis_service.py b/backend/app/services/flood_hazard_analysis_service.py index 4ba2011b..570bca88 100644 --- a/backend/app/services/flood_hazard_analysis_service.py +++ b/backend/app/services/flood_hazard_analysis_service.py @@ -2,8 +2,10 @@ from __future__ import annotations import io import math +from dataclasses import dataclass from datetime import UTC, datetime from pathlib import Path +from typing import Any from uuid import UUID from geoalchemy2.shape import to_shape @@ -13,6 +15,7 @@ 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.flood_hazard import ( FloodHazardMetric, @@ -25,6 +28,83 @@ from app.services.flood_hazard_acquisition_service import FloodHazardAcquisition from app.services.raster_partition_analysis_service import RasterPartitionAnalysisService +@dataclass(frozen=True) +class FloodHazardCellStatistics: + """Cell populations behind one flood-hazard selection. + + Three populations, deliberately kept apart: + + ``selected`` + every cell whose centre falls inside the drawn selection; + ``valid`` + the subset the VMM raster actually models — finite, not nodata; + ``inundated`` + the subset of valid cells with a positive modelled depth. + + Risk is a share of what was modelled. Dividing by the selected cells + instead silently reports "no data" as "no risk", which for a selection + reaching past the modelled extent understates the hazard by whatever + fraction of the rectangle the model never covered. + """ + + selected_cell_count: int + valid_cell_count: int + inundated_cell_count: int + depth_values: Any + + @property + def no_data_cell_count(self) -> int: + return max(0, self.selected_cell_count - self.valid_cell_count) + + @property + def data_coverage_ratio(self) -> float: + if self.selected_cell_count <= 0: + return 0.0 + return self.valid_cell_count / self.selected_cell_count + + @property + def inundated_fraction(self) -> float | None: + """``None`` when nothing was modelled: absence of data is not a zero.""" + + if self.valid_cell_count <= 0: + return None + return self.inundated_cell_count / self.valid_cell_count + + def inundated_area_ha(self, cell_area_m2: float) -> float: + return self.inundated_cell_count * cell_area_m2 / 10_000.0 + + def analysed_area_ha(self, cell_area_m2: float) -> float: + """Area the model actually covers inside the selection.""" + + return self.valid_cell_count * cell_area_m2 / 10_000.0 + + def selected_area_ha(self, cell_area_m2: float) -> float: + """Area of the selection as rasterised, model coverage aside.""" + + return self.selected_cell_count * cell_area_m2 / 10_000.0 + + @classmethod + def from_cells(cls, values: Any, selected: Any, *, nodata: float | None) -> "FloodHazardCellStatistics": + import numpy as np + + raw = np.asarray(values, dtype="float64") + selected_mask = np.asarray(selected, dtype=bool) + + has_data = selected_mask & np.isfinite(raw) + if nodata is not None: + has_data &= ~np.isclose(raw, float(nodata)) + # A modelled zero or negative depth is data: it says "dry here", which + # is a different statement from "not modelled here". + inundated = has_data & (raw > 0.0) + + return cls( + selected_cell_count=int(selected_mask.sum()), + valid_cell_count=int(has_data.sum()), + inundated_cell_count=int(inundated.sum()), + depth_values=raw[inundated], + ) + + class FloodHazardAnalysisService: UNSUPPORTED_METRICS = [ "bathymetry_depth_m", @@ -36,6 +116,77 @@ class FloodHazardAnalysisService: "gemodelleerde maxima op en is geen gelijktijdig opgeslagen watervolume, actuele waterstand of bathymetrie." ) + @staticmethod + def _coverage_metrics(stats: "FloodHazardCellStatistics", cell_area_m2: float, metric) -> list[FloodHazardMetric]: + """Headline metrics, each stating which population it is a share of. + + The analysed area is reported next to the drawn area so an operator can + see immediately how much of the rectangle the flood model covers. A + selection with no model data reports 0% coverage rather than 0% risk. + """ + + metrics = [ + metric( + "modelled_inundated_area_ha", + "Gemodelleerd overstroomd oppervlak", + stats.inundated_area_ha(cell_area_m2), + "ha", + "positive_depth_cells_times_cell_area", + ), + metric( + "modelled_inundated_share_pct", + "Aandeel gemodelleerd gebied met diepte", + 0.0 if stats.inundated_fraction is None else stats.inundated_fraction * 100.0, + "%", + "positive_depth_cells_divided_by_modelled_cells", + ), + metric( + "modelled_area_ha", + "Oppervlak met overstromingsmodel", + stats.analysed_area_ha(cell_area_m2), + "ha", + "modelled_cells_times_cell_area", + ), + metric( + "selection_area_ha", + "Oppervlak van de selectie", + stats.selected_area_ha(cell_area_m2), + "ha", + "selected_cells_times_cell_area", + ), + metric( + "model_coverage_pct", + "Deel van de selectie met een model", + stats.data_coverage_ratio * 100.0, + "%", + "modelled_cells_divided_by_selected_cells", + ), + ] + return metrics + + @staticmethod + def _combined_warning(stats: "FloodHazardCellStatistics", cell_selection_warning: str | None) -> str | None: + parts = [ + part + for part in (cell_selection_warning, FloodHazardAnalysisService._coverage_warning(stats)) + if part + ] + return " ".join(parts) if parts else None + + @staticmethod + def _coverage_warning(stats: "FloodHazardCellStatistics") -> str | None: + if stats.valid_cell_count <= 0: + return ( + "Voor deze selectie bestaat geen VMM-overstromingsmodel. Er is dus geen overstromingsrisico " + "gemeten; dit is geen bevestiging dat het risico nul is." + ) + if stats.data_coverage_ratio < 0.999: + return ( + f"Het VMM-model dekt {stats.data_coverage_ratio * 100:.1f}% van deze selectie. Percentages gelden " + "voor het gemodelleerde deel, niet voor de volledige selectie." + ) + return None + @staticmethod def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> Dataset: dataset = db.get(Dataset, dataset_id) @@ -110,17 +261,32 @@ class FloodHazardAnalysisService: details={"pixel_count": expected_cells, "max_pixels": resolved_settings.flood_hazard_max_pixels}, status_code=422, ) - clipped, clipped_transform = mask(source, [mapping(analysis_geometry)], crop=True, filled=False, indexes=[1]) + # ``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, + ) depth = np.ma.asarray(clipped[0], dtype="float64") raw = depth.filled(np.nan) - selected_cells = geometry_mask([mapping(analysis_geometry)], out_shape=depth.shape, transform=clipped_transform, invert=True) - nodata = source.nodata - valid = selected_cells & ~np.ma.getmaskarray(depth) & np.isfinite(raw) & (raw > 0.0) - if nodata is not None: - valid &= raw != float(nodata) - values = raw[valid] - selected_cell_count = int(selected_cells.sum()) - inundated_cell_count = int(values.size) + cell_selection = select_cells( + analysis_geometry, + out_shape=depth.shape, + transform=clipped_transform, + cell_area_m2=abs(float(source.res[0])) * abs(float(source.res[1])), + ) + selected_cells = cell_selection.mask + # A masked cell carries no model value, so fold the mask into + # the raw array before the populations are separated. + raw = np.where(np.ma.getmaskarray(depth), np.nan, raw) + stats = FloodHazardCellStatistics.from_cells(raw, selected_cells, nodata=source.nodata) + values = stats.depth_values resolution_x = abs(float(source.res[0])) resolution_y = abs(float(source.res[1])) cell_area_m2 = resolution_x * resolution_y @@ -143,18 +309,8 @@ class FloodHazardAnalysisService: 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 = FloodHazardAnalysisService._coverage_metrics(stats, cell_area_m2, metric) + if stats.inundated_cell_count: metrics.extend( [ metric("modelled_depth_mean_m", "Gemiddelde gemodelleerde maximumdiepte", values.mean(), "m", "mean_positive_depth_cells"), @@ -181,9 +337,14 @@ class FloodHazardAnalysisService: 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), + selected_cell_count=stats.selected_cell_count, + valid_cell_count=stats.valid_cell_count, + no_data_cell_count=stats.no_data_cell_count, + data_coverage_ratio=round(stats.data_coverage_ratio, 6), + inundated_cell_count=stats.inundated_cell_count, + inundated_fraction=( + None if stats.inundated_fraction is None else round(stats.inundated_fraction, 6) + ), resolution_m=round(max(resolution_x, resolution_y), 4), summary=FloodHazardSelectionSummary( metric_label=primary.metric_label, @@ -193,6 +354,7 @@ class FloodHazardAnalysisService: primary_metric_key=primary.metric_key, metrics=metrics, ), + coverage_warning=FloodHazardAnalysisService._combined_warning(stats, cell_selection.warning), unsupported_metrics=FloodHazardAnalysisService.UNSUPPORTED_METRICS, limitation_message=FloodHazardAnalysisService.LIMITATION, generated_at=datetime.now(UTC).isoformat(), @@ -236,16 +398,12 @@ class FloodHazardAnalysisService: status_code=503, ) from exc - raw = partition.values - valid = ( - partition.selected_cells - & np.isfinite(raw) - & (raw != FloodHazardAcquisitionService.NODATA) - & (raw > 0.0) + stats = FloodHazardCellStatistics.from_cells( + partition.values, + partition.selected_cells, + nodata=FloodHazardAcquisitionService.NODATA, ) - values = raw[valid] - selected_cell_count = int(partition.selected_cells.sum()) - inundated_cell_count = int(values.size) + values = stats.depth_values cell_area_m2 = partition.resolution_x * partition.resolution_y def metric(key: str, label: str, value: float, unit: str, method: str) -> FloodHazardMetric: @@ -257,24 +415,8 @@ class FloodHazardAnalysisService: 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 = FloodHazardAnalysisService._coverage_metrics(stats, cell_area_m2, metric) + if stats.inundated_cell_count: metrics.extend( [ metric("modelled_depth_mean_m", "Gemiddelde gemodelleerde maximumdiepte", values.mean(), "m", "mean_positive_depth_cells"), @@ -302,9 +444,14 @@ class FloodHazardAnalysisService: 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), + selected_cell_count=stats.selected_cell_count, + valid_cell_count=stats.valid_cell_count, + no_data_cell_count=stats.no_data_cell_count, + data_coverage_ratio=round(stats.data_coverage_ratio, 6), + inundated_cell_count=stats.inundated_cell_count, + inundated_fraction=( + None if stats.inundated_fraction is None else round(stats.inundated_fraction, 6) + ), resolution_m=round(max(partition.resolution_x, partition.resolution_y), 4), summary=FloodHazardSelectionSummary( metric_label=primary.metric_label, @@ -314,6 +461,7 @@ class FloodHazardAnalysisService: primary_metric_key=primary.metric_key, metrics=metrics, ), + coverage_warning=FloodHazardAnalysisService._combined_warning(stats, partition.cell_selection_warning), unsupported_metrics=FloodHazardAnalysisService.UNSUPPORTED_METRICS, limitation_message=( f"{FloodHazardAnalysisService.LIMITATION} De selectie werd exact berekend over " diff --git a/backend/app/services/raster_cell_selection.py b/backend/app/services/raster_cell_selection.py new file mode 100644 index 00000000..65f66ab6 --- /dev/null +++ b/backend/app/services/raster_cell_selection.py @@ -0,0 +1,85 @@ +"""Choosing which raster cells a drawn selection covers. + +``rasterio.features.geometry_mask`` selects a cell when the cell *centre* falls +inside the geometry. That is the right rule for a selection spanning many +cells, and the wrong one for a small selection: a rectangle smaller than a cell, +or one landing between four centres, selects nothing. The analysis then reports +zeros, which on screen reads as "we measured this area and found nothing" +rather than "this selection is finer than the source raster". + +Falling back to every touched cell keeps a small selection answerable, at the +cost of analysing more ground than was drawn. That trade is only honest if it +is stated, so the fallback is reported alongside the result. +""" + +from __future__ import annotations + +from dataclasses import dataclass +from typing import Any + + +@dataclass(frozen=True) +class RasterCellSelection: + mask: Any + mode: str + expanded_to_touched_cells: bool + warning: str | None + + +def select_cells( + geometry, + *, + out_shape: tuple[int, int], + transform, + cell_area_m2: float | None = None, +) -> RasterCellSelection: + """Mask the cells a selection covers, widening only when it covers none. + + ``geometry`` must already be in the raster's own CRS. + """ + + from rasterio.features import geometry_mask + from shapely.geometry import mapping + + shapes = [mapping(geometry)] + mask = geometry_mask(shapes, out_shape=out_shape, transform=transform, invert=True) + if mask.any(): + return RasterCellSelection( + mask=mask, + mode="cell_centre", + expanded_to_touched_cells=False, + warning=None, + ) + + touched = geometry_mask( + shapes, + out_shape=out_shape, + transform=transform, + invert=True, + all_touched=True, + ) + if not touched.any(): + # The selection does not reach the raster at all. Widening the rule + # must never manufacture coverage that is genuinely absent. + return RasterCellSelection( + mask=mask, + mode="cell_centre", + expanded_to_touched_cells=False, + warning=None, + ) + + cell_count = int(touched.sum()) + if cell_area_m2: + analysed_ha = cell_count * float(cell_area_m2) / 10_000.0 + detail = f"{cell_count} rastercel{'' if cell_count == 1 else 'len'} ({analysed_ha:.1f} ha)" + else: + detail = f"{cell_count} rastercel{'' if cell_count == 1 else 'len'}" + return RasterCellSelection( + mask=touched, + mode="all_touched", + expanded_to_touched_cells=True, + warning=( + f"De selectie is kleiner dan één rastercel van deze bron. Het resultaat geldt voor {detail} " + "die de selectie raken, dus voor een groter gebied dan getekend." + ), + ) diff --git a/backend/app/services/raster_partition_analysis_service.py b/backend/app/services/raster_partition_analysis_service.py index acdbd21d..a81df15b 100644 --- a/backend/app/services/raster_partition_analysis_service.py +++ b/backend/app/services/raster_partition_analysis_service.py @@ -12,6 +12,7 @@ from shapely.geometry import mapping from shapely.ops import transform as shapely_transform from app.core.errors import AppError +from app.services.raster_cell_selection import select_cells from app.models import Dataset @@ -22,6 +23,9 @@ class RasterPartitionSelection: selected_cells: Any resolution_x: float resolution_y: float + # Set when the selection is finer than one source cell and the analysis was + # widened to the touched cells, so it covers more ground than was drawn. + cell_selection_warning: str | None = None class RasterPartitionAnalysisService: @@ -185,18 +189,20 @@ class RasterPartitionAnalysisService: dtype="float32", ) values = np.asarray(mosaic[0], dtype="float64") - selected_cells = geometry_mask( - [mapping(selection_metric)], + cell_selection = select_cells( + selection_metric, out_shape=values.shape, transform=transform, - invert=True, + cell_area_m2=target_resolution * target_resolution, ) + selected_cells = cell_selection.mask return RasterPartitionSelection( datasets=datasets, values=values, selected_cells=selected_cells, resolution_x=target_resolution, resolution_y=target_resolution, + cell_selection_warning=cell_selection.warning, ) except AppError: raise diff --git a/backend/app/services/terrain_analysis_service.py b/backend/app/services/terrain_analysis_service.py index 5cb12f19..3e096368 100644 --- a/backend/app/services/terrain_analysis_service.py +++ b/backend/app/services/terrain_analysis_service.py @@ -13,6 +13,7 @@ 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.dhmv import ( TerrainMetric, @@ -171,12 +172,17 @@ class TerrainAnalysisService: }, 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, ) elevation = np.ma.asarray(clipped[0], dtype="float64") raw = elevation.filled(np.nan) @@ -184,12 +190,13 @@ class TerrainAnalysisService: invalid = ~np.isfinite(raw) if nodata is not None: invalid |= raw == float(nodata) - selected_cells = geometry_mask( - [mapping(analysis_geometry)], + cell_selection = select_cells( + analysis_geometry, out_shape=elevation.shape, transform=clipped_transform, - invert=True, + cell_area_m2=abs(float(source.res[0])) * abs(float(source.res[1])), ) + selected_cells = cell_selection.mask valid_mask = selected_cells & ~np.ma.getmaskarray(elevation) & ~invalid values = raw[valid_mask] if values.size == 0: @@ -319,6 +326,7 @@ class TerrainAnalysisService: sample_count=int(values.size), slope_sample_count=int(slope_values.size), coverage_ratio=round(float(values.size / max(1, selected_cell_count)), 6), + cell_selection_warning=cell_selection.warning, resolution_m=round(max(resolution_x, resolution_y), 4), vertical_reference=str( source_metadata.get("vertical_reference") @@ -516,6 +524,7 @@ class TerrainAnalysisService: sample_count=int(values.size), slope_sample_count=int(slope_values.size), coverage_ratio=round(float(values.size / max(1, selected_cell_count)), 6), + cell_selection_warning=partition.cell_selection_warning, resolution_m=round(max(partition.resolution_x, partition.resolution_y), 4), vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE, summary=TerrainSelectionSummary( diff --git a/backend/app/services/thematic_raster_analysis_service.py b/backend/app/services/thematic_raster_analysis_service.py index 62f9e2f7..d916486f 100644 --- a/backend/app/services/thematic_raster_analysis_service.py +++ b/backend/app/services/thematic_raster_analysis_service.py @@ -13,6 +13,7 @@ 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.thematic_raster import ( ThematicRasterMetric, @@ -118,10 +119,27 @@ class ThematicRasterAnalysisService: details={"pixel_count": expected_cells, "max_pixels": resolved_settings.thematic_raster_max_pixels}, status_code=422, ) - clipped, clipped_transform = mask(source, [mapping(analysis_geometry)], crop=True, filled=False, indexes=[1]) + # ``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) - selected = geometry_mask([mapping(analysis_geometry)], out_shape=band.shape, transform=clipped_transform, invert=True) + 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 = cell_selection.mask valid = selected & ~np.ma.getmaskarray(band) & np.isfinite(raw) if source.nodata is not None: valid &= ~np.isclose(raw, float(source.nodata)) @@ -195,6 +213,7 @@ class ThematicRasterAnalysisService: selected_cell_count=selected_cell_count, valid_cell_count=valid_cell_count, coverage_ratio=round(valid_cell_count / max(1, selected_cell_count), 6), + cell_selection_warning=cell_selection.warning, resolution_m=round(max(resolution_x, resolution_y), 4), observation_year=product.observation_year, summary=ThematicRasterSelectionSummary( diff --git a/backend/app/services/vector_feature_service.py b/backend/app/services/vector_feature_service.py index 4e2e7712..b47e11af 100644 --- a/backend/app/services/vector_feature_service.py +++ b/backend/app/services/vector_feature_service.py @@ -220,6 +220,73 @@ class VectorFeatureService: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} return isinstance(source_metadata.get("selection_aggregation"), dict) or VectorFeatureService._dataset_theme(dataset) is not None + @staticmethod + def deduplicate_rows(rows: list[Any]) -> list[Any]: + """Collapse rows that describe one source feature across partitions. + + Municipal partitions of one product overlap at their shared boundary, + so a rectangle drawn across it returns the same building from both. + An empty or missing ``source_feature_id`` is not a shared identity — + two rows without one are two features, not a duplicate pair. + """ + + seen: set[str] = set() + kept: list[Any] = [] + for row in rows: + source_feature_id = getattr(row, "source_feature_id", None) + identity = str(source_feature_id).strip() if source_feature_id is not None else "" + if not identity: + kept.append(row) + continue + if identity in seen: + continue + seen.add(identity) + kept.append(row) + return kept + + @staticmethod + def count_disclosure( + *, + total_feature_count: int, + fully_covered_feature_count: int | None, + ) -> dict[str, Any]: + """Describe how much of the counted population the selection cuts. + + A feature that merely touches the drawn rectangle is counted whole, + while ``intersection_area`` clips it. Reporting both numbers without + saying so puts two figures for different populations side by side. The + count stays whole-feature — that is what an operator expects from + "objecten" — but says how many of them the edge cuts, and is flagged as + an estimate when it does. + + ``fully_covered_feature_count`` is ``None`` when the selection covers a + pre-clipped whole work area, where no edge effect exists. + """ + + if fully_covered_feature_count is None: + return { + "partially_covered_feature_count": None, + "is_estimate": False, + "warning": None, + } + + partial = max(0, int(total_feature_count) - int(fully_covered_feature_count)) + if partial <= 0: + return { + "partially_covered_feature_count": 0, + "is_estimate": False, + "warning": None, + } + return { + "partially_covered_feature_count": partial, + "is_estimate": True, + "warning": ( + f"{partial} van de {int(total_feature_count)} objecten liggen deels buiten de selectie en zijn " + "aan de rand doorgesneden. Ze tellen volledig mee in het aantal; oppervlakte- en lengtematen " + "gebruiken alleen het deel binnen de selectie." + ), + } + @staticmethod def constrain_bbox_to_area( bbox: dict[str, Any], @@ -696,6 +763,11 @@ class VectorFeatureService: .limit(safe_limit + 1) .all() ) + if deduplicate_source_features: + # ``total_feature_count`` is already distinct; without this the map + # would draw a boundary feature once per partition and the returned + # count would exceed the headline number beside it. + rows = VectorFeatureService.deduplicate_rows(rows) truncated = total_feature_count > safe_limit selected_rows = rows[:safe_limit] features = [VectorFeatureService._row_to_geojson_feature(row) for row in selected_rows] @@ -767,6 +839,27 @@ class VectorFeatureService: if feature_count is None: feature_count = int(db.query(func.count(VectorFeature.id)).filter(*selection_filter).scalar() or 0) + # How many of the counted features the selection edge cuts. Skipped for + # a pre-clipped whole-area selection, which has no edge to cut against. + fully_covered_feature_count: int | None = None + if not full_dataset_area and feature_count: + try: + fully_covered_feature_count = int( + db.query(func.count(VectorFeature.id)) + .filter(*selection_filter) + .filter(func.ST_CoveredBy(VectorFeature.geometry, selection_shape)) + .scalar() + or 0 + ) + except Exception: + # Lightweight unit-test sessions do not implement every spatial + # predicate; the count then simply carries no edge disclosure. + fully_covered_feature_count = None + count_disclosure = VectorFeatureService.count_disclosure( + total_feature_count=feature_count, + fully_covered_feature_count=fully_covered_feature_count, + ) + source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} config = source_metadata.get("selection_aggregation") if not isinstance(config, dict): @@ -834,9 +927,17 @@ class VectorFeatureService: selection_shape=selection_shape, feature_count=feature_count, full_dataset_area=selection_is_preclipped, + partially_covered_feature_count=count_disclosure["partially_covered_feature_count"], ) for metric_config in metric_configs ] + # The whole-feature count carries the edge disclosure; a metric that + # already clips to the selection (area, length) does not need it. + for computed_metric in metrics: + if computed_metric["aggregation_method"] == "feature_count" and count_disclosure["is_estimate"]: + computed_metric["is_estimate"] = True + computed_metric["warning"] = computed_metric.get("warning") or count_disclosure["warning"] + primary_metric = metrics[0] return { "metric_label": primary_metric["metric_label"], @@ -845,6 +946,9 @@ class VectorFeatureService: "aggregation_method": primary_metric["aggregation_method"], "primary_metric_key": primary_metric["metric_key"], "feature_count": feature_count, + "fully_covered_feature_count": fully_covered_feature_count, + "partially_covered_feature_count": count_disclosure["partially_covered_feature_count"], + "selection_edge_warning": count_disclosure["warning"], "is_estimate": primary_metric["is_estimate"], "warning": primary_metric.get("warning"), "metrics": metrics, @@ -860,6 +964,7 @@ class VectorFeatureService: selection_shape: Any, feature_count: int, full_dataset_area: bool, + partially_covered_feature_count: int | None = None, ) -> dict[str, Any]: method = str(config.get("method") or "feature_count") unit = str(config.get("unit") or "objecten") @@ -950,12 +1055,16 @@ class VectorFeatureService: ) metric_value = float(aggregate_value or 0.0) if method == "area_weighted_sum" and not full_dataset_area: - partial_feature_count = ( - db.query(func.count(VectorFeature.id)) - .filter(*metric_filter) - .filter(~covered_by_selection) - .scalar() - ) + # The selection-edge count was already established for the + # feature count; a second query would ask the same question. + partial_feature_count = partially_covered_feature_count + if partial_feature_count is None: + partial_feature_count = ( + db.query(func.count(VectorFeature.id)) + .filter(*metric_filter) + .filter(~covered_by_selection) + .scalar() + ) is_estimate = bool(config.get("is_estimate", False)) or bool(partial_feature_count) if not is_estimate and config.get("warning_only_when_estimate", True): warning = None diff --git a/backend/tests/test_flood_hazard_selection_data_coverage.py b/backend/tests/test_flood_hazard_selection_data_coverage.py new file mode 100644 index 00000000..8f8c0eea --- /dev/null +++ b/backend/tests/test_flood_hazard_selection_data_coverage.py @@ -0,0 +1,103 @@ +"""Flood risk must be a share of what was modelled, not of what was drawn. + +The share and fraction divided the inundated cells by every cell whose centre +fell inside the selection, including cells where the VMM raster holds nodata +because the area lies outside the modelled extent. An operator drawing a +rectangle that reaches past the model coverage read "3% at risk" where the +honest answer is "of the 40% we have a model for, 7.5% is at risk, and for the +rest there is no model at all". + +Terrain, bathymetry and thematic raster analysis already divide by valid cells +and report a coverage ratio; this brings flood hazard in line. +""" + +from __future__ import annotations + +import pytest + +np = pytest.importorskip("numpy") + +from app.services.flood_hazard_analysis_service import FloodHazardCellStatistics + + +NODATA = -9999.0 + + +def _stats(values, selected) -> FloodHazardCellStatistics: + return FloodHazardCellStatistics.from_cells( + np.asarray(values, dtype="float64"), + np.asarray(selected, dtype=bool), + nodata=NODATA, + ) + + +def test_share_ignores_cells_the_model_does_not_cover() -> None: + # Ten selected cells: four modelled (one of them wet), six nodata. + values = [1.5, 0.0, 0.0, 0.0] + [NODATA] * 6 + selected = [True] * 10 + + stats = _stats(values, selected) + + assert stats.selected_cell_count == 10 + assert stats.valid_cell_count == 4 + assert stats.no_data_cell_count == 6 + assert stats.inundated_cell_count == 1 + # 1 of 4 modelled cells, not 1 of 10 drawn cells. + assert stats.inundated_fraction == pytest.approx(0.25) + assert stats.data_coverage_ratio == pytest.approx(0.4) + + +def test_cells_outside_the_drawn_selection_are_not_counted() -> None: + values = [1.5, 1.5, 0.0, 0.0] + selected = [True, False, True, False] + + stats = _stats(values, selected) + + assert stats.selected_cell_count == 2 + assert stats.valid_cell_count == 2 + assert stats.inundated_cell_count == 1 + assert stats.inundated_fraction == pytest.approx(0.5) + + +def test_a_selection_without_any_model_data_reports_zero_coverage() -> None: + stats = _stats([NODATA] * 4, [True] * 4) + + assert stats.valid_cell_count == 0 + assert stats.no_data_cell_count == 4 + assert stats.data_coverage_ratio == 0.0 + # No model, so no risk figure may be invented. + assert stats.inundated_fraction is None + + +def test_nan_is_treated_as_missing_model_data() -> None: + stats = _stats([float("nan"), 2.0], [True, True]) + + assert stats.valid_cell_count == 1 + assert stats.no_data_cell_count == 1 + assert stats.inundated_cell_count == 1 + + +def test_negative_depths_are_data_but_not_inundation() -> None: + """A modelled zero or negative depth means dry, not unknown.""" + + stats = _stats([0.0, 0.0, 3.0], [True, True, True]) + + assert stats.valid_cell_count == 3 + assert stats.inundated_cell_count == 1 + assert stats.inundated_fraction == pytest.approx(1 / 3) + + +def test_depth_statistics_use_only_inundated_cells() -> None: + stats = _stats([0.0, 2.0, 4.0, NODATA], [True] * 4) + + assert stats.depth_values.tolist() == [2.0, 4.0] + assert stats.depth_values.mean() == pytest.approx(3.0) + + +def test_areas_are_derived_from_the_matching_cell_populations() -> None: + stats = _stats([1.0, 1.0, 0.0, NODATA], [True] * 4) + + # 100 m2 cells: 2 inundated, 3 modelled, 4 drawn. + assert stats.inundated_area_ha(100.0) == pytest.approx(2 * 100.0 / 10_000.0) + assert stats.analysed_area_ha(100.0) == pytest.approx(3 * 100.0 / 10_000.0) + assert stats.selected_area_ha(100.0) == pytest.approx(4 * 100.0 / 10_000.0) diff --git a/backend/tests/test_partition_selection_deduplication.py b/backend/tests/test_partition_selection_deduplication.py new file mode 100644 index 00000000..9c386204 --- /dev/null +++ b/backend/tests/test_partition_selection_deduplication.py @@ -0,0 +1,61 @@ +"""A rectangle across a municipal boundary must not return the same object twice. + +Partitioned selection de-duplicated ``total_feature_count`` on +``source_feature_id`` but returned the raw rows. A feature present in two +municipal partitions was therefore drawn twice on the map and counted once in +the headline, so the number on the panel disagreed with the geometry beside it. +""" + +from __future__ import annotations + +from uuid import uuid4 + +from app.services.vector_feature_service import VectorFeatureService + + +class _Row: + def __init__(self, source_feature_id, row_id=None, dataset_id=None): + self.source_feature_id = source_feature_id + self.id = row_id or uuid4() + self.dataset_id = dataset_id or uuid4() + + +def _ids(rows): + return [row.source_feature_id or str(row.id) for row in rows] + + +def test_a_feature_in_two_partitions_is_returned_once() -> None: + shared = "grb-building-42" + rows = [_Row(shared), _Row("grb-building-7"), _Row(shared)] + + kept = VectorFeatureService.deduplicate_rows(rows) + + assert _ids(kept) == [shared, "grb-building-7"] + + +def test_the_first_occurrence_wins_so_the_result_is_stable() -> None: + first = _Row("dup") + second = _Row("dup") + + assert VectorFeatureService.deduplicate_rows([first, second])[0] is first + assert VectorFeatureService.deduplicate_rows([second, first])[0] is second + + +def test_rows_without_a_source_id_fall_back_to_their_own_identity() -> None: + """Two distinct rows with no source id are two distinct features.""" + + rows = [_Row(None), _Row(None)] + + assert len(VectorFeatureService.deduplicate_rows(rows)) == 2 + + +def test_an_empty_source_id_is_not_treated_as_a_shared_identity() -> None: + rows = [_Row(""), _Row("")] + + assert len(VectorFeatureService.deduplicate_rows(rows)) == 2 + + +def test_deduplication_leaves_a_clean_population_untouched() -> None: + rows = [_Row("a"), _Row("b"), _Row("c")] + + assert VectorFeatureService.deduplicate_rows(rows) == rows diff --git a/backend/tests/test_raster_cell_selection.py b/backend/tests/test_raster_cell_selection.py new file mode 100644 index 00000000..3c81c23c --- /dev/null +++ b/backend/tests/test_raster_cell_selection.py @@ -0,0 +1,84 @@ +"""A selection smaller than a raster cell must not silently read as zero. + +``geometry_mask`` selects a cell when its *centre* falls inside the geometry. +A rectangle smaller than one cell, or one that lands between four centres, +therefore selects nothing at all — and the analysis returned zeros, which on +screen is indistinguishable 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 every touched cell and says that it did, so the +value is readable as "at least one whole cell", not as an empty area. +""" + +from __future__ import annotations + +import pytest + +np = pytest.importorskip("numpy") +rasterio = pytest.importorskip("rasterio") + +from rasterio.transform import from_origin +from shapely.geometry import box + +from app.services.raster_cell_selection import select_cells + + +# 100 m cells, origin at the top-left corner of a 3x3 grid. +TRANSFORM = from_origin(200_000, 210_000, 100.0, 100.0) +SHAPE = (3, 3) + + +def test_a_normal_selection_uses_cell_centres() -> None: + selection = select_cells(box(200_000, 209_700, 200_300, 210_000), out_shape=SHAPE, transform=TRANSFORM) + + assert selection.mask.sum() == 9 + assert selection.mode == "cell_centre" + assert selection.expanded_to_touched_cells is False + assert selection.warning is None + + +def test_a_rectangle_smaller_than_one_cell_still_returns_that_cell() -> None: + selection = select_cells(box(200_010, 209_960, 200_050, 209_990), out_shape=SHAPE, transform=TRANSFORM) + + assert selection.mask.sum() == 1 + assert selection.mode == "all_touched" + assert selection.expanded_to_touched_cells is True + assert "cel" in selection.warning + + +def test_a_rectangle_between_four_cell_centres_returns_all_four() -> None: + selection = select_cells(box(200_080, 209_880, 200_120, 209_920), out_shape=SHAPE, transform=TRANSFORM) + + assert selection.mask.sum() == 4 + assert selection.expanded_to_touched_cells is True + + +def test_a_selection_entirely_off_the_raster_selects_nothing() -> None: + """Falling back must not invent coverage where the geometry does not reach.""" + + selection = select_cells(box(300_000, 300_000, 300_100, 300_100), out_shape=SHAPE, transform=TRANSFORM) + + assert selection.mask.sum() == 0 + assert selection.expanded_to_touched_cells is False + assert selection.mode == "cell_centre" + + +def test_the_warning_states_how_much_larger_the_analysed_area_is() -> None: + selection = select_cells( + box(200_010, 209_960, 200_050, 209_990), + out_shape=SHAPE, + transform=TRANSFORM, + cell_area_m2=100.0 * 100.0, + ) + + # One 100x100 m cell was analysed for a 40x30 m request. + assert "1 rastercel" in selection.warning + assert "1.0 ha" in selection.warning + + +def test_the_mask_shape_always_matches_the_raster_window() -> None: + selection = select_cells(box(200_010, 209_960, 200_050, 209_990), out_shape=SHAPE, transform=TRANSFORM) + + assert selection.mask.shape == SHAPE + assert selection.mask.dtype == np.bool_ diff --git a/backend/tests/test_small_selection_raster_analysis.py b/backend/tests/test_small_selection_raster_analysis.py new file mode 100644 index 00000000..ee5bc1a1 --- /dev/null +++ b/backend/tests/test_small_selection_raster_analysis.py @@ -0,0 +1,117 @@ +"""A selection finer than the source raster must answer, not return zero. + +End-to-end counterpart to ``test_raster_cell_selection``: the analysis reads a +real GeoTIFF, so it proves the fallback survives the clip/mask path the service +actually uses rather than only the helper in isolation. +""" + +from __future__ import annotations + +from pathlib import Path +from uuid import uuid4 + +import pytest + +np = pytest.importorskip("numpy") +rasterio = pytest.importorskip("rasterio") + +from pyproj import Transformer +from rasterio.transform import from_origin + +from app.core.config import Settings +from app.models import Dataset +from app.schemas.flood_hazard import FloodHazardSelectionRequest +from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService +from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService + + +TO_4326 = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) +PRODUCT_KEY = "fluviaal_current_t100" + + +class FakeSession: + def __init__(self, objects): + self.objects = objects + + def get(self, model, item_id): + return self.objects.get((model, item_id)) + + +def _write_raster(path: Path, *, resolution: float, depth: float) -> None: + values = np.full((4, 4), depth, dtype="float32") + with rasterio.open( + path, + "w", + driver="GTiff", + width=4, + height=4, + count=1, + dtype="float32", + crs="EPSG:31370", + transform=from_origin(200_000, 210_000, resolution, resolution), + nodata=-9999.0, + ) as output: + output.write(values, 1) + + +def _dataset(project_id, dataset_id, path: Path) -> Dataset: + return Dataset( + id=dataset_id, + project_id=project_id, + name="vmm-flood.tif", + dataset_type="raster", + source="vmm", + source_name=FloodHazardAcquisitionService.PROVIDER, + status="ready", + storage_path=str(path), + source_metadata={"product_key": PRODUCT_KEY, "normalized_value_unit": "m"}, + ) + + +def _bbox_for(min_x: float, min_y: float, max_x: float, max_y: float) -> dict: + left, bottom = TO_4326.transform(min_x, min_y) + right, top = TO_4326.transform(max_x, max_y) + return {"min_x": left, "min_y": bottom, "max_x": right, "max_y": top, "crs": "EPSG:4326"} + + +def _analyze(tmp_path: Path, bbox: dict, *, resolution: float = 100.0) -> dict: + project_id = uuid4() + dataset_id = uuid4() + path = tmp_path / "flood.tif" + _write_raster(path, resolution=resolution, depth=2.0) + dataset = _dataset(project_id, dataset_id, path) + db = FakeSession({(Dataset, dataset_id): dataset}) + + return FloodHazardAnalysisService.analyze( + db, + project_id, + dataset_id, + FloodHazardSelectionRequest(bbox=bbox), + settings=Settings(_env_file=None), + ) + + +def test_a_selection_smaller_than_one_cell_reports_the_cell_it_touches(tmp_path: Path) -> None: + # A 40 x 30 m rectangle wholly inside one 100 m cell: no cell centre falls + # inside it, so the centre rule alone would report an empty selection. + result = _analyze(tmp_path, _bbox_for(200_010, 209_960, 200_050, 209_990)) + + assert result["inundated_cell_count"] == 1 + assert result["inundated_fraction"] == pytest.approx(1.0) + assert "kleiner dan één rastercel" in result["coverage_warning"] + + +def test_a_normal_selection_is_unaffected(tmp_path: Path) -> None: + result = _analyze(tmp_path, _bbox_for(200_000, 209_700, 200_300, 210_000)) + + assert result["inundated_cell_count"] >= 9 + assert result["coverage_warning"] is None + + +def test_the_reported_area_matches_the_cells_that_were_analysed(tmp_path: Path) -> None: + result = _analyze(tmp_path, _bbox_for(200_010, 209_960, 200_050, 209_990)) + metrics = {item["metric_key"]: item["metric_value"] for item in result["summary"]["metrics"]} + + # One 100 x 100 m cell, not the 0.12 ha that was drawn. + assert metrics["modelled_inundated_area_ha"] == pytest.approx(1.0) + assert metrics["selection_area_ha"] == pytest.approx(1.0) diff --git a/backend/tests/test_sprint187_temporal_map_foundation.py b/backend/tests/test_sprint187_temporal_map_foundation.py index 962473f6..3c25ad36 100644 --- a/backend/tests/test_sprint187_temporal_map_foundation.py +++ b/backend/tests/test_sprint187_temporal_map_foundation.py @@ -292,13 +292,13 @@ def test_population_area_weighting_is_exact_for_full_features_and_estimated_for_ bbox = {"min_x": 5.0, "min_y": 51.1, "max_x": 5.2, "max_y": 51.3, "crs": "EPSG:4326"} full = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([38_675.0, 0]), + SequenceScalarSession([49, 38_675.0]), dataset=dataset, bbox=bbox, total_feature_count=49, ) partial = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([1_250.5, 2]), + SequenceScalarSession([1, 1_250.5]), dataset=dataset, bbox=bbox, total_feature_count=3, diff --git a/backend/tests/test_sprint201_semantic_selection_metrics.py b/backend/tests/test_sprint201_semantic_selection_metrics.py index 2aa22054..3b9024f0 100644 --- a/backend/tests/test_sprint201_semantic_selection_metrics.py +++ b/backend/tests/test_sprint201_semantic_selection_metrics.py @@ -24,8 +24,15 @@ class ScalarQuery: class SequenceScalarSession: - def __init__(self, values: list[float]): - self.values = iter(values) + """Answers scalar queries in the order the summary issues them. + + The first query is the fully-covered feature count that produces the + selection-edge disclosure; ``covered_count`` defaults to the full + population, i.e. a selection that cuts nothing. + """ + + def __init__(self, values: list[float], covered_count: float | None = None): + self.values = iter(([covered_count] if covered_count is not None else []) + values) def query(self, *args): # noqa: ANN002, ARG002 return ScalarQuery(next(self.values)) @@ -53,7 +60,7 @@ def themed_dataset(theme: str, *, method: str = "feature_count") -> Dataset: def test_building_selection_promotes_footprint_area_and_retains_object_count() -> None: result = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([125_000.0]), + SequenceScalarSession([125_000.0], covered_count=40), dataset=themed_dataset("buildings"), bbox=BBOX, total_feature_count=40, @@ -73,7 +80,7 @@ def test_building_selection_promotes_footprint_area_and_retains_object_count() - def test_water_selection_reports_surface_length_and_honest_volume_limitation() -> None: result = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([52_500.0, 12_750.0]), + SequenceScalarSession([52_500.0, 12_750.0], covered_count=23), dataset=themed_dataset("water"), bbox=BBOX, total_feature_count=23, @@ -95,7 +102,7 @@ def test_population_keeps_configured_metric_and_adds_sector_count() -> None: {"metric_key": "population", "property": "population_total", "label": "Inwoners", "unit": "inwoners"} ) result = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([86_458.0]), + SequenceScalarSession([86_458.0], covered_count=733), dataset=dataset, bbox=BBOX, total_feature_count=733, @@ -132,7 +139,7 @@ def test_station_measurement_uses_numeric_mean_without_area_extrapolation() -> N ) result = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([30.455]), + SequenceScalarSession([30.455], covered_count=1), dataset=dataset, bbox=BBOX, total_feature_count=1, @@ -152,7 +159,7 @@ def test_regional_historical_polygons_do_not_emit_irrelevant_line_metrics() -> N dataset.provenance_metadata = {"operator_tool": "provision_regional_historical_landuse.py"} result = VectorFeatureService.summarize_features_by_bbox( - SequenceScalarSession([52_500.0]), + SequenceScalarSession([52_500.0], covered_count=23), dataset=dataset, bbox=BBOX, total_feature_count=23, diff --git a/backend/tests/test_sprint208_vmm_flood_hazard.py b/backend/tests/test_sprint208_vmm_flood_hazard.py index 44311b1e..ede20cf1 100644 --- a/backend/tests/test_sprint208_vmm_flood_hazard.py +++ b/backend/tests/test_sprint208_vmm_flood_hazard.py @@ -287,8 +287,20 @@ def test_flood_hazard_analysis_reports_scenario_metrics_without_claiming_waterbo metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]} assert result["inundated_cell_count"] == 200 - assert result["inundated_fraction"] == pytest.approx(0.5) + # The fixture models the left half and marks the right half nodata. All of + # the modelled half is wet, and the model covers half the selection. The + # earlier 0.5 conflated "not modelled" with "modelled dry" and reported + # half the risk that the model actually describes. + assert result["valid_cell_count"] == 200 + assert result["no_data_cell_count"] == 200 + assert result["inundated_fraction"] == pytest.approx(1.0) + assert result["data_coverage_ratio"] == pytest.approx(0.5) + assert "50.0%" in result["coverage_warning"] + assert metrics["modelled_inundated_share_pct"]["metric_value"] == pytest.approx(100.0) + assert metrics["model_coverage_pct"]["metric_value"] == pytest.approx(50.0) assert metrics["modelled_inundated_area_ha"]["metric_value"] == pytest.approx(0.5) + assert metrics["modelled_area_ha"]["metric_value"] == pytest.approx(0.5) + assert metrics["selection_area_ha"]["metric_value"] == pytest.approx(1.0) assert metrics["modelled_depth_mean_m"]["metric_value"] == pytest.approx(1.0) assert metrics["modelled_max_depth_area_integral_m3"]["metric_value"] == pytest.approx(5000.0) assert "concurrent_flood_volume_m3" in result["unsupported_metrics"] diff --git a/backend/tests/test_vector_selection_partial_features.py b/backend/tests/test_vector_selection_partial_features.py new file mode 100644 index 00000000..e2a8424f --- /dev/null +++ b/backend/tests/test_vector_selection_partial_features.py @@ -0,0 +1,161 @@ +"""The object count and the area metric must describe the same selection. + +``intersection_area`` clips a feature to the drawn rectangle, but the object +count treated any feature that merely touches the rectangle as wholly inside. +For a rectangle across a built-up area that overstates the count at every +edge, and the two headline numbers on the same panel then describe different +populations: "1.000 gebouwen" next to the clipped area of rather fewer. + +The count now reports how many features lie entirely inside and how many are +cut by the selection edge, and is marked as an estimate when any are. +""" + +from __future__ import annotations + +from app.services.vector_feature_service import VectorFeatureService + + +def test_a_count_without_partial_features_is_exact() -> None: + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=120, + fully_covered_feature_count=120, + ) + + assert disclosure["partially_covered_feature_count"] == 0 + assert disclosure["is_estimate"] is False + assert disclosure["warning"] is None + + +def test_features_cut_by_the_selection_edge_are_reported() -> None: + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=120, + fully_covered_feature_count=98, + ) + + assert disclosure["partially_covered_feature_count"] == 22 + assert disclosure["is_estimate"] is True + assert "22" in disclosure["warning"] + assert "rand" in disclosure["warning"] + + +def test_a_selection_of_only_partial_features_is_still_coherent() -> None: + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=3, + fully_covered_feature_count=0, + ) + + assert disclosure["partially_covered_feature_count"] == 3 + assert disclosure["is_estimate"] is True + + +def test_an_empty_selection_makes_no_claim() -> None: + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=0, + fully_covered_feature_count=0, + ) + + assert disclosure["partially_covered_feature_count"] == 0 + assert disclosure["is_estimate"] is False + assert disclosure["warning"] is None + + +def test_a_preclipped_full_area_selection_has_no_edge_effect() -> None: + """Selecting the whole work area cuts nothing; the count is exact.""" + + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=500, + fully_covered_feature_count=None, + ) + + assert disclosure["partially_covered_feature_count"] is None + assert disclosure["is_estimate"] is False + assert disclosure["warning"] is None + + +def test_an_inconsistent_covered_count_never_produces_a_negative() -> None: + disclosure = VectorFeatureService.count_disclosure( + total_feature_count=10, + fully_covered_feature_count=14, + ) + + assert disclosure["partially_covered_feature_count"] == 0 + assert disclosure["is_estimate"] is False + + +class _ScalarQuery: + def __init__(self, value): + self.value = value + + def filter(self, *args): # noqa: ANN002, ARG002 + return self + + def scalar(self): + return self.value + + +class _SequenceSession: + """Answers the summary's scalar queries in order: covered count, then metrics.""" + + def __init__(self, values): + self.values = iter(values) + + def query(self, *args): # noqa: ANN002, ARG002 + return _ScalarQuery(next(self.values)) + + +def _buildings_dataset(): + from uuid import uuid4 + + from app.models import Dataset + + return Dataset( + id=uuid4(), + project_id=uuid4(), + name="grb-buildings.geojson", + dataset_type="vector", + dataset_role="reference", + source_name="grb", + reference_layer_name="buildings", + source_metadata={"theme": "buildings"}, + ) + + +BBOX = {"min_x": 5.0, "min_y": 51.1, "max_x": 5.2, "max_y": 51.3, "crs": "EPSG:4326"} + + +def test_summary_reports_the_edge_cut_next_to_the_object_count() -> None: + summary = VectorFeatureService.summarize_features_by_bbox( + _SequenceSession([88, 125_000.0]), + dataset=_buildings_dataset(), + bbox=BBOX, + total_feature_count=100, + ) + + assert summary["feature_count"] == 100 + assert summary["fully_covered_feature_count"] == 88 + assert summary["partially_covered_feature_count"] == 12 + assert "12 van de 100" in summary["selection_edge_warning"] + + count_metric = next( + item for item in summary["metrics"] if item["aggregation_method"] == "feature_count" + ) + assert count_metric["is_estimate"] is True + assert "doorgesneden" in count_metric["warning"] + + # The clipped area metric is exact and must not inherit the count's caveat. + area_metric = next( + item for item in summary["metrics"] if item["aggregation_method"] == "intersection_area" + ) + assert area_metric["is_estimate"] is False + + +def test_summary_stays_exact_when_the_selection_cuts_nothing() -> None: + summary = VectorFeatureService.summarize_features_by_bbox( + _SequenceSession([100, 125_000.0]), + dataset=_buildings_dataset(), + bbox=BBOX, + total_feature_count=100, + ) + + assert summary["partially_covered_feature_count"] == 0 + assert summary["selection_edge_warning"] is None