diff --git a/.env.example b/.env.example index 99044329..c80c35e8 100644 --- a/.env.example +++ b/.env.example @@ -12,6 +12,14 @@ ORTHOPHOTO_RESOLUTION_M=1.0 ORTHOPHOTO_MIN_SIDE_M=128 ORTHOPHOTO_MAX_SIDE_M=1024 ORTHOPHOTO_CACHE_TTL_HOURS=24 +DHMV_ENABLED=true +DHMV_WCS_URL=https://geo.api.vlaanderen.be/DHMV/wcs +DHMV_RESOLUTION_M=5.0 +DHMV_MIN_SIDE_M=10 +DHMV_MAX_SIDE_M=20000 +DHMV_MAX_PIXELS=12000000 +DHMV_TIMEOUT_SECONDS=300 +DHMV_MAX_RESPONSE_MB=160 YOLO_ENABLED=false YOLO_MODELS_DIR=/app/models YOLO_MODEL_PATH= diff --git a/CHANGELOG.md b/CHANGELOG.md index be5c3d58..45012190 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,22 @@ # Changelog +## Sprint 207 Governed DHMV II terrain foundation (2026-07-15) + +- Added a fixed official Digitaal Vlaanderen DHMV II DTM/DSM WCS registry, + bounded multipart GeoTIFF acquisition and exact persisted-Area clipping. +- Validates and retains native 1 m product identity, 5 m analysis resolution, + EPSG:31370, Float32 nodata, TAW, acquisition period and source/output + checksums through ordinary Dataset, DatasetVersion and Job persistence. +- Added exact raster-selection metrics for mean/min/max/P10/P90 height, relief + and slope, with explicit valid-cell coverage and computation methods. +- Added a Mol operator and a `Hoogte & reliëf` map theme with colour-relief + MapLibre overlay over the existing OpenStreetMap context. +- Explicitly keeps drainage as a future derived analysis and prohibits + presenting DHMV terrain/surface height as water depth or volume. +- Added runtime/Unraid configuration and focused acquisition, GIS formula, + persistence, API, frontend and packaging tests without a migration. + ## Sprint 206 Governed Buildings and Addresses Register snapshot (2026-07-15) - Added an explicit operator for the current official Digitaal Vlaanderen diff --git a/backend/README.md b/backend/README.md index 77923b4a..4b68db8c 100644 --- a/backend/README.md +++ b/backend/README.md @@ -1184,6 +1184,32 @@ Settings: `ORTHOPHOTO_ENABLED`, `ORTHOPHOTO_WMS_URL`, `ORTHOPHOTO_CACHE_TTL_HOURS`. Keep the official HTTPS URL and 1 m profile unless a separately verified deployment/model profile requires a change. +## Governed DHMV terrain acquisition + +`GET /api/v1/projects/{project_id}/datasets/dhmv/products` exposes the fixed +official DTM/DSM registry. `POST .../datasets/dhmv/acquire` requests only +`DHMVII_DTM_1m` or `DHMVII_DSM_1m` from the production Digitaal Vlaanderen WCS. +The default 5 m analysis copy keeps complete-Mol processing bounded while +retaining native 1 m resolution, EPSG:31370, TAW, `-9999` nodata and the +2013-2015 acquisition period in provenance. + +Run the complete Mol operator after the regional workspace and Mol Area exist: + +```bash +docker exec geointel python /app/scripts/provision_mol_dhmv.py +``` + +The operator acquires DTM and DSM, clips each raster to the exact persisted +Area, validates checksums and calls the terrain selection endpoint as a smoke. +Use `--products dtm_1m`, `--resolution-m 5` or `--force` when explicitly +needed. `POST .../raster/terrain/select` returns height in m TAW, relief in +metres and slope in degrees. `GET .../raster/terrain/image` returns the +constrained MapLibre PNG. Water depth, volume and drainage remain unavailable. + +Settings: `DHMV_ENABLED`, `DHMV_WCS_URL`, `DHMV_RESOLUTION_M`, +`DHMV_MIN_SIDE_M`, `DHMV_MAX_SIDE_M`, `DHMV_MAX_PIXELS`, +`DHMV_TIMEOUT_SECONDS` and `DHMV_MAX_RESPONSE_MB`. + ## Waterinfo station histories Run the explicit operator after the regional workspace and Mol Area exist: diff --git a/backend/app/api/routes/datasets.py b/backend/app/api/routes/datasets.py index 325738bf..7a3faf16 100644 --- a/backend/app/api/routes/datasets.py +++ b/backend/app/api/routes/datasets.py @@ -22,6 +22,8 @@ from app.schemas import ( RasterNdwiRequest, RasterNdbiRequest, OrthophotoAcquireRequest, + DhmvAcquireRequest, + TerrainSelectionRequest, VectorBBoxResponse, VectorBufferRequest, VectorClipRequest, @@ -40,6 +42,8 @@ from app.services.vector_operations_service import VectorOperationsService from app.services.vector_feature_service import VectorFeatureService from app.services.dataset_service import DatasetService from app.services.orthophoto_acquisition_service import OrthophotoAcquisitionService +from app.services.dhmv_acquisition_service import DhmvAcquisitionService +from app.services.terrain_analysis_service import TerrainAnalysisService from app.utils.response import envelope router = APIRouter(prefix="/projects/{project_id}", tags=["datasets"]) @@ -151,6 +155,30 @@ def list_orthophoto_products(project_id: UUID, db: Session = Depends(get_db)): return envelope({"items": items, "total": len(items)}) +@router.post("/datasets/dhmv/acquire", response_model=dict) +def acquire_bounded_dhmv( + project_id: UUID, + payload: DhmvAcquireRequest, + db: Session = Depends(get_db), +): + job = JobService.run_sync_job( + db=db, + project_id=project_id, + job_type="raster.dhmv.acquire", + parameters=payload.model_dump(mode="json"), + operation=lambda: DhmvAcquisitionService.acquire(db, project_id, payload), + ) + return envelope(job) + + +@router.get("/datasets/dhmv/products", response_model=dict) +def list_dhmv_products(project_id: UUID, db: Session = Depends(get_db)): + if not db.get(Project, project_id): + raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404) + items = DhmvAcquisitionService.list_products() + return envelope({"items": items, "total": len(items)}) + + @router.get("/datasets", response_model=dict) def list_datasets( project_id: UUID, @@ -448,6 +476,30 @@ def raster_orthophoto_image( ) +@router.post("/datasets/{dataset_id}/raster/terrain/select", response_model=dict) +def raster_terrain_selection( + project_id: UUID, + dataset_id: UUID, + payload: TerrainSelectionRequest, + db: Session = Depends(get_db), +): + return envelope(TerrainAnalysisService.analyze(db, project_id, dataset_id, payload)) + + +@router.get("/datasets/{dataset_id}/raster/terrain/image") +def raster_terrain_image( + project_id: UUID, + dataset_id: UUID, + db: Session = Depends(get_db), +): + content = TerrainAnalysisService.render_png(db, project_id, dataset_id) + return Response( + content=content, + media_type="image/png", + headers={"Cache-Control": "private, max-age=86400"}, + ) + + @router.get("/datasets/{dataset_id}/raster/stats", response_model=dict) def raster_stats( project_id: UUID, diff --git a/backend/app/core/config.py b/backend/app/core/config.py index bbf0c9ab..0d0dfd9e 100644 --- a/backend/app/core/config.py +++ b/backend/app/core/config.py @@ -31,6 +31,17 @@ class Settings(BaseSettings): orthophoto_timeout_seconds: int = Field(default=120, ge=1, validation_alias="ORTHOPHOTO_TIMEOUT_SECONDS") orthophoto_max_response_mb: int = Field(default=32, ge=1, validation_alias="ORTHOPHOTO_MAX_RESPONSE_MB") orthophoto_cache_ttl_hours: int = Field(default=24, ge=0, validation_alias="ORTHOPHOTO_CACHE_TTL_HOURS") + dhmv_enabled: bool = Field(default=True, validation_alias="DHMV_ENABLED") + dhmv_wcs_url: str = Field( + default="https://geo.api.vlaanderen.be/DHMV/wcs", + validation_alias="DHMV_WCS_URL", + ) + dhmv_resolution_m: float = Field(default=5.0, ge=1.0, le=10.0, validation_alias="DHMV_RESOLUTION_M") + dhmv_min_side_m: float = Field(default=10.0, gt=0, validation_alias="DHMV_MIN_SIDE_M") + dhmv_max_side_m: float = Field(default=20_000.0, gt=0, validation_alias="DHMV_MAX_SIDE_M") + dhmv_max_pixels: int = Field(default=12_000_000, ge=1, validation_alias="DHMV_MAX_PIXELS") + dhmv_timeout_seconds: int = Field(default=300, ge=1, validation_alias="DHMV_TIMEOUT_SECONDS") + dhmv_max_response_mb: int = Field(default=160, ge=1, validation_alias="DHMV_MAX_RESPONSE_MB") redis_url: str | None = Field(default=None, validation_alias="REDIS_URL") log_level: str = Field(default="INFO", validation_alias="GEOINTEL_LOG_LEVEL") database_statement_timeout_ms: int = Field(default=5_000, validation_alias="DATABASE_STATEMENT_TIMEOUT_MS") diff --git a/backend/app/schemas/__init__.py b/backend/app/schemas/__init__.py index d56c90d2..49f200b4 100644 --- a/backend/app/schemas/__init__.py +++ b/backend/app/schemas/__init__.py @@ -33,6 +33,15 @@ from .segmentation import ( from .health import HealthResponse, SystemCapabilities from .job import JobCreate, JobList, JobRead, JobStatus from .orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult, OrthophotoProductRead +from .dhmv import ( + DhmvAcquireRequest, + DhmvAcquisitionResult, + DhmvProductRead, + TerrainMetric, + TerrainSelectionRequest, + TerrainSelectionResponse, + TerrainSelectionSummary, +) from .external import ( ExternalFetchRequest, ExternalFetchResponse, @@ -135,6 +144,13 @@ __all__ = [ "OrthophotoAcquireRequest", "OrthophotoAcquisitionResult", "OrthophotoProductRead", + "DhmvAcquireRequest", + "DhmvAcquisitionResult", + "DhmvProductRead", + "TerrainMetric", + "TerrainSelectionRequest", + "TerrainSelectionResponse", + "TerrainSelectionSummary", "VectorBBoxResponse", "VectorClipRequest", "VectorBufferRequest", diff --git a/backend/app/schemas/dhmv.py b/backend/app/schemas/dhmv.py new file mode 100644 index 00000000..fb96f35b --- /dev/null +++ b/backend/app/schemas/dhmv.py @@ -0,0 +1,91 @@ +from __future__ import annotations + +from uuid import UUID + +from pydantic import BaseModel, Field + +from .operations import VectorSelectionBBox + + +class DhmvAcquireRequest(BaseModel): + bbox: VectorSelectionBBox + area_id: UUID | None = None + product_key: str = "dtm_1m" + resolution_m: float | None = Field(default=None, ge=1.0, le=10.0) + force_refresh: bool = False + + +class DhmvProductRead(BaseModel): + key: str + display_name: str + surface_model: str + coverage_id: str + native_resolution_m: float + source_crs: str + vertical_reference: str + acquisition_period: str + catalog_url: str + attribution: str + limitation_message: str + + +class DhmvAcquisitionResult(BaseModel): + output_dataset_id: UUID + reused: bool + provider: str + product_key: str + display_name: str + surface_model: str + coverage_id: str + native_resolution_m: float + resolution_m: float + width: int + height: int + valid_pixel_count: int + nodata_value: float + bbox_epsg4326: list[float] + bbox_epsg31370: list[float] + vertical_reference: str + acquisition_period: str + attribution: str + limitation_message: str + + +class TerrainSelectionRequest(BaseModel): + bbox: VectorSelectionBBox + area_id: UUID | None = None + + +class TerrainMetric(BaseModel): + metric_key: str + metric_label: str + metric_value: float + metric_unit: str + aggregation_method: str + derived: bool = True + + +class TerrainSelectionSummary(BaseModel): + metric_label: str + metric_value: float + metric_unit: str + aggregation_method: str + primary_metric_key: str + metrics: list[TerrainMetric] + + +class TerrainSelectionResponse(BaseModel): + dataset_id: UUID + product_key: str + surface_model: str + selection_bbox: VectorSelectionBBox + selection_area_id: UUID | None = None + sample_count: int + slope_sample_count: int + coverage_ratio: float + resolution_m: float + vertical_reference: str + summary: TerrainSelectionSummary + unsupported_metrics: list[str] + limitation_message: str + generated_at: str diff --git a/backend/app/services/dhmv_acquisition_service.py b/backend/app/services/dhmv_acquisition_service.py new file mode 100644 index 00000000..b0c4128e --- /dev/null +++ b/backend/app/services/dhmv_acquisition_service.py @@ -0,0 +1,501 @@ +from __future__ import annotations + +import hashlib +import json +import math +from dataclasses import dataclass +from datetime import UTC, datetime +from email.parser import BytesParser +from email.policy import default +from pathlib import Path +from typing import Any, Callable +from urllib.error import HTTPError, URLError +from urllib.parse import urlencode +from urllib.request import Request, urlopen +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, Project +from app.schemas.dhmv import DhmvAcquireRequest, DhmvAcquisitionResult, DhmvProductRead +from app.services.dataset_service import DatasetService + + +@dataclass(frozen=True) +class DhmvProduct: + key: str + display_name: str + surface_model: str + coverage_id: str + native_resolution_m: float + catalog_url: str + limitation_message: str + + +class DhmvAcquisitionService: + PROVIDER = "digitaal_vlaanderen_dhmv" + SOURCE_CRS = "EPSG:31370" + VERTICAL_REFERENCE = "TAW (Tweede Algemene Waterpassing)" + ACQUISITION_PERIOD = "2013-2015" + SOURCE_VERSION = "DHMV II 2014.01" + NODATA = -9999.0 + ATTRIBUTION = "Bron: Digitaal Vlaanderen, Digitaal Hoogtemodel Vlaanderen II" + LICENSE_NOTE = "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen." + DTM_CATALOG_URL = ( + "https://www.vlaanderen.be/datavindplaats/catalogus/" + "digitaal-hoogtemodel-vlaanderen-ii-dtm-raster-1-m" + ) + DSM_CATALOG_URL = ( + "https://www.vlaanderen.be/datavindplaats/catalogus/" + "digitaal-hoogtemodel-vlaanderen-ii-dsm-raster-1-m" + ) + + @staticmethod + def _products() -> dict[str, DhmvProduct]: + products = ( + DhmvProduct( + key="dtm_1m", + display_name="DHMV II terreinmodel (DTM)", + surface_model="terrain", + coverage_id="DHMVII_DTM_1m", + native_resolution_m=1.0, + catalog_url=DhmvAcquisitionService.DTM_CATALOG_URL, + limitation_message=( + "Maaiveldhoogte uit de opnameperiode 2013-2015. Gebouwen en andere objecten zijn verwijderd. " + "Afstroming is een afgeleide interpretatie; dit product bevat geen waterdiepte." + ), + ), + DhmvProduct( + key="dsm_1m", + display_name="DHMV II oppervlaktemodel (DSM)", + surface_model="surface", + coverage_id="DHMVII_DSM_1m", + native_resolution_m=1.0, + catalog_url=DhmvAcquisitionService.DSM_CATALOG_URL, + limitation_message=( + "Oppervlaktehoogte uit de opnameperiode 2013-2015, inclusief gebouwen en vegetatie. " + "Dit is geen maaiveldmodel, waterdiepte of rechtstreeks gebouwhoogteproduct." + ), + ), + ) + return {product.key: product for product in products} + + @staticmethod + def list_products() -> list[dict[str, Any]]: + return [ + DhmvProductRead( + key=product.key, + display_name=product.display_name, + surface_model=product.surface_model, + coverage_id=product.coverage_id, + native_resolution_m=product.native_resolution_m, + source_crs=DhmvAcquisitionService.SOURCE_CRS, + vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE, + acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD, + catalog_url=product.catalog_url, + attribution=DhmvAcquisitionService.ATTRIBUTION, + limitation_message=product.limitation_message, + ).model_dump() + for product in DhmvAcquisitionService._products().values() + ] + + @staticmethod + def _product(product_key: str) -> DhmvProduct: + product = DhmvAcquisitionService._products().get(product_key.strip().lower()) + if product is None: + raise AppError( + code="DHMV_PRODUCT_NOT_SUPPORTED", + message="Select DTM or DSM from the governed DHMV II product registry", + details={"product_key": product_key}, + status_code=422, + ) + return product + + @staticmethod + def _prepared_request(payload: DhmvAcquireRequest, settings: Settings) -> dict[str, Any]: + if not settings.dhmv_enabled: + raise AppError(code="DHMV_NOT_CONFIGURED", message="DHMV acquisition is disabled", status_code=503) + product = DhmvAcquisitionService._product(payload.product_key) + resolution_m = float(payload.resolution_m or settings.dhmv_resolution_m) + if resolution_m < product.native_resolution_m or resolution_m > 10.0: + raise AppError( + code="DHMV_RESOLUTION_NOT_SUPPORTED", + message="DHMV analysis resolution must be between the native 1 metre and 10 metres", + status_code=422, + ) + values = (payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y) + if not all(math.isfinite(value) for value in values) or payload.bbox.min_x >= payload.bbox.max_x or payload.bbox.min_y >= payload.bbox.max_y: + raise AppError(code="INVALID_BBOX", message="DHMV selection must be a finite non-empty rectangle", status_code=400) + + transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.SOURCE_CRS, always_xy=True) + lambert_bounds = transformer.transform_bounds(*values, densify_pts=21) + width_m = float(lambert_bounds[2] - lambert_bounds[0]) + height_m = float(lambert_bounds[3] - lambert_bounds[1]) + if width_m < settings.dhmv_min_side_m or height_m < settings.dhmv_min_side_m: + raise AppError( + code="DHMV_SELECTION_TOO_SMALL", + message=f"Select an area of at least {settings.dhmv_min_side_m:g} by {settings.dhmv_min_side_m:g} metres", + status_code=422, + ) + if width_m > settings.dhmv_max_side_m or height_m > settings.dhmv_max_side_m: + raise AppError( + code="DHMV_SELECTION_TOO_LARGE", + message=f"Select an area no larger than {settings.dhmv_max_side_m:g} by {settings.dhmv_max_side_m:g} metres", + details={"width_m": width_m, "height_m": height_m}, + status_code=422, + ) + width = max(1, math.ceil(width_m / resolution_m)) + height = max(1, math.ceil(height_m / resolution_m)) + if width * height > settings.dhmv_max_pixels: + raise AppError( + code="DHMV_SELECTION_TOO_LARGE", + message="DHMV selection exceeds the configured raster cell limit", + details={"pixel_count": width * height, "max_pixels": settings.dhmv_max_pixels}, + status_code=422, + ) + + bbox_4326 = [float(value) for value in values] + bbox_31370 = [float(value) for value in lambert_bounds] + request_identity = { + "provider": DhmvAcquisitionService.PROVIDER, + "coverage_id": product.coverage_id, + "bbox_epsg4326": [round(value, 8) for value in bbox_4326], + "bbox_epsg31370": [round(value, 3) for value in bbox_31370], + "resolution_m": resolution_m, + "area_id": str(payload.area_id) if payload.area_id else None, + } + request_hash = hashlib.sha256(json.dumps(request_identity, sort_keys=True).encode()).hexdigest() + params = { + "SERVICE": "WCS", + "VERSION": "2.0.1", + "REQUEST": "GetCoverage", + "COVERAGEID": product.coverage_id, + "FORMAT": "image/tiff", + "SUBSET": [ + f"x({bbox_31370[0]:.3f},{bbox_31370[2]:.3f})", + f"y({bbox_31370[1]:.3f},{bbox_31370[3]:.3f})", + ], + "SCALEFACTOR": f"{resolution_m / product.native_resolution_m:g}", + } + query = [ + ("SERVICE", params["SERVICE"]), + ("VERSION", params["VERSION"]), + ("REQUEST", params["REQUEST"]), + ("COVERAGEID", params["COVERAGEID"]), + ("FORMAT", params["FORMAT"]), + ("SUBSET", params["SUBSET"][0]), + ("SUBSET", params["SUBSET"][1]), + ("SCALEFACTOR", params["SCALEFACTOR"]), + ] + return { + **request_identity, + "product": product, + "request_hash": request_hash, + "request_url": f"{settings.dhmv_wcs_url}?{urlencode(query)}", + "params": params, + "bbox_epsg4326": bbox_4326, + "bbox_epsg31370": bbox_31370, + "width": width, + "height": height, + } + + @staticmethod + def _scope_geometry(db, project_id: UUID, area_id: UUID | None, bbox_epsg4326: list[float]): + if not db.get(Project, project_id): + raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404) + selection = box(*bbox_epsg4326) + if area_id is None: + return selection + area = db.get(Area, 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) + intersection = to_shape(area.geometry).intersection(selection) + if intersection.is_empty or intersection.area <= 0: + raise AppError( + code="DHMV_SELECTION_OUTSIDE_AREA", + message="The DHMV selection does not overlap the selected work area", + status_code=422, + ) + return intersection + + @staticmethod + def _cached_dataset(db, project_id: UUID, filename: str) -> Dataset | None: + candidate = ( + db.query(Dataset) + .filter( + Dataset.project_id == project_id, + Dataset.name == filename, + Dataset.source_name == DhmvAcquisitionService.PROVIDER, + Dataset.status == "ready", + ) + .order_by(Dataset.imported_at.desc()) + .first() + ) + if candidate and candidate.storage_path and Path(candidate.storage_path).is_file(): + return candidate + return None + + @staticmethod + def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]: + request = Request(request_url, headers={"User-Agent": "GeoIntel/0.1 bounded-dhmv-acquisition"}) + max_bytes = settings.dhmv_max_response_mb * 1024 * 1024 + try: + with (opener or urlopen)(request, timeout=settings.dhmv_timeout_seconds) as response: + content_type = str(response.headers.get("Content-Type", "")) + content_length = response.headers.get("Content-Length") + if content_length and int(content_length) > max_bytes: + raise AppError(code="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502) + content = response.read(max_bytes + 1) + except AppError: + raise + except (HTTPError, URLError, TimeoutError, OSError) as exc: + raise AppError( + code="DHMV_PROVIDER_UNAVAILABLE", + message="The official DHMV WCS could not complete the bounded request", + details={"reason": str(exc)}, + status_code=502, + ) from exc + if len(content) > max_bytes: + raise AppError(code="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502) + return content, content_type + + @staticmethod + def _extract_geotiff(content: bytes, content_type: str) -> bytes: + if content.startswith((b"II*\x00", b"MM\x00*")): + return content + if "multipart" not in content_type.lower(): + preview = content[:300].decode("utf-8", errors="replace") + raise AppError( + code="DHMV_PROVIDER_INVALID_RESPONSE", + message="The official DHMV service did not return a GeoTIFF coverage", + details={"content_type": content_type, "response_preview": preview}, + status_code=502, + ) + message = BytesParser(policy=default).parsebytes( + f"Content-Type: {content_type}\r\nMIME-Version: 1.0\r\n\r\n".encode() + content + ) + for part in message.iter_parts(): + payload = part.get_payload(decode=True) or b"" + if part.get_content_type() == "image/tiff" and payload.startswith((b"II*\x00", b"MM\x00*")): + return payload + raise AppError( + code="DHMV_PROVIDER_INVALID_RESPONSE", + message="The official DHMV multipart response contains no valid GeoTIFF coverage", + status_code=502, + ) + + @staticmethod + def _normalize_raster(content: bytes, scope_geometry_4326, prepared: dict[str, Any]) -> tuple[bytes, dict[str, Any]]: + try: + import numpy as np + import rasterio + from rasterio.io import MemoryFile + from rasterio.mask import mask + except ImportError as exc: + raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for DHMV validation", status_code=503) from exc + + try: + with MemoryFile(content) as source_memory, source_memory.open() as source: + if source.crs is None or source.crs.to_epsg() != 31370: + raise AppError(code="DHMV_INVALID_CRS", message="DHMV coverage must use EPSG:31370", status_code=502) + if source.count != 1: + raise AppError(code="DHMV_INVALID_BANDS", message="DHMV coverage must contain exactly one elevation band", status_code=502) + resolution = max(abs(float(source.res[0])), abs(float(source.res[1]))) + if not math.isclose(resolution, prepared["resolution_m"], rel_tol=0.02, abs_tol=0.05): + raise AppError( + code="DHMV_INVALID_RESOLUTION", + message="DHMV coverage resolution differs from the governed request", + details={"expected_m": prepared["resolution_m"], "actual_m": resolution}, + status_code=502, + ) + transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.SOURCE_CRS, always_xy=True) + scope_metric = shapely_transform(transformer.transform, scope_geometry_4326) + clipped, transform = mask( + source, + [mapping(scope_metric)], + crop=True, + filled=False, + indexes=[1], + ) + band = np.ma.asarray(clipped[0], dtype="float32") + nodata = float(source.nodata if source.nodata is not None else DhmvAcquisitionService.NODATA) + invalid = ~np.isfinite(np.asarray(band.filled(np.nan), dtype="float64")) + combined_mask = np.ma.getmaskarray(band) | invalid | (np.asarray(band) == nodata) + normalized = np.ma.array(np.asarray(band, dtype="float32"), mask=combined_mask) + valid_pixel_count = int(normalized.count()) + if valid_pixel_count == 0: + raise AppError(code="DHMV_NO_VALID_DATA", message="DHMV coverage contains no valid elevation cells in this selection", status_code=422) + profile = source.profile.copy() + profile.pop("blockxsize", None) + profile.pop("blockysize", None) + profile.update( + driver="GTiff", + width=int(normalized.shape[1]), + height=int(normalized.shape[0]), + count=1, + dtype="float32", + crs=DhmvAcquisitionService.SOURCE_CRS, + transform=transform, + nodata=DhmvAcquisitionService.NODATA, + compress="deflate", + predictor=3, + ) + with MemoryFile() as output_memory: + with output_memory.open(**profile) as output: + output.write(normalized.filled(DhmvAcquisitionService.NODATA), 1) + normalized_content = output_memory.read() + valid_values = normalized.compressed().astype("float64") + return normalized_content, { + "width": int(normalized.shape[1]), + "height": int(normalized.shape[0]), + "valid_pixel_count": valid_pixel_count, + "nodata_value": DhmvAcquisitionService.NODATA, + "resolution_m": resolution, + "minimum_m_taw": float(valid_values.min()), + "maximum_m_taw": float(valid_values.max()), + } + except AppError: + raise + except Exception as exc: + raise AppError( + code="DHMV_RASTER_INVALID", + message="The official DHMV response could not be validated as a georeferenced elevation raster", + details={"reason": str(exc)}, + status_code=502, + ) from exc + + @staticmethod + def acquire( + db, + project_id: UUID, + payload: DhmvAcquireRequest, + *, + settings: Settings | None = None, + opener: Callable[..., Any] | None = None, + ) -> dict[str, Any]: + resolved_settings = settings or get_settings() + prepared = DhmvAcquisitionService._prepared_request(payload, resolved_settings) + product: DhmvProduct = prepared["product"] + scope_geometry = DhmvAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"]) + resolution_token = f"{prepared['resolution_m']:g}".replace(".", "p") + filename = f"dhmvii_{product.surface_model}_{resolution_token}m_{prepared['request_hash'][:12]}.tif" + if not payload.force_refresh: + cached = DhmvAcquisitionService._cached_dataset(db, project_id, filename) + if cached is not None: + source_metadata = cached.source_metadata or {} + raster_metadata = cached.metadata_json or {} + return DhmvAcquisitionResult( + output_dataset_id=cached.id, + reused=True, + provider=DhmvAcquisitionService.PROVIDER, + product_key=product.key, + display_name=product.display_name, + surface_model=product.surface_model, + coverage_id=product.coverage_id, + native_resolution_m=product.native_resolution_m, + resolution_m=float(source_metadata.get("analysis_resolution_m", prepared["resolution_m"])), + width=int(raster_metadata.get("width", prepared["width"])), + height=int(raster_metadata.get("height", prepared["height"])), + valid_pixel_count=int(source_metadata.get("valid_pixel_count", 0)), + nodata_value=float(raster_metadata.get("nodata", DhmvAcquisitionService.NODATA)), + bbox_epsg4326=prepared["bbox_epsg4326"], + bbox_epsg31370=prepared["bbox_epsg31370"], + vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE, + acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD, + attribution=DhmvAcquisitionService.ATTRIBUTION, + limitation_message=product.limitation_message, + ).model_dump(mode="json") + + raw_content, content_type = DhmvAcquisitionService._fetch(prepared["request_url"], resolved_settings, opener) + coverage_content = DhmvAcquisitionService._extract_geotiff(raw_content, content_type) + normalized_content, validation = DhmvAcquisitionService._normalize_raster(coverage_content, scope_geometry, prepared) + acquired_at = datetime.now(UTC) + dataset = DatasetService.import_raster_bytes( + db, + project_id=project_id, + area_id=payload.area_id, + filename=filename, + content=normalized_content, + source=f"Digitaal Vlaanderen WCS {product.coverage_id}", + source_name=DhmvAcquisitionService.PROVIDER, + temporal_series_key=f"digitaal-vlaanderen:dhmvii:{product.key}:{prepared['request_hash'][:24]}", + observed_at=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC), + valid_from=datetime(2013, 1, 1, tzinfo=UTC), + valid_to=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC), + temporal_granularity="period", + source_version=DhmvAcquisitionService.SOURCE_VERSION, + content_type="image/tiff", + source_metadata={ + "provider": DhmvAcquisitionService.PROVIDER, + "service": "WCS", + "service_version": "2.0.1", + "product_key": product.key, + "product_display_name": product.display_name, + "surface_model": product.surface_model, + "coverage_id": product.coverage_id, + "native_resolution_m": product.native_resolution_m, + "analysis_resolution_m": validation["resolution_m"], + "source_crs": DhmvAcquisitionService.SOURCE_CRS, + "vertical_reference": DhmvAcquisitionService.VERTICAL_REFERENCE, + "vertical_unit": "m", + "acquisition_period": DhmvAcquisitionService.ACQUISITION_PERIOD, + "observation_date_precision": "period", + "nodata_value": validation["nodata_value"], + "valid_pixel_count": validation["valid_pixel_count"], + "minimum_m_taw": validation["minimum_m_taw"], + "maximum_m_taw": validation["maximum_m_taw"], + "bbox_epsg4326": prepared["bbox_epsg4326"], + "bbox_epsg31370": prepared["bbox_epsg31370"], + "catalog_url": product.catalog_url, + "attribution": DhmvAcquisitionService.ATTRIBUTION, + "license_note": DhmvAcquisitionService.LICENSE_NOTE, + "theme": "elevation", + "coverage_scope": "municipality" if payload.area_id else "bounded_selection", + }, + provenance_metadata={ + "acquisition": "explicit_bounded_wcs_coverage", + "acquired_at": acquired_at.isoformat(), + "request_hash": prepared["request_hash"], + "request_url": prepared["request_url"], + "response_content_type": content_type, + "response_sha256": hashlib.sha256(raw_content).hexdigest(), + "coverage_sha256": hashlib.sha256(coverage_content).hexdigest(), + "normalized_sha256": hashlib.sha256(normalized_content).hexdigest(), + "bbox_epsg4326": prepared["bbox_epsg4326"], + "bbox_epsg31370": prepared["bbox_epsg31370"], + "requested_resolution_m": prepared["resolution_m"], + "clipped_to_area_id": str(payload.area_id) if payload.area_id else None, + "validation": validation, + "limitation_message": product.limitation_message, + "water_depth_available": False, + "water_volume_available": False, + }, + ) + return DhmvAcquisitionResult( + output_dataset_id=dataset.id, + reused=False, + provider=DhmvAcquisitionService.PROVIDER, + product_key=product.key, + display_name=product.display_name, + surface_model=product.surface_model, + coverage_id=product.coverage_id, + native_resolution_m=product.native_resolution_m, + resolution_m=validation["resolution_m"], + width=validation["width"], + height=validation["height"], + valid_pixel_count=validation["valid_pixel_count"], + nodata_value=validation["nodata_value"], + bbox_epsg4326=prepared["bbox_epsg4326"], + bbox_epsg31370=prepared["bbox_epsg31370"], + vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE, + acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD, + attribution=DhmvAcquisitionService.ATTRIBUTION, + limitation_message=product.limitation_message, + ).model_dump(mode="json") diff --git a/backend/app/services/terrain_analysis_service.py b/backend/app/services/terrain_analysis_service.py new file mode 100644 index 00000000..a756d585 --- /dev/null +++ b/backend/app/services/terrain_analysis_service.py @@ -0,0 +1,254 @@ +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.dhmv import TerrainMetric, TerrainSelectionRequest, TerrainSelectionResponse, TerrainSelectionSummary +from app.services.dhmv_acquisition_service import DhmvAcquisitionService + + +class TerrainAnalysisService: + UNSUPPORTED_METRICS = ["water_depth_m", "water_volume_m3"] + LIMITATION = ( + "Hoogte, reliëf en helling zijn afgeleid uit DHMV II. Afstroming vraagt bijkomende hydrologische modellering. " + "Waterdiepte en watervolume zijn niet beschikbaar uit DTM/DSM alleen." + ) + + @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 != DhmvAcquisitionService.PROVIDER: + raise AppError( + code="INVALID_TERRAIN_DATASET", + message="Terrain analysis requires a governed DHMV 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 DHMV raster file is unavailable", status_code=404) + return dataset + + @staticmethod + def _selection_geometry(db, project_id: UUID, payload: TerrainSelectionRequest): + 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="TERRAIN_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: TerrainSelectionRequest, + *, + settings: Settings | None = None, + ) -> dict: + resolved_settings = settings or get_settings() + dataset = TerrainAnalysisService._load_dataset(db, project_id, dataset_id) + selection_4326 = TerrainAnalysisService._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 terrain analysis", status_code=503) from exc + + source_metadata = dataset.source_metadata or {} + product_key = str(source_metadata.get("product_key") or "") + surface_model = str(source_metadata.get("surface_model") or "") + if product_key not in DhmvAcquisitionService._products() or surface_model not in {"terrain", "surface"}: + raise AppError(code="INVALID_TERRAIN_METADATA", message="DHMV product provenance is incomplete", status_code=409) + + try: + with rasterio.open(dataset.storage_path) as source: + if source.crs is None: + raise AppError(code="INVALID_DATASET_CRS", message="DHMV raster CRS is missing", status_code=409) + transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True) + selection_metric = shapely_transform(transformer.transform, selection_4326) + source_extent = box(*source.bounds) + analysis_geometry = selection_metric.intersection(source_extent) + if analysis_geometry.is_empty or analysis_geometry.area <= 0: + raise AppError( + code="TERRAIN_SELECTION_OUTSIDE_DATASET", + message="Selection does not overlap the persisted DHMV 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.dhmv_max_pixels: + raise AppError( + code="TERRAIN_SELECTION_TOO_LARGE", + message="Terrain analysis exceeds the configured raster cell limit", + details={"pixel_count": expected_cells, "max_pixels": resolved_settings.dhmv_max_pixels}, + status_code=422, + ) + clipped, clipped_transform = mask( + source, + [mapping(analysis_geometry)], + crop=True, + filled=False, + indexes=[1], + ) + elevation = np.ma.asarray(clipped[0], dtype="float64") + raw = elevation.filled(np.nan) + nodata = source.nodata + invalid = ~np.isfinite(raw) + if nodata is not None: + invalid |= raw == float(nodata) + selected_cells = geometry_mask( + [mapping(analysis_geometry)], + out_shape=elevation.shape, + transform=clipped_transform, + invert=True, + ) + valid_mask = selected_cells & ~np.ma.getmaskarray(elevation) & ~invalid + values = raw[valid_mask] + if values.size == 0: + raise AppError(code="TERRAIN_NO_VALID_DATA", message="No valid DHMV height cells occur in this selection", status_code=422) + + resolution_x = abs(float(source.res[0])) + resolution_y = abs(float(source.res[1])) + slope_values = np.asarray([], dtype="float64") + if raw.shape[0] >= 2 and raw.shape[1] >= 2: + surface = np.where(valid_mask, raw, np.nan) + gradient_y, gradient_x = np.gradient(surface, resolution_y, resolution_x) + slope = np.degrees(np.arctan(np.hypot(gradient_x, gradient_y))) + slope_values = slope[np.isfinite(slope) & valid_mask] + except AppError: + raise + except Exception as exc: + raise AppError( + code="TERRAIN_ANALYSIS_FAILED", + message="The persisted DHMV 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) -> TerrainMetric: + return TerrainMetric( + metric_key=key, + metric_label=label, + metric_value=round(float(value), 4), + metric_unit=unit, + aggregation_method=method, + ) + + prefix = "terrain" if surface_model == "terrain" else "surface" + elevation_label = "Gemiddelde maaiveldhoogte" if surface_model == "terrain" else "Gemiddelde oppervlaktehoogte" + metrics = [ + metric(f"{prefix}_elevation_mean_m", elevation_label, values.mean(), "m TAW", "mean_valid_cells"), + metric(f"{prefix}_elevation_min_m", "Laagste hoogte", values.min(), "m TAW", "minimum_valid_cells"), + metric(f"{prefix}_elevation_max_m", "Hoogste hoogte", values.max(), "m TAW", "maximum_valid_cells"), + metric(f"{prefix}_elevation_p10_m", "10e percentiel hoogte", np.percentile(values, 10), "m TAW", "percentile_10_valid_cells"), + metric(f"{prefix}_elevation_p90_m", "90e percentiel hoogte", np.percentile(values, 90), "m TAW", "percentile_90_valid_cells"), + metric("relief_m", "Reliëfverschil", values.max() - values.min(), "m", "maximum_minus_minimum"), + ] + if slope_values.size: + metrics.extend( + [ + metric("slope_mean_deg", "Gemiddelde helling", slope_values.mean(), "°", "mean_finite_gradient"), + metric("slope_p90_deg", "90e percentiel helling", np.percentile(slope_values, 90), "°", "percentile_90_finite_gradient"), + metric("slope_max_deg", "Steilste helling", slope_values.max(), "°", "maximum_finite_gradient"), + ] + ) + primary = metrics[0] + selected_cell_count = int(selected_cells.sum()) + response = TerrainSelectionResponse( + dataset_id=dataset.id, + product_key=product_key, + surface_model=surface_model, + selection_bbox=payload.bbox, + selection_area_id=payload.area_id, + sample_count=int(values.size), + slope_sample_count=int(slope_values.size), + coverage_ratio=round(float(values.size / max(1, selected_cell_count)), 6), + resolution_m=round(max(resolution_x, resolution_y), 4), + vertical_reference=str(source_metadata.get("vertical_reference") or DhmvAcquisitionService.VERTICAL_REFERENCE), + summary=TerrainSelectionSummary( + 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=TerrainAnalysisService.UNSUPPORTED_METRICS, + limitation_message=TerrainAnalysisService.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 = TerrainAnalysisService._load_dataset(db, project_id, dataset_id) + 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 terrain 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 not valid.any(): + raise AppError(code="TERRAIN_NO_VALID_DATA", message="DHMV 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) + stops = np.asarray([0.0, 0.25, 0.5, 0.75, 1.0]) + colors = np.asarray( + [ + [30, 94, 91], + [79, 139, 102], + [194, 183, 105], + [173, 121, 79], + [105, 94, 108], + ], + 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, 225, 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="TERRAIN_PREVIEW_FAILED", + message="The persisted DHMV raster could not be rendered", + details={"reason": str(exc)}, + status_code=500, + ) from exc diff --git a/backend/tests/test_sprint205_dhmv_terrain.py b/backend/tests/test_sprint205_dhmv_terrain.py new file mode 100644 index 00000000..d8db9946 --- /dev/null +++ b/backend/tests/test_sprint205_dhmv_terrain.py @@ -0,0 +1,404 @@ +from __future__ import annotations + +from pathlib import Path +from uuid import uuid4 + +import numpy as np +import pytest +import rasterio +from fastapi.testclient import TestClient +from geoalchemy2.shape import from_shape +from pyproj import Transformer +from rasterio.io import MemoryFile +from rasterio.transform import from_origin +from shapely.geometry import MultiPolygon, box + +from app.core.config import Settings +from app.core.errors import AppError +from app.db.session import get_db +from app.main import app +from app.models import Area, Dataset, DatasetVersion, Job, Project +from app.schemas.dhmv import DhmvAcquireRequest, TerrainSelectionRequest +from app.services.dhmv_acquisition_service import DhmvAcquisitionService +from app.services.terrain_analysis_service import TerrainAnalysisService + + +ROOT = Path(__file__).resolve().parents[2] + + +class FakeQuery: + def __init__(self, result=None): + self.result = result + + def filter(self, *_args): + return self + + def order_by(self, *_args): + return self + + def first(self): + return self.result + + +class FakeSession: + def __init__(self, rows=None, query_result=None): + self.rows = rows or {} + self.query_result = query_result + self.added = [] + + def get(self, model, row_id): + row = self.rows.get((model, row_id)) + if row is not None: + return row + return next((item for item in self.added if isinstance(item, model) and item.id == row_id), None) + + def add(self, row): + self.added.append(row) + + def commit(self): + return None + + def rollback(self): + return None + + def refresh(self, row): + return row + + def query(self, _model): + return FakeQuery(self.query_result) + + +class FakeResponse: + def __init__(self, content: bytes, content_type: str): + self.content = content + self.headers = {"Content-Type": content_type, "Content-Length": str(len(content))} + + def __enter__(self): + return self + + def __exit__(self, *_args): + return None + + def read(self, limit: int): + return self.content[:limit] + + +def lambert_bbox_payload(*, side_m: float = 100.0, product_key: str = "dtm_1m", area_id=None) -> DhmvAcquireRequest: + west, south = 200_000.0, 210_000.0 + transformer = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) + min_x, min_y = transformer.transform(west, south) + max_x, max_y = transformer.transform(west + side_m, south + side_m) + return DhmvAcquireRequest( + bbox={"min_x": min_x, "min_y": min_y, "max_x": max_x, "max_y": max_y, "crs": "EPSG:4326"}, + area_id=area_id, + product_key=product_key, + resolution_m=5.0, + force_refresh=True, + ) + + +def elevation_tiff(*, left: float, top: float, width: int, height: int, resolution: float = 5.0) -> bytes: + rows, columns = np.indices((height, width)) + values = (20.0 + columns * 0.5 + rows * 1.0).astype("float32") + with MemoryFile() as memory: + with memory.open( + driver="GTiff", + width=width, + height=height, + count=1, + dtype="float32", + crs="EPSG:31370", + transform=from_origin(left, top, resolution, resolution), + nodata=-9999.0, + ) as output: + output.write(values, 1) + return memory.read() + + +def multipart_tiff(content: bytes) -> tuple[bytes, str]: + boundary = "wcs-test" + payload = ( + f"--{boundary}\r\nContent-Type: text/xml\r\nContent-ID: GML-Part\r\n\r\n\r\n" + f"--{boundary}\r\nContent-Type: image/tiff\r\nContent-ID: coverage.tif\r\n\r\n" + ).encode() + content + f"\r\n--{boundary}--\r\n".encode() + return payload, f'multipart/mixed; boundary="{boundary}"' + + +def test_dhmv_registry_is_governed_and_semantically_explicit() -> None: + products = DhmvAcquisitionService.list_products() + + assert [item["key"] for item in products] == ["dtm_1m", "dsm_1m"] + assert {item["coverage_id"] for item in products} == {"DHMVII_DTM_1m", "DHMVII_DSM_1m"} + assert all(item["native_resolution_m"] == 1.0 for item in products) + assert all(item["source_crs"] == "EPSG:31370" for item in products) + assert all("TAW" in item["vertical_reference"] for item in products) + assert all(item["acquisition_period"] == "2013-2015" for item in products) + assert "waterdiepte" in products[0]["limitation_message"] + + +def test_dhmv_request_uses_bounded_official_wcs_scaling() -> None: + prepared = DhmvAcquisitionService._prepared_request(lambert_bbox_payload(), Settings(_env_file=None)) + + assert prepared["coverage_id"] == "DHMVII_DTM_1m" + assert prepared["params"]["SCALEFACTOR"] == "5" + assert prepared["params"]["SUBSET"][0].startswith("x(") + assert prepared["params"]["SUBSET"][1].startswith("y(") + assert "geo.api.vlaanderen.be%2FDHMV" not in prepared["request_url"] + assert prepared["request_url"].startswith("https://geo.api.vlaanderen.be/DHMV/wcs?") + assert prepared["width"] * prepared["height"] <= 12_000_000 + assert len(prepared["request_hash"]) == 64 + + with pytest.raises(AppError) as exc_info: + DhmvAcquisitionService._prepared_request( + lambert_bbox_payload(product_key="arbitrary"), + Settings(_env_file=None), + ) + assert exc_info.value.code == "DHMV_PRODUCT_NOT_SUPPORTED" + + +def test_dhmv_request_rejects_unsafe_size_and_resolution() -> None: + with pytest.raises(AppError) as exc_info: + DhmvAcquisitionService._prepared_request(lambert_bbox_payload(side_m=5.0), Settings(_env_file=None)) + assert exc_info.value.code == "DHMV_SELECTION_TOO_SMALL" + + payload = lambert_bbox_payload() + payload.resolution_m = 0.5 + with pytest.raises(Exception): + DhmvAcquireRequest.model_validate(payload.model_dump()) + + +def test_dhmv_multipart_geotiff_is_extracted_and_invalid_response_fails_closed() -> None: + tiff = elevation_tiff(left=200_000, top=210_100, width=20, height=20) + multipart, content_type = multipart_tiff(tiff) + + assert DhmvAcquisitionService._extract_geotiff(multipart, content_type) == tiff + with pytest.raises(AppError) as exc_info: + DhmvAcquisitionService._extract_geotiff(b"", "text/xml") + assert exc_info.value.code == "DHMV_PROVIDER_INVALID_RESPONSE" + + +def test_dhmv_acquisition_clips_validates_and_persists_via_dataset_service(tmp_path) -> None: + project_id = uuid4() + area_id = uuid4() + payload = lambert_bbox_payload(area_id=area_id) + area_geometry = MultiPolygon([box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y)]) + db = FakeSession( + { + (Project, project_id): Project(id=project_id, name="Mol"), + (Area, area_id): Area( + id=area_id, + project_id=project_id, + name="Gemeente Mol", + geometry=from_shape(area_geometry, srid=4326), + ), + } + ) + settings = Settings(_env_file=None, storage_root=str(tmp_path), dhmv_resolution_m=5.0) + prepared = DhmvAcquisitionService._prepared_request(payload, settings) + tiff = elevation_tiff( + left=prepared["bbox_epsg31370"][0], + top=prepared["bbox_epsg31370"][3], + width=prepared["width"], + height=prepared["height"], + ) + multipart, content_type = multipart_tiff(tiff) + + result = DhmvAcquisitionService.acquire( + db, + project_id, + payload, + settings=settings, + opener=lambda *_args, **_kwargs: FakeResponse(multipart, content_type), + ) + + dataset = next(item for item in db.added if isinstance(item, Dataset)) + version = next(item for item in db.added if isinstance(item, DatasetVersion)) + assert result["output_dataset_id"] == str(dataset.id) + assert dataset.source_name == "digitaal_vlaanderen_dhmv" + assert dataset.area_id == area_id + assert dataset.dataset_type == "raster" + assert dataset.crs == "EPSG:31370" + assert dataset.checksum_sha256 == version.checksum_sha256 + assert dataset.source_metadata["surface_model"] == "terrain" + assert dataset.source_metadata["native_resolution_m"] == 1.0 + assert dataset.source_metadata["analysis_resolution_m"] == 5.0 + assert dataset.source_metadata["nodata_value"] == -9999.0 + assert dataset.provenance_metadata["water_depth_available"] is False + assert dataset.provenance_metadata["water_volume_available"] is False + assert len(dataset.provenance_metadata["response_sha256"]) == 64 + with rasterio.open(dataset.storage_path) as stored: + assert stored.crs.to_epsg() == 31370 + assert stored.count == 1 + assert stored.nodata == -9999.0 + assert stored.res == pytest.approx((5.0, 5.0)) + + +def test_terrain_analysis_returns_governed_elevation_relief_and_slope(tmp_path) -> None: + project_id = uuid4() + dataset_id = uuid4() + path = tmp_path / "terrain.tif" + path.write_bytes(elevation_tiff(left=200_000, top=210_100, width=20, height=20)) + to_wgs84 = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) + min_x, min_y = to_wgs84.transform(200_000, 210_000) + max_x, max_y = to_wgs84.transform(200_100, 210_100) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="dhmvii_terrain_5m.tif", + dataset_type="raster", + source="official WCS", + source_name="digitaal_vlaanderen_dhmv", + source_metadata={ + "product_key": "dtm_1m", + "surface_model": "terrain", + "vertical_reference": "TAW (Tweede Algemene Waterpassing)", + }, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + payload = TerrainSelectionRequest( + bbox={"min_x": min_x, "min_y": min_y, "max_x": max_x, "max_y": max_y, "crs": "EPSG:4326"} + ) + + result = TerrainAnalysisService.analyze(db, project_id, dataset_id, payload, settings=Settings(_env_file=None)) + metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]} + + assert result["sample_count"] > 300 + assert result["coverage_ratio"] > 0.99 + assert result["resolution_m"] == 5.0 + assert result["summary"]["metric_unit"] == "m TAW" + assert metrics["relief_m"]["metric_value"] > 20 + assert metrics["slope_mean_deg"]["metric_value"] == pytest.approx(12.6044, abs=0.01) + assert result["unsupported_metrics"] == ["water_depth_m", "water_volume_m3"] + assert "Waterdiepte" in result["limitation_message"] + + +def test_terrain_analysis_rejects_non_dhmv_raster(tmp_path) -> None: + project_id = uuid4() + dataset_id = uuid4() + path = tmp_path / "other.tif" + path.write_bytes(elevation_tiff(left=200_000, top=210_100, width=20, height=20)) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="other.tif", + dataset_type="raster", + source="manual", + source_name="manual", + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + + with pytest.raises(AppError) as exc_info: + TerrainAnalysisService.analyze(db, project_id, dataset_id, TerrainSelectionRequest(bbox=lambert_bbox_payload().bbox)) + assert exc_info.value.code == "INVALID_TERRAIN_DATASET" + + +def test_terrain_renderer_returns_browser_png(tmp_path) -> None: + project_id = uuid4() + dataset_id = uuid4() + path = tmp_path / "terrain.tif" + path.write_bytes(elevation_tiff(left=200_000, top=210_100, width=20, height=20)) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="terrain.tif", + dataset_type="raster", + source="official", + source_name="digitaal_vlaanderen_dhmv", + source_metadata={"product_key": "dtm_1m", "surface_model": "terrain"}, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + + assert TerrainAnalysisService.render_png(db, project_id, dataset_id).startswith(b"\x89PNG\r\n\x1a\n") + + +def test_dhmv_endpoints_use_canonical_envelopes(monkeypatch) -> None: + project_id = uuid4() + output_dataset_id = uuid4() + db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")}) + monkeypatch.setattr( + DhmvAcquisitionService, + "acquire", + lambda *_args, **_kwargs: { + "output_dataset_id": str(output_dataset_id), + "provider": "digitaal_vlaanderen_dhmv", + "reused": False, + }, + ) + monkeypatch.setattr( + TerrainAnalysisService, + "analyze", + lambda *_args, **_kwargs: { + "dataset_id": str(output_dataset_id), + "sample_count": 100, + "summary": {"metric_value": 25.0, "metric_unit": "m TAW", "metrics": []}, + "unsupported_metrics": ["water_depth_m", "water_volume_m3"], + }, + ) + app.dependency_overrides[get_db] = lambda: db + try: + products = TestClient(app).get(f"/api/v1/projects/{project_id}/datasets/dhmv/products") + acquisition = TestClient(app).post( + f"/api/v1/projects/{project_id}/datasets/dhmv/acquire", + json=lambert_bbox_payload().model_dump(mode="json"), + ) + terrain = TestClient(app).post( + f"/api/v1/projects/{project_id}/datasets/{output_dataset_id}/raster/terrain/select", + json={"bbox": lambert_bbox_payload().bbox.model_dump()}, + ) + finally: + app.dependency_overrides.clear() + + assert products.status_code == 200 + assert set(products.json()) == {"data"} + assert products.json()["data"]["total"] == 2 + assert acquisition.status_code == 200 + assert set(acquisition.json()) == {"data"} + assert acquisition.json()["data"]["job_type"] == "raster.dhmv.acquire" + assert acquisition.json()["data"]["output_dataset_id"] == str(output_dataset_id) + assert terrain.status_code == 200 + assert set(terrain.json()) == {"data"} + assert terrain.json()["data"]["sample_count"] == 100 + assert terrain.json()["data"]["unsupported_metrics"] == ["water_depth_m", "water_volume_m3"] + assert any(isinstance(item, Job) for item in db.added) + + +def test_frontend_and_runtime_expose_dhmv_workflow() -> None: + app_source = (ROOT / "frontend" / "src" / "App.tsx").read_text(encoding="utf-8") + map_source = (ROOT / "frontend" / "src" / "components" / "map" / "MapWorkspace.tsx").read_text(encoding="utf-8") + hook_source = (ROOT / "frontend" / "src" / "hooks" / "useMapSelectionExtract.ts").read_text(encoding="utf-8") + service_source = (ROOT / "frontend" / "src" / "services" / "api" / "datasets.ts").read_text(encoding="utf-8") + + assert "digitaal_vlaanderen_dhmv" in app_source + assert "Hoogte & reliëf" in map_source + assert "terrainImageUrl" in map_source + assert "selectTerrain" in hook_source + assert "/raster/terrain/select" in service_source + for path in ( + ROOT / ".env.example", + ROOT / "docker-compose.yml", + ROOT / "docker-compose.unraid.yml", + ROOT / "deploy" / "unraid" / "run-dockerman-container.sh", + ROOT / "deploy" / "unraid" / "geointel-unraid-template.xml", + ): + content = path.read_text(encoding="utf-8") + assert "DHMV_ENABLED" in content + assert "DHMV_RESOLUTION_M" in content + assert "DHMV_MAX_PIXELS" in content + + +def test_dhmv_operator_is_packaged_and_release_checked() -> None: + operator = (ROOT / "scripts" / "provision_mol_dhmv.py").read_text(encoding="utf-8") + readiness = (ROOT / "scripts" / "run_readiness_check.sh").read_text(encoding="utf-8") + dockerfile = (ROOT / "backend" / "Dockerfile").read_text(encoding="utf-8") + + assert "/datasets/dhmv/acquire" in operator + assert "/raster/terrain/select" in operator + assert "water_depth_m" in operator + assert "py_compile scripts/provision_mol_dhmv.py" in readiness + assert "COPY . /app" in dockerfile diff --git a/deploy/unraid/geointel-unraid-template.xml b/deploy/unraid/geointel-unraid-template.xml index 6ef56a45..61f69dc4 100644 --- a/deploy/unraid/geointel-unraid-template.xml +++ b/deploy/unraid/geointel-unraid-template.xml @@ -34,6 +34,11 @@ https://geo.api.vlaanderen.be/OMWRGBMRVL/wms 1.0 1024 + true + https://geo.api.vlaanderen.be/DHMV/wcs + 5.0 + 20000 + 12000000 true http://host.docker.internal:11434 qwen3.5:9b diff --git a/deploy/unraid/geointel.env.example b/deploy/unraid/geointel.env.example index 1a18e6b8..536ece29 100644 --- a/deploy/unraid/geointel.env.example +++ b/deploy/unraid/geointel.env.example @@ -32,6 +32,14 @@ ORTHOPHOTO_RESOLUTION_M=1.0 ORTHOPHOTO_MIN_SIDE_M=128 ORTHOPHOTO_MAX_SIDE_M=1024 ORTHOPHOTO_CACHE_TTL_HOURS=24 +DHMV_ENABLED=true +DHMV_WCS_URL=https://geo.api.vlaanderen.be/DHMV/wcs +DHMV_RESOLUTION_M=5.0 +DHMV_MIN_SIDE_M=10 +DHMV_MAX_SIDE_M=20000 +DHMV_MAX_PIXELS=12000000 +DHMV_TIMEOUT_SECONDS=300 +DHMV_MAX_RESPONSE_MB=160 # Optional configured-YOLO runtime. Keep disabled unless a local model is mounted. GEOINTEL_INSTALL_AI=false diff --git a/deploy/unraid/run-dockerman-container.sh b/deploy/unraid/run-dockerman-container.sh index a733b8aa..f66b2ae4 100644 --- a/deploy/unraid/run-dockerman-container.sh +++ b/deploy/unraid/run-dockerman-container.sh @@ -27,6 +27,14 @@ ORTHOPHOTO_RESOLUTION_M="${ORTHOPHOTO_RESOLUTION_M:-1.0}" ORTHOPHOTO_MIN_SIDE_M="${ORTHOPHOTO_MIN_SIDE_M:-128}" ORTHOPHOTO_MAX_SIDE_M="${ORTHOPHOTO_MAX_SIDE_M:-1024}" ORTHOPHOTO_CACHE_TTL_HOURS="${ORTHOPHOTO_CACHE_TTL_HOURS:-24}" +DHMV_ENABLED="${DHMV_ENABLED:-true}" +DHMV_WCS_URL="${DHMV_WCS_URL:-https://geo.api.vlaanderen.be/DHMV/wcs}" +DHMV_RESOLUTION_M="${DHMV_RESOLUTION_M:-5.0}" +DHMV_MIN_SIDE_M="${DHMV_MIN_SIDE_M:-10}" +DHMV_MAX_SIDE_M="${DHMV_MAX_SIDE_M:-20000}" +DHMV_MAX_PIXELS="${DHMV_MAX_PIXELS:-12000000}" +DHMV_TIMEOUT_SECONDS="${DHMV_TIMEOUT_SECONDS:-300}" +DHMV_MAX_RESPONSE_MB="${DHMV_MAX_RESPONSE_MB:-160}" YOLO_ENABLED="${YOLO_ENABLED:-false}" YOLO_MODELS_DIR="${YOLO_MODELS_DIR:-/app/models}" YOLO_MODEL_PATH="${YOLO_MODEL_PATH:-}" @@ -103,6 +111,14 @@ docker run -d \ -e ORTHOPHOTO_MIN_SIDE_M="$ORTHOPHOTO_MIN_SIDE_M" \ -e ORTHOPHOTO_MAX_SIDE_M="$ORTHOPHOTO_MAX_SIDE_M" \ -e ORTHOPHOTO_CACHE_TTL_HOURS="$ORTHOPHOTO_CACHE_TTL_HOURS" \ + -e DHMV_ENABLED="$DHMV_ENABLED" \ + -e DHMV_WCS_URL="$DHMV_WCS_URL" \ + -e DHMV_RESOLUTION_M="$DHMV_RESOLUTION_M" \ + -e DHMV_MIN_SIDE_M="$DHMV_MIN_SIDE_M" \ + -e DHMV_MAX_SIDE_M="$DHMV_MAX_SIDE_M" \ + -e DHMV_MAX_PIXELS="$DHMV_MAX_PIXELS" \ + -e DHMV_TIMEOUT_SECONDS="$DHMV_TIMEOUT_SECONDS" \ + -e DHMV_MAX_RESPONSE_MB="$DHMV_MAX_RESPONSE_MB" \ -e YOLO_ENABLED="$YOLO_ENABLED" \ -e YOLO_MODELS_DIR="$YOLO_MODELS_DIR" \ -e YOLO_MODEL_PATH="$YOLO_MODEL_PATH" \ diff --git a/docker-compose.unraid.yml b/docker-compose.unraid.yml index 770a3632..f068e25f 100644 --- a/docker-compose.unraid.yml +++ b/docker-compose.unraid.yml @@ -25,6 +25,14 @@ services: ORTHOPHOTO_MIN_SIDE_M: ${ORTHOPHOTO_MIN_SIDE_M:-128} ORTHOPHOTO_MAX_SIDE_M: ${ORTHOPHOTO_MAX_SIDE_M:-1024} ORTHOPHOTO_CACHE_TTL_HOURS: ${ORTHOPHOTO_CACHE_TTL_HOURS:-24} + DHMV_ENABLED: ${DHMV_ENABLED:-true} + DHMV_WCS_URL: ${DHMV_WCS_URL:-https://geo.api.vlaanderen.be/DHMV/wcs} + DHMV_RESOLUTION_M: ${DHMV_RESOLUTION_M:-5.0} + DHMV_MIN_SIDE_M: ${DHMV_MIN_SIDE_M:-10} + DHMV_MAX_SIDE_M: ${DHMV_MAX_SIDE_M:-20000} + DHMV_MAX_PIXELS: ${DHMV_MAX_PIXELS:-12000000} + DHMV_TIMEOUT_SECONDS: ${DHMV_TIMEOUT_SECONDS:-300} + DHMV_MAX_RESPONSE_MB: ${DHMV_MAX_RESPONSE_MB:-160} YOLO_ENABLED: ${YOLO_ENABLED:-false} YOLO_MODELS_DIR: ${YOLO_MODELS_DIR:-/app/models} YOLO_MODEL_PATH: ${YOLO_MODEL_PATH:-} diff --git a/docker-compose.yml b/docker-compose.yml index 1c3eab88..7845b3f6 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -30,6 +30,14 @@ services: ORTHOPHOTO_MIN_SIDE_M: ${ORTHOPHOTO_MIN_SIDE_M:-128} ORTHOPHOTO_MAX_SIDE_M: ${ORTHOPHOTO_MAX_SIDE_M:-1024} ORTHOPHOTO_CACHE_TTL_HOURS: ${ORTHOPHOTO_CACHE_TTL_HOURS:-24} + DHMV_ENABLED: ${DHMV_ENABLED:-true} + DHMV_WCS_URL: ${DHMV_WCS_URL:-https://geo.api.vlaanderen.be/DHMV/wcs} + DHMV_RESOLUTION_M: ${DHMV_RESOLUTION_M:-5.0} + DHMV_MIN_SIDE_M: ${DHMV_MIN_SIDE_M:-10} + DHMV_MAX_SIDE_M: ${DHMV_MAX_SIDE_M:-20000} + DHMV_MAX_PIXELS: ${DHMV_MAX_PIXELS:-12000000} + DHMV_TIMEOUT_SECONDS: ${DHMV_TIMEOUT_SECONDS:-300} + DHMV_MAX_RESPONSE_MB: ${DHMV_MAX_RESPONSE_MB:-160} YOLO_ENABLED: ${YOLO_ENABLED:-false} YOLO_MODELS_DIR: ${YOLO_MODELS_DIR:-/app/models} YOLO_MODEL_PATH: ${YOLO_MODEL_PATH:-} diff --git a/docs/API_CONTRACTS.md b/docs/API_CONTRACTS.md index 018d47fe..ca81eccc 100644 --- a/docs/API_CONTRACTS.md +++ b/docs/API_CONTRACTS.md @@ -219,6 +219,49 @@ dimensions, EPSG:4326/EPSG:31370 bounds, sampling resolution, attribution, cache reuse and limitation text. Historical products also persist their observation/validity period and a spatially scoped temporal-series key. +### GET `/api/v1/projects/{project_id}/datasets/dhmv/products` + +Returns the fixed official DHMV II product registry in the canonical envelope. +The registry contains only `dtm_1m` (`DHMVII_DTM_1m`) and `dsm_1m` +(`DHMVII_DSM_1m`). Each item records native 1 m resolution, EPSG:31370, TAW, +the 2013-2015 acquisition period, attribution, catalogue and limitations. + +### POST `/api/v1/projects/{project_id}/datasets/dhmv/acquire` + +Runs a bounded WCS 2.0.1 `GetCoverage` request behind the existing synchronous +Job abstraction. Arbitrary coverage identifiers are rejected. + +```json +{ + "bbox": {"min_x": 5.0, "min_y": 51.1, "max_x": 5.2, "max_y": 51.3, "crs": "EPSG:4326"}, + "area_id": "optional-project-area-uuid", + "product_key": "dtm_1m", + "resolution_m": 5.0, + "force_refresh": false +} +``` + +The service extracts the GeoTIFF from the official multipart response, clips +to the exact persisted Area when supplied, validates EPSG:31370, one band, +resolution, nodata and valid cells, then persists through `DatasetService`. +The default 5 m file is an analysis copy of the retained 1 m source product; +both resolutions and all request/response/output checksums remain provenance. + +### POST `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/terrain/select` + +Accepts an EPSG:4326 rectangle and optional Area id. It reads only a governed, +ready DHMV Dataset and returns a canonical envelope with valid-cell coverage, +mean/min/max/P10/P90 height in `m TAW`, relief in metres and mean/P90/max slope +in degrees. Area geometry is an exact mask, not only a bounding box. + +The response always lists `water_depth_m` and `water_volume_m3` under +`unsupported_metrics`. Drainage is not calculated by this endpoint. + +### GET `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/terrain/image` + +Returns a browser-safe PNG colour relief for the persisted governed DHMV +Dataset. This binary MapLibre source never accepts an arbitrary file path. + Safety contract: - every side must measure between 128 m and 1,024 m in EPSG:31370; diff --git a/docs/CODEX_EXECUTION_LOG.md b/docs/CODEX_EXECUTION_LOG.md index 2fb51ccc..29f280f9 100644 --- a/docs/CODEX_EXECUTION_LOG.md +++ b/docs/CODEX_EXECUTION_LOG.md @@ -8684,3 +8684,36 @@ Next: - Implement P4 Digitaal Hoogtemodel Vlaanderen with a governed DTM/DSM product, vertical reference, bounded raster storage and measured elevation/slope statistics. Do not infer water depth or volume from terrain height alone. + +## Sprint 207 - Governed DHMV II terrain foundation (2026-07-15) + +Implemented: +- Added a fixed `dtm_1m`/`dsm_1m` product registry for the official Digitaal + Vlaanderen production WCS and rejected arbitrary coverage identifiers. +- Added bounded WCS scaling, multipart GeoTIFF extraction, exact persisted-Area + clipping and validation of EPSG:31370, one Float32 band, resolution, nodata + and valid cells before canonical DatasetService persistence. +- Retained response, coverage and normalized-output SHA256 evidence plus native + 1 m resolution, default 5 m analysis resolution, TAW and period 2013-2015. +- Added exact masked terrain selection with governed height, relief and slope + metrics. DTM and DSM semantics remain separate; water depth/volume remain + explicit unsupported metrics and drainage is not calculated. +- Added the `Hoogte & reliëf` map theme, MapLibre colour-relief overlay, + readable period/source labels and raster-aware rectangle/full-Area queries. +- Added `scripts/provision_mol_dhmv.py`, Docker/Unraid settings and documentation. + +Initial validation: +- Focused acquisition, raster, GIS metric, API, frontend and packaging tests + pass with deterministic synthetic GeoTIFF fixtures. +- No migration, direct raster database write, LiDAR point-cloud processing or + water-volume inference was introduced. + +Known limitations: +- DHMV II represents acquisition period 2013-2015 and is not a current or + annual height series. Exact flight-day contours are not yet joined. +- The 5 m analysis copy is operationally bounded; sub-5 m detail requires an + explicit smaller acquisition. DSM-DTM building height remains future work. + +Next: +- Complete the full readiness and live Mol operator/browser validation, then + audit P5 bathymetry sources before enabling any water depth or volume metric. diff --git a/docs/DATA_SOURCES.md b/docs/DATA_SOURCES.md index afe8fc32..7be8570e 100644 --- a/docs/DATA_SOURCES.md +++ b/docs/DATA_SOURCES.md @@ -220,9 +220,10 @@ until a governed operator import, provenance record and validation pass exist. - Buildings and Addresses Register (Digitaal Vlaanderen): continuously updated building status, life cycle and address linkage; complementary to GRB geometry and not yet imported. -- DHMV II DTM/DSM (Digitaal Vlaanderen): 1 m/5 m elevation based on 2013-2015 - LiDAR, suitable for elevation, slope and drainage. It does not provide water - depth. +- DHMV II DTM/DSM (Digitaal Vlaanderen): governed bounded WCS acquisition is + implemented below. The official 1 m source and 5 m analysis copy are suitable + for measured elevation, relief and slope. Drainage remains an interpretation; + the products do not provide water depth. ## Governed BWK and Natura 2000 state 2025 @@ -322,10 +323,22 @@ quality metrics. - Naam: Digitaal Hoogtemodel Vlaanderen - Type: raster/hoogte -- Gebruik: DEM, DSM, helling, laagste punten -- Toegang: Vlaamse open data, download/WCS nader te bepalen -- Cache: raster storage -- Prioriteit: V3 +- Gebruik: DTM/DSM hoogte, reliëf, helling en laagste punten +- Toegang: productie-WCS `https://geo.api.vlaanderen.be/DHMV/wcs` +- Coverages: `DHMVII_DTM_1m`, `DHMVII_DSM_1m` +- Native raster: 1 m Float32, EPSG:31370, nodata `-9999`, hoogte in TAW +- Opnameperiode: 2013-2015; geen uniforme recente peildatum +- Cache: canonical raster Dataset plus WCS request/response/output checksums +- Operator: `scripts/provision_mol_dhmv.py` +- Prioriteit: P4 uitgevoerd voor Mol + +The operator requests a bounded 5 m analysis copy by default so a complete +municipality remains operationally manageable while retaining the official +1 m native resolution in source metadata. The response is clipped to the exact +persisted Area before `DatasetService` stores it. DTM is bare-earth terrain; +DSM includes buildings and vegetation. Neither is exposed as water depth, +water volume or a directly measured building-height product. A drainage model +would require a separately governed hydrological processing pass. ## Gebouwenregister diff --git a/docs/DATA_SPECIFICATION.md b/docs/DATA_SPECIFICATION.md index dbe17977..c56ce8d5 100644 --- a/docs/DATA_SPECIFICATION.md +++ b/docs/DATA_SPECIFICATION.md @@ -159,12 +159,23 @@ Hoogtedata voor terrein- en watergevoeligheidsanalyse. ### Gebruik -- DEM -- DSM -- helling -- laagste punten -- hoogteprofiel -- gebouwhoogte-inschatting indien DSM + gebouwpolygonen beschikbaar zijn +- DTM-maaiveldhoogte in meter TAW +- DSM-oppervlaktehoogte in meter TAW +- reliëf en helling in graden uit geldige rastercellen +- laagste punten als terreinindicator +- gebouwhoogte-inschatting alleen in een latere, gevalideerde DSM-DTM/gebouwketen + +### Governed product + +- Digitaal Vlaanderen DHMV II, DTM/DSM raster 1 m +- WCS coverages `DHMVII_DTM_1m` en `DHMVII_DSM_1m` +- EPSG:31370, Float32, nodata `-9999`, verticale referentie TAW +- opnameperiode 2013-2015 +- standaard GeoIntel-analysekopie 5 m, exact geclipt op Area + +Waterdiepte, waterinhoud, actuele toestand en afstroming zijn geen directe +DHMV-metingen. Zij blijven onbeschikbaar tot een afzonderlijke bron en +gevalideerde methode bestaan. ### Type @@ -172,7 +183,7 @@ Raster / afgeleid van LiDAR. ### Prioriteit -V3. +P4 operationeel voor Mol; regionale uitrol volgt dezelfde operatorgrenzen. ## LAS/LAZ LiDAR diff --git a/docs/STORAGE_ARCHITECTURE.md b/docs/STORAGE_ARCHITECTURE.md index 5504c2b2..4a457399 100644 --- a/docs/STORAGE_ARCHITECTURE.md +++ b/docs/STORAGE_ARCHITECTURE.md @@ -109,6 +109,15 @@ product/layer, request/spatial hash, temporal validity and limitations are held in source/provenance metadata. Browser PNG rendering is derived on request and does not replace the stored GeoTIFF. +DHMV II DTM/DSM outputs are also normal raster Dataset files. The provider WCS +returns multipart coverage data; GeoIntel retains response and extracted +coverage SHA256 values in provenance, then stores one normalized, compressed, +Area-clipped GeoTIFF with its ordinary Dataset/DatasetVersion checksum. Source +metadata records the 1 m native product, 5 m default analysis grid, EPSG:31370, +`-9999` nodata, TAW and acquisition period 2013-2015. Colour-relief PNGs are +derived browser views and are never authoritative. No raster binary is stored +in PostgreSQL and no DHMV file is treated as water depth or volume. + BWK/Natura 2000 evidence lives under `storage/operator-evidence/bwk-natura2000-2025/mol/`. The `raw/` directory contains immutable WFS pages; the adjacent manifest records their URLs, diff --git a/docs/TODO.md b/docs/TODO.md index 3385175e..41d40d92 100644 --- a/docs/TODO.md +++ b/docs/TODO.md @@ -60,10 +60,10 @@ pass live Mol validation before regional expansion. ### P4 - Digitaal Hoogtemodel Vlaanderen -- [ ] Select official DTM/DSM product, service and native resolution; retain acquisition date and vertical reference. -- [ ] Add bounded raster acquisition/clip storage with nodata, CRS, resolution and checksum validation. -- [ ] Implement governed elevation, relief and slope statistics in metres/degrees; keep drainage interpretation explicitly derived. -- [ ] Do not label terrain/surface height as water depth and do not enable volume from DHMV alone. +- [x] Select official DTM/DSM product, service and native resolution; retain acquisition date and vertical reference. +- [x] Add bounded raster acquisition/clip storage with nodata, CRS, resolution and checksum validation. +- [x] Implement governed elevation, relief and slope statistics in metres/degrees; keep drainage interpretation explicitly derived. +- [x] Do not label terrain/surface height as water depth and do not enable volume from DHMV alone. ### P5 - Water depth / bathymetry diff --git a/frontend/README.md b/frontend/README.md index b72403f1..3598cb2d 100644 --- a/frontend/README.md +++ b/frontend/README.md @@ -492,6 +492,14 @@ another municipality or the complete Kempen scope falls back to the complete regional GRB building layer. This avoids presenting a Mol-only snapshot as regional coverage. +When a governed DHMV Dataset is loaded, the Map explorer adds `Hoogte & +reliëf`. The active 5 m analysis raster is rendered as a colour-relief +MapLibre image over the ordinary OpenStreetMap context. A drawn rectangle or +the exact selected Area returns measured DTM/DSM height, relief and slope with +TAW and the 2013-2015 source period visible. Raster cells are not presented as +objects. GeoJSON export is disabled for this raster-only result, and the UI +explicitly states that water depth and volume cannot be derived from DHMV. + ## Useful repository scripts - `bash scripts/frontend_install.sh` diff --git a/frontend/src/App.tsx b/frontend/src/App.tsx index 27caea52..5ace1680 100644 --- a/frontend/src/App.tsx +++ b/frontend/src/App.tsx @@ -182,7 +182,13 @@ function App(): JSX.Element { ) const availableVectorDatasets = useMemo(() => datasets.filter((item) => isVectorDatasetType(item.dataset_type)), [datasets]) const availableMapDatasets = useMemo( - () => datasets.filter((dataset) => isVectorDatasetType(dataset.dataset_type) && dataset.status === 'ready'), + () => { + const vectors = datasets.filter((dataset) => isVectorDatasetType(dataset.dataset_type) && dataset.status === 'ready') + const terrain = datasets.filter( + (dataset) => dataset.dataset_type === 'raster' && dataset.source_name === 'digitaal_vlaanderen_dhmv' && dataset.status === 'ready', + ) + return [...vectors, ...terrain] + }, [datasets], ) const referenceDatasets = useMemo( @@ -1023,7 +1029,7 @@ function App(): JSX.Element { orthophotoResult={mapOrthophotoAnalysis.lastResult} orthophotoImageUrl={mapOrthophotoAnalysis.imageUrl} availableMapDatasets={availableMapDatasets} - selectedMapDatasetId={selectedDataset && isVectorDatasetType(selectedDataset.dataset_type) ? selectedDataset.id : ''} + selectedMapDatasetId={selectedDataset && availableMapDatasets.some((dataset) => dataset.id === selectedDataset.id) ? selectedDataset.id : ''} selectedFeature={selectedMapFeature} onSelectMapArea={setSelectedMapAreaId} onOpenDatasetInMap={openDatasetInMap} diff --git a/frontend/src/components/datasets/SourceCatalogPanel.tsx b/frontend/src/components/datasets/SourceCatalogPanel.tsx index 457dd782..8ef90257 100644 --- a/frontend/src/components/datasets/SourceCatalogPanel.tsx +++ b/frontend/src/components/datasets/SourceCatalogPanel.tsx @@ -12,6 +12,7 @@ const THEME_LABELS: Record = { nature_value: 'Natuurwaarde', agriculture: 'Landbouw', water: 'Water', + elevation: 'Hoogte en reliëf', roads: 'Wegen en transport', parcels: 'Percelen', } @@ -54,7 +55,7 @@ const AVAILABLE_SOURCES = [ name: 'Digitaal Hoogtemodel Vlaanderen II', owner: 'Digitaal Vlaanderen', coverage: 'LiDAR-opname 2013-2015, DTM/DSM 1 m en 5 m', - value: 'Hoogte, reliëf, helling en afstroming; geen waterdiepte', + value: 'Hoogte, reliëf en helling; afstroming is afgeleid, geen waterdiepte', url: 'https://www.vlaanderen.be/digitaal-vlaanderen/onze-diensten-en-platformen/earth-observation-data-science-eodas/het-digitaal-hoogtemodel/digitaal-hoogtemodel-vlaanderen-ii', }, { @@ -111,6 +112,7 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E const buildingsRegisterDatasets = ready.filter( (dataset) => dataset.source_name === 'digitaal_vlaanderen_buildings_addresses_register', ) + const dhmvDatasets = ready.filter((dataset) => dataset.source_name === 'digitaal_vlaanderen_dhmv') const latestBuildingsRegister = [...buildingsRegisterDatasets].sort( (left, right) => new Date(right.observed_at ?? 0).getTime() - new Date(left.observed_at ?? 0).getTime(), )[0] @@ -122,6 +124,7 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E if (source.key === 'bwk') return bwkDatasets.length === 0 if (source.key === 'agriculture') return agricultureDatasets.length === 0 if (source.key === 'buildings_register') return buildingsRegisterDatasets.length === 0 + if (source.key === 'elevation') return dhmvDatasets.length === 0 return true }) const themes = Object.keys(THEME_LABELS).map((theme) => { @@ -181,7 +184,7 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E ))} - {waterinfoDatasets.length > 0 || historicalOrthophotos.length > 0 || bwkDatasets.length > 0 || agricultureDatasets.length > 0 || buildingsRegisterDatasets.length > 0 ? ( + {waterinfoDatasets.length > 0 || historicalOrthophotos.length > 0 || bwkDatasets.length > 0 || agricultureDatasets.length > 0 || buildingsRegisterDatasets.length > 0 || dhmvDatasets.length > 0 ? (
{waterinfoDatasets.length > 0 ? (
@@ -222,6 +225,13 @@ export function SourceCatalogPanel({ datasets }: SourceCatalogPanelProps): JSX.E

Registerstatus en geaggregeerde koppelingen voor Mol; adreslabels en persoonsgegevens worden niet in de kaartlaag getoond.

) : null} + {dhmvDatasets.length > 0 ? ( +
+ Digitaal Hoogtemodel Vlaanderen II + {dhmvDatasets.length} rasterproducten · DTM/DSM · analyse op 5 m +

Hoogte in TAW, reliëf en helling uit opnameperiode 2013-2015. Geen waterdiepte of watervolume.

+
+ ) : null}
) : null} diff --git a/frontend/src/components/map/MapWorkspace.tsx b/frontend/src/components/map/MapWorkspace.tsx index 5168192c..32868433 100644 --- a/frontend/src/components/map/MapWorkspace.tsx +++ b/frontend/src/components/map/MapWorkspace.tsx @@ -6,6 +6,7 @@ import { useMapThemeSelectionInsights } from '../../hooks/useMapThemeSelectionIn import { useTemporalComparison } from '../../hooks/useTemporalComparison' import { getDatasetDisplayName, getDatasetSourceDisplayName } from '../../lib/datasetDisplay' import { TemporalTrendChart } from './TemporalTrendChart' +import { terrainImageUrl } from '../../lib/terrainImage' const DEFAULT_SELECTED_FEATURE_FILENAME = 'selected-feature.geojson' const DEFAULT_AREA_SELECTION_FILENAME = 'area-selection.geojson' @@ -13,7 +14,7 @@ const EMPTY_TEMPORAL_SERIES: DatasetCreateResponse[] = [] const MOL_PROJECT_NAME = 'Mol Municipality Workbench' const KEMPEN_PROJECT_NAME = 'Kempen Regional Workbench' -type DataThemeId = 'buildings' | 'population' | 'forest' | 'nature_value' | 'agriculture' | 'water' | 'roads' | 'parcels' +type DataThemeId = 'buildings' | 'population' | 'forest' | 'nature_value' | 'agriculture' | 'water' | 'elevation' | 'roads' | 'parcels' interface DataTheme { id: DataThemeId @@ -72,6 +73,13 @@ const DATA_THEMES: DataTheme[] = [ description: 'Waterlopen, grachten, kanalen en wateroppervlakken.', tokens: ['waterways', 'waterway', 'water', 'hydro', 'river', 'stream', 'canal', 'waterloop'], }, + { + id: 'elevation', + label: 'Hoogte & reliëf', + shortLabel: 'Hoogte', + description: 'Maaiveld- of oppervlaktehoogte, reliëf en helling uit DHMV II.', + tokens: ['dhmv', 'elevation', 'height', 'hoogte', 'terrain', 'surface', 'dtm', 'dsm', 'reliëf'], + }, { id: 'roads', label: 'Wegen', @@ -95,10 +103,19 @@ const DATA_THEME_MAP_STYLES: Record nature_value: { fill: '#9a4f64', line: '#74364a' }, agriculture: { fill: '#7b8f32', line: '#53671d' }, water: { fill: '#2676a8', line: '#155b85' }, + elevation: { fill: '#a57a4b', line: '#315f59' }, roads: { fill: '#6b7280', line: '#4b5563' }, parcels: { fill: '#a7792f', line: '#7d571f' }, } +function datasetAvailabilityLabel(dataset: DatasetCreateResponse): string { + if (dataset.dataset_type === 'raster' && dataset.source_name === 'digitaal_vlaanderen_dhmv') { + const resolution = Number(dataset.source_metadata?.['analysis_resolution_m']) + return `${Number.isFinite(resolution) ? `${resolution.toLocaleString('nl-BE')} m` : 'Raster'} hoogtegrid beschikbaar` + } + return `${(dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0).toLocaleString('nl-BE')} objecten beschikbaar` +} + function datasetSearchText(dataset: DatasetCreateResponse): string { return [ dataset.name, @@ -144,6 +161,8 @@ function pickThemeDataset( (dataset.source_name === 'inbo_bwk_natura2000' ? 95_000 : 0) + (dataset.source_name === 'agentschap_landbouw_zeevisserij_agricultural_parcels' ? 98_000 : 0) + (dataset.source_name === 'digitaal_vlaanderen_buildings_addresses_register' ? 5_000_000 : 0) + + (dataset.source_name === 'digitaal_vlaanderen_dhmv' ? 5_000_000 : 0) + + (dataset.source_metadata?.['product_key'] === 'dtm_1m' ? 1_000_000 : 0) + (dataset.dataset_role === 'reference' ? 10_000 : 0) + (dataset.observed_at ? new Date(dataset.observed_at).getTime() / 100_000_000 : 0) + (dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0) @@ -200,6 +219,14 @@ function formatObservationDate(value: string | null | undefined): string { return new Intl.DateTimeFormat('nl-BE', { year: 'numeric', month: 'short', day: 'numeric' }).format(new Date(value)) } +function formatDatasetObservation(dataset: DatasetCreateResponse): string { + const period = dataset.source_metadata?.['acquisition_period'] + if (typeof period === 'string' && period.trim()) { + return `opnameperiode ${period}` + } + return formatObservationDate(dataset.observed_at) +} + function operationalScopeProjectLabel(project: ProjectRead): string { if (project.name === MOL_PROJECT_NAME) { return 'Mol' @@ -668,6 +695,18 @@ export function MapWorkspace({ } : null const activeThemeDataset = themeDatasetMap[activeTheme.id] + const terrainBounds = activeThemeDataset?.source_name === 'digitaal_vlaanderen_dhmv' + ? activeThemeDataset.source_metadata?.['bbox_epsg4326'] + : null + const terrainImageOverlay = activeTheme.id === 'elevation' && activeThemeDataset && selectedProjectId && Array.isArray(terrainBounds) && terrainBounds.length === 4 + ? { + url: terrainImageUrl(selectedProjectId, activeThemeDataset.id), + bbox: terrainBounds.map(Number) as [number, number, number, number], + label: getDatasetDisplayName(activeThemeDataset), + opacity: 0.82, + } + : null + const activeImageOverlay = terrainImageOverlay ?? orthophotoImageOverlay const activeScopeProject = projects.find((project) => project.id === selectedProjectId) ?? null const activeScopeLabel = activeScopeProject ? operationalScopeProjectLabel(activeScopeProject) : 'Werkgebied' const municipalityAreaCount = areas.filter((area) => /^Gemeente\s/i.test(area.name)).length @@ -709,11 +748,14 @@ export function MapWorkspace({ const activeSupportingMetrics = (activeSelectionResult?.summary?.metrics ?? []).filter( (metric) => metric.metric_key !== activeSelectionResult?.summary?.primary_metric_key, ) - const activeSecondaryMetric = selectedAreaSquareMetres && selectedAreaSquareMetres > 0 - ? activeMetricUnit === 'ha' + const terrainReliefMetric = activeSupportingMetrics.find((metric) => metric.metric_key === 'relief_m') + const activeSecondaryMetric = activeMetricUnit === 'm TAW' + ? terrainReliefMetric ? `${terrainReliefMetric.metric_value.toLocaleString('nl-BE', { maximumFractionDigits: 2 })} m reliëf` : null + : selectedAreaSquareMetres && selectedAreaSquareMetres > 0 + ? activeMetricUnit === 'ha' ? `${((activeMetricValue * 10_000) / selectedAreaSquareMetres * 100).toLocaleString('nl-BE', { maximumFractionDigits: 1 })}% dekking` : `${(activeMetricValue / (selectedAreaSquareMetres / 1_000_000)).toLocaleString('nl-BE', { maximumFractionDigits: 1 })} ${activeMetricUnit} / km2` - : null + : null const selectedResultProperties = useMemo(() => { const keys = new Map>() for (const feature of activeSelectionResult?.geojson.features ?? []) { @@ -1094,7 +1136,7 @@ export function MapWorkspace({ ? 'Alleen huidige toestand' : 'Bron nog niet ingeladen' : dataset - ? `${(dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0).toLocaleString('nl-BE')} objecten beschikbaar` + ? datasetAvailabilityLabel(dataset) : 'Bron nog niet ingeladen'} @@ -1125,7 +1167,7 @@ export function MapWorkspace({ ? `${activeTemporalSeries.length} officiële meetmomenten · ${formatObservationDate(activeTemporalSeries[0].observed_at)} tot ${formatObservationDate(activeTemporalSeries[activeTemporalSeries.length - 1].observed_at)}` : 'Voor dit thema is nog geen tweede officieel meetmoment beschikbaar.' : activeThemeDataset - ? `${getDatasetSourceDisplayName(activeThemeDataset)} · ${formatObservationDate(activeThemeDataset.observed_at)}` + ? `${getDatasetSourceDisplayName(activeThemeDataset)} · ${formatDatasetObservation(activeThemeDataset)}` : activeTheme.description} @@ -1225,7 +1267,7 @@ export function MapWorkspace({ areaData={areaFeatureCollection} selectedFeature={selectedFeature} selectionData={analysisMode === 'current' ? mapSelectionResult?.geojson ?? null : null} - imageOverlay={orthophotoImageOverlay} + imageOverlay={activeImageOverlay} selectionBbox={mapSelectionBbox} bboxSelectionMode={bboxSelectionMode} visible={mapLayerVisible} @@ -1241,7 +1283,7 @@ export function MapWorkspace({ />
Werkgebied - {orthophotoImageOverlay ? {orthophotoImageOverlay.label} : null} + {activeImageOverlay ? {activeImageOverlay.label} : null} {analysisOverlayActive ? ( <> AI-kandidaten @@ -1444,7 +1486,7 @@ export function MapWorkspace({ {activeSelectionResult ? resultMetricLabel(activeSelectionResult) : 'Geen resultaat'}
- {activeMetricUnit === 'ha' ? 'Aandeel selectie' : 'Dichtheid'} + {activeMetricUnit === 'ha' ? 'Aandeel selectie' : activeMetricUnit === 'm TAW' ? 'Reliëf' : 'Dichtheid'} {activeSecondaryMetric ?? (selectedDensity === null ? 'n.v.t.' : `${selectedDensity.toLocaleString('nl-BE', { maximumFractionDigits: 1 })} / km2`)}
@@ -1517,7 +1559,7 @@ export function MapWorkspace({ {analysisMode === 'current' ? (
- +
) : null} diff --git a/frontend/src/hooks/useMapSelectionExtract.ts b/frontend/src/hooks/useMapSelectionExtract.ts index 23cf6110..56659f4b 100644 --- a/frontend/src/hooks/useMapSelectionExtract.ts +++ b/frontend/src/hooks/useMapSelectionExtract.ts @@ -2,6 +2,7 @@ import { useEffect, useRef, useState } from 'react' import { datasetsApi } from '../services/api' import { formatError } from '../lib/formatError' import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types' +import { terrainSelectionToMapSelection } from '../lib/terrainSelection' interface MapSelectionExtractOptions { selectedProjectId: string | null @@ -38,8 +39,9 @@ export function useMapSelectionExtract({ setMapSelectionError('Open a vector dataset before extracting a map area.') return null } - if (!isVectorDatasetType(selectedDataset.dataset_type)) { - setMapSelectionError('Area extraction requires an active vector dataset.') + const terrainDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'digitaal_vlaanderen_dhmv' + if (!isVectorDatasetType(selectedDataset.dataset_type) && !terrainDataset) { + setMapSelectionError('Gebiedsanalyse ondersteunt een vectorlaag of een beheerd DHMV-hoogtemodel.') return null } @@ -49,11 +51,16 @@ export function useMapSelectionExtract({ setMapSelectionError(null) setMapSelectionBbox(bbox) try { - const response = await datasetsApi.selectVectorFeatures(selectedProjectId, selectedDataset.id, { - bbox: { ...bbox, crs: 'EPSG:4326' }, - area_id: areaId, - limit: 1000, - }) + const response = terrainDataset + ? terrainSelectionToMapSelection(await datasetsApi.selectTerrain(selectedProjectId, selectedDataset.id, { + bbox: { ...bbox, crs: 'EPSG:4326' }, + area_id: areaId, + })) + : await datasetsApi.selectVectorFeatures(selectedProjectId, selectedDataset.id, { + bbox: { ...bbox, crs: 'EPSG:4326' }, + area_id: areaId, + limit: 1000, + }) if (requestSequence.current !== sequence) { return null } diff --git a/frontend/src/hooks/useMapThemeSelectionInsights.ts b/frontend/src/hooks/useMapThemeSelectionInsights.ts index 09b6460f..6e083e87 100644 --- a/frontend/src/hooks/useMapThemeSelectionInsights.ts +++ b/frontend/src/hooks/useMapThemeSelectionInsights.ts @@ -2,6 +2,7 @@ import { useEffect, useRef, useState } from 'react' import { formatError } from '../lib/formatError' import { datasetsApi } from '../services/api/datasets' import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types' +import { terrainSelectionToMapSelection } from '../lib/terrainSelection' export interface MapThemeQuery { themeId: TThemeId @@ -54,11 +55,16 @@ export function useMapThemeSelectionInsights( queries.map(async ({ themeId, dataset }) => ({ themeId, dataset, - result: await datasetsApi.selectVectorFeatures(selectedProjectId, dataset.id, { - bbox, - area_id: areaId, - limit: 1000, - }), + result: dataset.dataset_type === 'raster' && dataset.source_name === 'digitaal_vlaanderen_dhmv' + ? terrainSelectionToMapSelection(await datasetsApi.selectTerrain(selectedProjectId, dataset.id, { + bbox, + area_id: areaId, + })) + : await datasetsApi.selectVectorFeatures(selectedProjectId, dataset.id, { + bbox, + area_id: areaId, + limit: 1000, + }), })), ) const successful = settled.flatMap((item) => (item.status === 'fulfilled' ? [item.value] : [])) diff --git a/frontend/src/lib/datasetDisplay.ts b/frontend/src/lib/datasetDisplay.ts index 91d3ab27..cb1cbe51 100644 --- a/frontend/src/lib/datasetDisplay.ts +++ b/frontend/src/lib/datasetDisplay.ts @@ -9,6 +9,7 @@ const DATASET_LABEL_BY_LAYER: Record = { forest: 'Bos en groen', nature_value: 'Natuurwaarde', agriculture: 'Landbouwgebruikspercelen', + elevation: 'Hoogte en reliëf', building_registry: 'Gebouwenregister', regional_boundary: 'Grens vervoerregio Kempen', municipality_boundaries: 'Gemeentegrenzen Kempen', @@ -19,6 +20,7 @@ const DATASET_SOURCE_LABELS: Record = { agentschap_landbouw_zeevisserij_agricultural_parcels: 'Agentschap Landbouw en Zeevisserij', digitaal_vlaanderen_orthophoto: 'Digitaal Vlaanderen', digitaal_vlaanderen_buildings_addresses_register: 'Digitaal Vlaanderen', + digitaal_vlaanderen_dhmv: 'Digitaal Vlaanderen', grb: 'GRB', historical_landuse: 'Digitaal Vlaanderen', inbo_bwk_natura2000: 'INBO', @@ -33,6 +35,10 @@ export function getDatasetSourceDisplayName(dataset: DatasetCreateResponse): str } export function getDatasetDisplayName(dataset: DatasetCreateResponse): string { + if (dataset.source_name === 'digitaal_vlaanderen_dhmv') { + const productName = dataset.source_metadata?.['product_display_name'] + return typeof productName === 'string' && productName.trim() ? productName : 'DHMV II hoogtemodel' + } const layer = (dataset.reference_layer_name ?? dataset.source_metadata?.layer_name ?? dataset.source_metadata?.layer_type ?? '') .toString() .toLowerCase() diff --git a/frontend/src/lib/terrainImage.ts b/frontend/src/lib/terrainImage.ts new file mode 100644 index 00000000..14b7cec0 --- /dev/null +++ b/frontend/src/lib/terrainImage.ts @@ -0,0 +1,3 @@ +export function terrainImageUrl(projectId: string, datasetId: string): string { + return `/api/v1/projects/${projectId}/datasets/${datasetId}/raster/terrain/image` +} diff --git a/frontend/src/lib/terrainSelection.ts b/frontend/src/lib/terrainSelection.ts new file mode 100644 index 00000000..7a4188e8 --- /dev/null +++ b/frontend/src/lib/terrainSelection.ts @@ -0,0 +1,23 @@ +import type { TerrainSelectionResponse, VectorSelectionResponse } from '../types' + +export function terrainSelectionToMapSelection(result: TerrainSelectionResponse): VectorSelectionResponse { + return { + selection_bbox: result.selection_bbox, + selection_area_id: result.selection_area_id, + feature_count: 0, + total_feature_count: 0, + limit: 0, + truncated: false, + geojson: { type: 'FeatureCollection', features: [] }, + summary: { + ...result.summary, + feature_count: result.sample_count, + is_estimate: false, + warning: result.limitation_message, + metrics: result.summary.metrics.map((metric) => ({ + ...metric, + is_estimate: false, + })), + }, + } +} diff --git a/frontend/src/services/api/datasets.ts b/frontend/src/services/api/datasets.ts index 8a8bf13e..fb1c0cab 100644 --- a/frontend/src/services/api/datasets.ts +++ b/frontend/src/services/api/datasets.ts @@ -18,6 +18,9 @@ import type { RasterNdbiRequest, OrthophotoAcquireRequest, OrthophotoProductRead, + DhmvAcquireRequest, + DhmvProductRead, + TerrainSelectionResponse, } from '../../types' const DATASET_PAGE_SIZE = 200 @@ -115,6 +118,16 @@ export const datasetsApi = { apiGet<{ items: OrthophotoProductRead[]; total: number }>(`/api/v1/projects/${projectId}/datasets/orthophoto/products`), orthophotoImageUrl: (projectId: string, datasetId: string): string => `/api/v1/projects/${projectId}/datasets/${datasetId}/raster/image`, + acquireDhmv: (projectId: string, payload: DhmvAcquireRequest): Promise => + apiPost(`/api/v1/projects/${projectId}/datasets/dhmv/acquire`, payload), + listDhmvProducts: (projectId: string): Promise<{ items: DhmvProductRead[]; total: number }> => + apiGet<{ items: DhmvProductRead[]; total: number }>(`/api/v1/projects/${projectId}/datasets/dhmv/products`), + selectTerrain: ( + projectId: string, + datasetId: string, + payload: { bbox: VectorSelectionRequest['bbox']; area_id?: string }, + ): Promise => + apiPost(`/api/v1/projects/${projectId}/datasets/${datasetId}/raster/terrain/select`, payload), refreshMetadata: (projectId: string, datasetId: string): Promise => apiPost(`/api/v1/projects/${projectId}/datasets/${datasetId}/metadata/refresh`, {}), inspectRaster: (projectId: string, datasetId: string): Promise => diff --git a/frontend/src/styles/app.css b/frontend/src/styles/app.css index aaecc6da..efe1d0f7 100644 --- a/frontend/src/styles/app.css +++ b/frontend/src/styles/app.css @@ -5756,6 +5756,7 @@ section { .geo-theme-symbol-nature_value { background: #9a4f64; } .geo-theme-symbol-agriculture { background: #7b8f32; } .geo-theme-symbol-water { background: #2676a8; } +.geo-theme-symbol-elevation { background: #a57a4b; } .geo-theme-symbol-roads { background: #6b7280; } .geo-theme-symbol-parcels { background: #a7792f; } @@ -5940,6 +5941,11 @@ section { background: rgba(38, 118, 168, 0.24); } +.geo-map-legend .geo-legend-layer-elevation { + border-color: #315f59; + background: rgba(165, 122, 75, 0.26); +} + .geo-map-legend .geo-legend-layer-roads { border-color: #4b5563; background: rgba(107, 114, 128, 0.24); diff --git a/frontend/src/types.ts b/frontend/src/types.ts index a4ee877c..536ce974 100644 --- a/frontend/src/types.ts +++ b/frontend/src/types.ts @@ -336,6 +336,52 @@ export interface OrthophotoAcquisitionResult { limitation_message: string } +export interface DhmvAcquireRequest { + bbox: VectorSelectionBBox + area_id?: string | null + product_key?: 'dtm_1m' | 'dsm_1m' + resolution_m?: number | null + force_refresh?: boolean +} + +export interface DhmvProductRead { + key: 'dtm_1m' | 'dsm_1m' + display_name: string + surface_model: 'terrain' | 'surface' + coverage_id: string + native_resolution_m: number + source_crs: string + vertical_reference: string + acquisition_period: string + catalog_url: string + attribution: string + limitation_message: string +} + +export interface TerrainSelectionResponse { + dataset_id: string + product_key: string + surface_model: 'terrain' | 'surface' + selection_bbox: VectorSelectionBBox + selection_area_id?: string | null + sample_count: number + slope_sample_count: number + coverage_ratio: number + resolution_m: number + vertical_reference: string + summary: { + metric_label: string + metric_value: number + metric_unit: string + aggregation_method: string + primary_metric_key: string + metrics: VectorSelectionMetric[] + } + unsupported_metrics: string[] + limitation_message: string + generated_at: string +} + export interface MapImageOverlay { url: string bbox: [number, number, number, number] diff --git a/scripts/README.md b/scripts/README.md index 556fbe72..3bc4ac63 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -1513,6 +1513,31 @@ Only aggregate unit/address counts enter the queryable building layer. Review manifest before accepting a broader import. Raw address response pages are operator evidence and must not be published. +## Mol DHMV terrain rasters + +Acquire and validate the official DHMV II DTM and DSM for the exact persisted +Mol Area: + +```bash +docker exec geointel python /app/scripts/provision_mol_dhmv.py +``` + +The operator resolves project and Area through the API, derives the bounded +EPSG:4326 request rectangle and calls the canonical DHMV endpoints. The backend +requests the fixed official WCS coverages, extracts multipart GeoTIFF, clips to +the exact Area, validates EPSG:31370/resolution/nodata/valid cells and stores +through DatasetService. It then runs a full-Area terrain metric smoke. + +Useful safe overrides: + +```bash +docker exec geointel python /app/scripts/provision_mol_dhmv.py --products dtm_1m +docker exec geointel python /app/scripts/provision_mol_dhmv.py --resolution-m 5 --force +``` + +Do not use DHMV output as water depth or water volume. The command fails when +the API no longer reports those metrics as explicitly unsupported. + ## Tower deployment Push the local branch to Gitea, then rebuild the Unraid/Tower Docker runtime: diff --git a/scripts/provision_mol_dhmv.py b/scripts/provision_mol_dhmv.py new file mode 100644 index 00000000..672554d5 --- /dev/null +++ b/scripts/provision_mol_dhmv.py @@ -0,0 +1,157 @@ +"""Provision governed DHMV II terrain/surface rasters for the persisted Mol Area. + +The operator calls the canonical GeoIntel DHMV acquisition API. The backend +performs the bounded official WCS request, exact Area clipping, validation, +checksum storage and Dataset/Job persistence. No raster rows are written +directly by this script. +""" + +from __future__ import annotations + +import argparse +import json +import os +from typing import Any, Iterable + +import requests + + +DEFAULT_API_URL = "http://127.0.0.1:8000" +DEFAULT_PROJECT_NAME = "Kempen Regional Workbench" +DEFAULT_AREA_FRAGMENT = "Gemeente Mol" +PRODUCTS = ("dtm_1m", "dsm_1m") + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Provision official DHMV II DTM/DSM rasters for Mol.") + parser.add_argument("--base-url", default=os.environ.get("GEOINTEL_INTERNAL_API_URL", DEFAULT_API_URL)) + parser.add_argument("--project-name", default=DEFAULT_PROJECT_NAME) + parser.add_argument("--area-name", default=DEFAULT_AREA_FRAGMENT) + parser.add_argument("--products", default=",".join(PRODUCTS), help="Comma-separated governed product keys.") + parser.add_argument("--resolution-m", type=float, default=5.0) + parser.add_argument("--timeout", type=int, default=1800) + parser.add_argument("--force", action="store_true") + return parser.parse_args() + + +def unwrap(response: requests.Response) -> Any: + response.raise_for_status() + payload = response.json() + if not isinstance(payload, dict) or "data" not in payload: + raise RuntimeError(f"Non-canonical API response from {response.url}") + return payload["data"] + + +def coordinates(geometry: dict[str, Any]) -> Iterable[tuple[float, float]]: + def walk(value: Any): + if isinstance(value, list) and len(value) >= 2 and all(isinstance(item, (int, float)) for item in value[:2]): + yield float(value[0]), float(value[1]) + return + if isinstance(value, list): + for child in value: + yield from walk(child) + + yield from walk(geometry.get("coordinates", [])) + + +def geometry_bbox(geometry: dict[str, Any]) -> dict[str, float | str]: + points = list(coordinates(geometry)) + if not points: + raise RuntimeError("Persisted Area geometry contains no coordinates") + xs = [point[0] for point in points] + ys = [point[1] for point in points] + return { + "min_x": min(xs), + "min_y": min(ys), + "max_x": max(xs), + "max_y": max(ys), + "crs": "EPSG:4326", + } + + +def main() -> int: + args = parse_args() + base_url = args.base_url.rstrip("/") + session = requests.Session() + session.headers.update({"User-Agent": "GeoIntel-DHMV-Operator/1.0"}) + + projects = unwrap(session.get(f"{base_url}/api/v1/projects", params={"limit": 200, "offset": 0}, timeout=60))["items"] + project = next((item for item in projects if item["name"] == args.project_name), None) + if project is None: + raise RuntimeError(f"Project {args.project_name!r} was not found") + + areas = unwrap( + session.get( + f"{base_url}/api/v1/projects/{project['id']}/areas", + params={"limit": 200, "offset": 0}, + timeout=60, + ) + )["items"] + fragment = args.area_name.casefold() + area = next((item for item in areas if fragment in item["name"].casefold()), None) + if area is None: + raise RuntimeError(f"Area containing {args.area_name!r} was not found") + bbox = geometry_bbox(area["geometry"]) + + requested_products = [item.strip() for item in args.products.split(",") if item.strip()] + invalid = sorted(set(requested_products) - set(PRODUCTS)) + if invalid: + raise RuntimeError(f"Unsupported DHMV product keys: {', '.join(invalid)}") + + results = [] + for product_key in requested_products: + job = unwrap( + session.post( + f"{base_url}/api/v1/projects/{project['id']}/datasets/dhmv/acquire", + json={ + "bbox": bbox, + "area_id": area["id"], + "product_key": product_key, + "resolution_m": args.resolution_m, + "force_refresh": args.force, + }, + timeout=args.timeout, + ) + ) + if job.get("status") != "success" or not job.get("output_dataset_id"): + raise RuntimeError(f"DHMV acquisition failed for {product_key}: {job.get('error_message') or job}") + analysis = unwrap( + session.post( + f"{base_url}/api/v1/projects/{project['id']}/datasets/{job['output_dataset_id']}/raster/terrain/select", + json={"bbox": bbox, "area_id": area["id"]}, + timeout=args.timeout, + ) + ) + if sorted(analysis.get("unsupported_metrics", [])) != ["water_depth_m", "water_volume_m3"]: + raise RuntimeError("DHMV terrain contract must explicitly keep water depth and volume unavailable") + results.append( + { + "product_key": product_key, + "dataset_id": job["output_dataset_id"], + "reused": bool((job.get("result_json") or {}).get("reused")), + "resolution_m": analysis["resolution_m"], + "sample_count": analysis["sample_count"], + "coverage_ratio": analysis["coverage_ratio"], + "metrics": analysis["summary"]["metrics"], + } + ) + + print( + json.dumps( + { + "status": "ok", + "project_id": project["id"], + "area_id": area["id"], + "area_name": area["name"], + "bbox": bbox, + "products": results, + }, + ensure_ascii=False, + indent=2, + ) + ) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/run_readiness_check.sh b/scripts/run_readiness_check.sh index 79f0fd5c..1a0c1f61 100755 --- a/scripts/run_readiness_check.sh +++ b/scripts/run_readiness_check.sh @@ -51,6 +51,7 @@ ${PYTHON_BIN} -m py_compile scripts/provision_waterinfo_station_history.py ${PYTHON_BIN} -m py_compile scripts/provision_mol_bwk_natura2000.py ${PYTHON_BIN} -m py_compile scripts/provision_agricultural_parcel_history.py ${PYTHON_BIN} -m py_compile scripts/provision_buildings_addresses_register.py +${PYTHON_BIN} -m py_compile scripts/provision_mol_dhmv.py ${PYTHON_BIN} -m py_compile scripts/provision_regional_timeseries.py ${PYTHON_BIN} -m py_compile scripts/geographic_scopes.py ${PYTHON_BIN} -m py_compile scripts/provision_geographic_scope.py