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>
291 lines
16 KiB
Python
291 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.thematic_raster import (
|
|
ThematicRasterMetric,
|
|
ThematicRasterSelectionRequest,
|
|
ThematicRasterSelectionResponse,
|
|
ThematicRasterSelectionSummary,
|
|
)
|
|
from app.services.thematic_raster_acquisition_service import (
|
|
ThematicRasterAcquisitionService,
|
|
ThematicRasterProduct,
|
|
)
|
|
|
|
|
|
class ThematicRasterAnalysisService:
|
|
@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 != ThematicRasterAcquisitionService.PROVIDER:
|
|
raise AppError(
|
|
code="INVALID_THEMATIC_RASTER_DATASET",
|
|
message="Thematic analysis requires a governed Departement Omgeving 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 thematic raster file is unavailable", status_code=404)
|
|
return dataset
|
|
|
|
@staticmethod
|
|
def _product(dataset: Dataset) -> ThematicRasterProduct:
|
|
source_metadata = dataset.source_metadata or {}
|
|
product = ThematicRasterAcquisitionService._products().get(str(source_metadata.get("product_key") or ""))
|
|
if product is None or source_metadata.get("coverage_id") != product.coverage_id:
|
|
raise AppError(code="INVALID_THEMATIC_RASTER_METADATA", message="Thematic raster provenance is incomplete", status_code=409)
|
|
return product
|
|
|
|
@staticmethod
|
|
def _selection_geometry(db, project_id: UUID, payload: ThematicRasterSelectionRequest):
|
|
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="THEMATIC_RASTER_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422)
|
|
return selection
|
|
|
|
@staticmethod
|
|
def _unsupported_metrics(product: ThematicRasterProduct) -> list[str]:
|
|
if product.metric_kind == "binary_area":
|
|
if product.theme == "forest":
|
|
return ["tree_count", "canopy_cover", "timber_volume", "legal_forest_boundary"]
|
|
if product.theme == "agriculture":
|
|
return ["declared_parcel_area", "crop_declaration", "ownership", "cadastral_area"]
|
|
return ["object_count", "parcel_area", "current_land_use"]
|
|
if product.metric_kind == "population_density":
|
|
return ["current_population", "household_count", "address_level_population"]
|
|
if product.metric_kind == "index_score":
|
|
return ["travel_time_minutes", "current_timetable", "stop_count"]
|
|
return ["facility_count", "opening_hours", "current_service_availability"]
|
|
|
|
@staticmethod
|
|
def analyze(
|
|
db,
|
|
project_id: UUID,
|
|
dataset_id: UUID,
|
|
payload: ThematicRasterSelectionRequest,
|
|
*,
|
|
settings: Settings | None = None,
|
|
) -> dict:
|
|
resolved_settings = settings or get_settings()
|
|
dataset = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id)
|
|
product = ThematicRasterAnalysisService._product(dataset)
|
|
selection_4326 = ThematicRasterAnalysisService._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 thematic raster 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() != 31370:
|
|
raise AppError(code="INVALID_DATASET_CRS", message="Thematic raster CRS must be EPSG:31370", 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="THEMATIC_RASTER_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted thematic 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.thematic_raster_max_pixels:
|
|
raise AppError(
|
|
code="THEMATIC_RASTER_SELECTION_TOO_LARGE",
|
|
message="Thematic raster analysis exceeds the configured cell limit",
|
|
details={"pixel_count": expected_cells, "max_pixels": resolved_settings.thematic_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 = 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))
|
|
values = raw[valid]
|
|
ThematicRasterAcquisitionService._validate_values(values, product)
|
|
selected_cell_count = int(selected.sum())
|
|
valid_cell_count = int(values.size)
|
|
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="THEMATIC_RASTER_ANALYSIS_FAILED",
|
|
message="The persisted thematic 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, *, estimate: bool = True) -> ThematicRasterMetric:
|
|
return ThematicRasterMetric(
|
|
metric_key=key,
|
|
metric_label=label,
|
|
metric_value=round(float(value), 4),
|
|
metric_unit=unit,
|
|
aggregation_method=method,
|
|
is_estimate=estimate,
|
|
)
|
|
|
|
if product.metric_kind == "binary_area":
|
|
positive_count = int(np.count_nonzero(values >= 0.5))
|
|
positive_area_ha = positive_count * cell_area_m2 / 10_000.0
|
|
positive_share = positive_count / max(1, valid_cell_count) * 100.0
|
|
label = {
|
|
"space_occupation": "Ruimtebeslag",
|
|
"open_space": "Open ruimte",
|
|
"forest": "Bos",
|
|
"agriculture": "Akker en landbouwgrasland",
|
|
}[product.theme]
|
|
metrics = [
|
|
metric(f"{product.theme}_area_ha", f"{label} in selectie", positive_area_ha, "ha", "positive_source_cells_times_cell_area"),
|
|
metric(f"{product.theme}_share_pct", f"Aandeel {label.lower()}", positive_share, "%", "positive_source_cells_divided_by_valid_selected_cells"),
|
|
metric("valid_raster_area_ha", "Rasteroppervlakte met bronwaarde", valid_cell_count * cell_area_m2 / 10_000.0, "ha", "valid_selected_cells_times_cell_area"),
|
|
]
|
|
elif product.metric_kind == "population_density":
|
|
estimated_population = float(values.sum() * (cell_area_m2 / 10_000.0))
|
|
metrics = [
|
|
metric("estimated_inhabitants", "Geraamd aantal inwoners (2019)", estimated_population, "inwoners", "sum_density_times_selected_cell_area_hectares"),
|
|
metric("population_density_mean_per_ha", "Gemiddelde inwonersdichtheid", values.mean(), "inwoners/ha", "mean_valid_one_hectare_source_cells"),
|
|
metric("population_density_p90_per_ha", "90e percentiel inwonersdichtheid", np.percentile(values, 90), "inwoners/ha", "percentile_90_valid_source_cells"),
|
|
]
|
|
else:
|
|
unit = "score" if product.metric_kind == "index_score" else "score (0-1)"
|
|
label = "Knooppuntwaarde" if product.metric_kind == "index_score" else "Voorzieningenniveau"
|
|
metrics = [
|
|
metric(f"{product.theme}_mean", f"Gemiddelde {label.lower()}", values.mean(), unit, "mean_valid_source_cells"),
|
|
metric(f"{product.theme}_p10", f"10e percentiel {label.lower()}", np.percentile(values, 10), unit, "percentile_10_valid_source_cells"),
|
|
metric(f"{product.theme}_median", f"Mediaan {label.lower()}", np.percentile(values, 50), unit, "median_valid_source_cells"),
|
|
metric(f"{product.theme}_p90", f"90e percentiel {label.lower()}", np.percentile(values, 90), unit, "percentile_90_valid_source_cells"),
|
|
]
|
|
|
|
primary = metrics[0]
|
|
response = ThematicRasterSelectionResponse(
|
|
dataset_id=dataset.id,
|
|
product_key=product.key,
|
|
theme=product.theme,
|
|
metric_kind=product.metric_kind,
|
|
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(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(
|
|
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=ThematicRasterAnalysisService._unsupported_metrics(product),
|
|
limitation_message=product.limitation_message,
|
|
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 = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id)
|
|
product = ThematicRasterAnalysisService._product(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 thematic raster rendering", status_code=503) from exc
|
|
palettes = {
|
|
"space_occupation": np.asarray([[251, 231, 211], [190, 62, 51]], dtype="float64"),
|
|
"open_space": np.asarray([[221, 238, 219], [38, 122, 70]], dtype="float64"),
|
|
"forest": np.asarray([[223, 237, 226], [43, 117, 72]], dtype="float64"),
|
|
"agriculture": np.asarray([[245, 237, 204], [166, 122, 35]], dtype="float64"),
|
|
"population": np.asarray([[238, 231, 246], [103, 58, 151]], dtype="float64"),
|
|
"accessibility": np.asarray([[233, 241, 244], [15, 118, 110]], dtype="float64"),
|
|
"services": np.asarray([[255, 244, 191], [182, 109, 22]], dtype="float64"),
|
|
}
|
|
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))
|
|
resampling = Resampling.nearest if product.metric_kind == "binary_area" else Resampling.bilinear
|
|
data = source.read(1, out_shape=(height, width), masked=True, resampling=resampling)
|
|
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 product.metric_kind == "binary_area":
|
|
valid &= values >= 0.5
|
|
normalized = np.where(valid, 1.0, 0.0)
|
|
else:
|
|
source_metadata = dataset.source_metadata or {}
|
|
lower = float(source_metadata.get("render_min_value", np.nanpercentile(values[valid], 2) if valid.any() else 0.0))
|
|
upper = float(source_metadata.get("render_max_value", np.nanpercentile(values[valid], 98) if valid.any() else 1.0))
|
|
if upper <= lower:
|
|
upper = lower + 1.0
|
|
normalized = np.clip((values - lower) / (upper - lower), 0.0, 1.0)
|
|
colors = palettes[product.theme]
|
|
rgba = np.zeros((height, width, 4), dtype="uint8")
|
|
for channel in range(3):
|
|
rgba[:, :, channel] = (colors[0, channel] + normalized * (colors[1, channel] - colors[0, channel])).astype("uint8")
|
|
rgba[:, :, 3] = np.where(valid, 205, 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="THEMATIC_RASTER_PREVIEW_FAILED",
|
|
message="The persisted thematic raster could not be rendered",
|
|
details={"reason": str(exc)},
|
|
status_code=500,
|
|
) from exc
|