diff --git a/.env.example b/.env.example index fa8d40be..276e0181 100644 --- a/.env.example +++ b/.env.example @@ -28,6 +28,13 @@ FLOOD_HAZARD_MAX_SIDE_M=20000 FLOOD_HAZARD_MAX_PIXELS=12000000 FLOOD_HAZARD_TIMEOUT_SECONDS=300 FLOOD_HAZARD_MAX_RESPONSE_MB=160 +THEMATIC_RASTER_ENABLED=true +THEMATIC_RASTER_WCS_URL=https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs +THEMATIC_RASTER_MIN_SIDE_M=100 +THEMATIC_RASTER_MAX_SIDE_M=20000 +THEMATIC_RASTER_MAX_PIXELS=12000000 +THEMATIC_RASTER_TIMEOUT_SECONDS=300 +THEMATIC_RASTER_MAX_RESPONSE_MB=160 YOLO_ENABLED=false YOLO_MODELS_DIR=/app/models YOLO_MODEL_PATH= diff --git a/CHANGELOG.md b/CHANGELOG.md index 871c5429..fbdf810f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,30 @@ # Changelog +## Sprint 213-214 Cross-domain area profile (2026-07-16) + +- Implemented one allowlisted MercatorNet WCS registry for official Flemish + space occupation 2025, open space 2022, population density 2019, public + transport node value 2022 and total service level 2022 rasters. +- Added bounded tiled acquisition, exact Area clipping in EPSG:31370, raster + value/unit validation, checksummed Dataset/DatasetVersion provenance and + source-correct selection metrics without accepting arbitrary service URLs or + coverage identifiers. +- Added MapLibre image overlays, legends and current-state selection for all + five products. Raster cell values are presented as hectares, an explicitly + estimated population total/density or source scores, never as object counts. +- Grounded local Ollama answers in the persisted thematic measurements and + retained unsupported-current-count, live-timetable and causal limitations. +- Added the official DOV digital soil map as an explicit Mol operator. It + paginates all bounded `bodemkaart:bodemtypes` features, stores checksummed raw + evidence, clips exactly in Lambert 72 and imports through DatasetService. +- Added a Soil map theme with mapped hectares and inspectable soil type, + texture and drainage fields. The 1949-1971 survey period and 1:20,000 scale + remain visible; current drainage is never inferred. +- Added focused acquisition, GIS, analysis, API, AI-context, operator, + packaging and frontend-contract tests. No migration or direct database write + was introduced. + ## Sprint 212 Platform-wide official source portfolio (2026-07-16) - Rebalanced the source strategy across six user-facing domains: space and diff --git a/backend/README.md b/backend/README.md index 25ef7824..dc28c5df 100644 --- a/backend/README.md +++ b/backend/README.md @@ -1289,6 +1289,43 @@ Settings: `FLOOD_HAZARD_ENABLED`, `FLOOD_HAZARD_WCS_URL`, `FLOOD_HAZARD_MAX_SIDE_M`, `FLOOD_HAZARD_MAX_PIXELS`, `FLOOD_HAZARD_TIMEOUT_SECONDS` and `FLOOD_HAZARD_MAX_RESPONSE_MB`. +## Cross-domain thematic rasters and DOV soil + +The governed thematic registry exposes five fixed MercatorNet products through +`GET .../datasets/thematic-raster/products`. Acquisition uses +`POST .../datasets/thematic-raster/acquire`; selection and PNG rendering use +`POST .../raster/thematic/select` and `GET .../raster/thematic/image`. + +Provision every product for the exact persisted Mol Area: + +```bash +docker exec geointel python /app/scripts/provision_thematic_rasters.py +``` + +Inspect the complete 28-municipality matrix without writes, then run it after +the Mol source/runtime gate passes: + +```bash +docker exec geointel python /app/scripts/provision_thematic_rasters.py \ + --project-name "Kempen Regional Workbench" --all-municipalities --dry-run +``` + +Settings: `THEMATIC_RASTER_ENABLED`, `THEMATIC_RASTER_WCS_URL`, +`THEMATIC_RASTER_MIN_SIDE_M`, `THEMATIC_RASTER_MAX_SIDE_M`, +`THEMATIC_RASTER_MAX_PIXELS`, `THEMATIC_RASTER_TIMEOUT_SECONDS` and +`THEMATIC_RASTER_MAX_RESPONSE_MB`. + +Provision the official DOV soil polygons for Mol through the existing vector +upload path: + +```bash +docker exec geointel python /app/scripts/provision_mol_soil_map.py +``` + +Use `--fetch-only` to retain and validate source evidence without importing. +The operator never writes directly to PostGIS. Soil drainage and related map +classes represent the 1949-1971 survey and are not current observations. + ## 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 e52bb70f..830c53f9 100644 --- a/backend/app/api/routes/datasets.py +++ b/backend/app/api/routes/datasets.py @@ -26,6 +26,8 @@ from app.schemas import ( TerrainSelectionRequest, FloodHazardAcquireRequest, FloodHazardSelectionRequest, + ThematicRasterAcquireRequest, + ThematicRasterSelectionRequest, VectorBBoxResponse, VectorBufferRequest, VectorClipRequest, @@ -48,6 +50,8 @@ from app.services.dhmv_acquisition_service import DhmvAcquisitionService from app.services.terrain_analysis_service import TerrainAnalysisService from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService +from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService +from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService from app.utils.response import envelope router = APIRouter(prefix="/projects/{project_id}", tags=["datasets"]) @@ -207,6 +211,30 @@ def list_flood_hazard_products(project_id: UUID, db: Session = Depends(get_db)): return envelope({"items": items, "total": len(items)}) +@router.post("/datasets/thematic-raster/acquire", response_model=dict) +def acquire_bounded_thematic_raster( + project_id: UUID, + payload: ThematicRasterAcquireRequest, + db: Session = Depends(get_db), +): + job = JobService.run_sync_job( + db=db, + project_id=project_id, + job_type="raster.thematic.acquire", + parameters=payload.model_dump(mode="json"), + operation=lambda: ThematicRasterAcquisitionService.acquire(db, project_id, payload), + ) + return envelope(job) + + +@router.get("/datasets/thematic-raster/products", response_model=dict) +def list_thematic_raster_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 = ThematicRasterAcquisitionService.list_products() + return envelope({"items": items, "total": len(items)}) + + @router.get("/datasets", response_model=dict) def list_datasets( project_id: UUID, @@ -560,6 +588,30 @@ def raster_flood_hazard_image( ) +@router.post("/datasets/{dataset_id}/raster/thematic/select", response_model=dict) +def raster_thematic_selection( + project_id: UUID, + dataset_id: UUID, + payload: ThematicRasterSelectionRequest, + db: Session = Depends(get_db), +): + return envelope(ThematicRasterAnalysisService.analyze(db, project_id, dataset_id, payload)) + + +@router.get("/datasets/{dataset_id}/raster/thematic/image") +def raster_thematic_image( + project_id: UUID, + dataset_id: UUID, + db: Session = Depends(get_db), +): + content = ThematicRasterAnalysisService.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 1e57c526..24e258d4 100644 --- a/backend/app/core/config.py +++ b/backend/app/core/config.py @@ -53,6 +53,16 @@ class Settings(BaseSettings): flood_hazard_max_pixels: int = Field(default=12_000_000, ge=1, validation_alias="FLOOD_HAZARD_MAX_PIXELS") flood_hazard_timeout_seconds: int = Field(default=300, ge=1, validation_alias="FLOOD_HAZARD_TIMEOUT_SECONDS") flood_hazard_max_response_mb: int = Field(default=160, ge=1, validation_alias="FLOOD_HAZARD_MAX_RESPONSE_MB") + thematic_raster_enabled: bool = Field(default=True, validation_alias="THEMATIC_RASTER_ENABLED") + thematic_raster_wcs_url: str = Field( + default="https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs", + validation_alias="THEMATIC_RASTER_WCS_URL", + ) + thematic_raster_min_side_m: float = Field(default=100.0, gt=0, validation_alias="THEMATIC_RASTER_MIN_SIDE_M") + thematic_raster_max_side_m: float = Field(default=20_000.0, gt=0, validation_alias="THEMATIC_RASTER_MAX_SIDE_M") + thematic_raster_max_pixels: int = Field(default=12_000_000, ge=1, validation_alias="THEMATIC_RASTER_MAX_PIXELS") + thematic_raster_timeout_seconds: int = Field(default=300, ge=1, validation_alias="THEMATIC_RASTER_TIMEOUT_SECONDS") + thematic_raster_max_response_mb: int = Field(default=160, ge=1, validation_alias="THEMATIC_RASTER_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 2d2fb04f..043faafe 100644 --- a/backend/app/schemas/__init__.py +++ b/backend/app/schemas/__init__.py @@ -51,6 +51,15 @@ from .flood_hazard import ( FloodHazardSelectionResponse, FloodHazardSelectionSummary, ) +from .thematic_raster import ( + ThematicRasterAcquireRequest, + ThematicRasterAcquisitionResult, + ThematicRasterMetric, + ThematicRasterProductRead, + ThematicRasterSelectionRequest, + ThematicRasterSelectionResponse, + ThematicRasterSelectionSummary, +) from .external import ( ExternalFetchRequest, ExternalFetchResponse, @@ -167,6 +176,13 @@ __all__ = [ "FloodHazardSelectionRequest", "FloodHazardSelectionResponse", "FloodHazardSelectionSummary", + "ThematicRasterAcquireRequest", + "ThematicRasterAcquisitionResult", + "ThematicRasterMetric", + "ThematicRasterProductRead", + "ThematicRasterSelectionRequest", + "ThematicRasterSelectionResponse", + "ThematicRasterSelectionSummary", "VectorBBoxResponse", "VectorClipRequest", "VectorBufferRequest", diff --git a/backend/app/schemas/thematic_raster.py b/backend/app/schemas/thematic_raster.py new file mode 100644 index 00000000..fddb2218 --- /dev/null +++ b/backend/app/schemas/thematic_raster.py @@ -0,0 +1,96 @@ +from __future__ import annotations + +from uuid import UUID + +from pydantic import BaseModel + +from .operations import VectorSelectionBBox + + +class ThematicRasterAcquireRequest(BaseModel): + bbox: VectorSelectionBBox + area_id: UUID | None = None + product_key: str + force_refresh: bool = False + + +class ThematicRasterProductRead(BaseModel): + key: str + display_name: str + theme: str + metric_kind: str + coverage_id: str + native_resolution_m: float + source_crs: str + source_value_unit: str + observation_year: int + source_version: str + catalog_url: str + attribution: str + license_note: str + legend_min_label: str + legend_max_label: str + limitation_message: str + + +class ThematicRasterAcquisitionResult(BaseModel): + output_dataset_id: UUID + reused: bool + provider: str + product_key: str + display_name: str + theme: str + metric_kind: str + coverage_id: str + resolution_m: float + width: int + height: int + valid_pixel_count: int + bbox_epsg4326: list[float] + bbox_epsg31370: list[float] + observation_year: int + source_value_unit: str + attribution: str + limitation_message: str + + +class ThematicRasterSelectionRequest(BaseModel): + bbox: VectorSelectionBBox + area_id: UUID | None = None + + +class ThematicRasterMetric(BaseModel): + metric_key: str + metric_label: str + metric_value: float + metric_unit: str + aggregation_method: str + derived: bool = True + is_estimate: bool = True + + +class ThematicRasterSelectionSummary(BaseModel): + metric_label: str + metric_value: float + metric_unit: str + aggregation_method: str + primary_metric_key: str + metrics: list[ThematicRasterMetric] + + +class ThematicRasterSelectionResponse(BaseModel): + dataset_id: UUID + product_key: str + theme: str + metric_kind: str + selection_bbox: VectorSelectionBBox + selection_area_id: UUID | None = None + selected_cell_count: int + valid_cell_count: int + coverage_ratio: float + resolution_m: float + observation_year: int + summary: ThematicRasterSelectionSummary + unsupported_metrics: list[str] + limitation_message: str + generated_at: str diff --git a/backend/app/services/geo_assistant_service.py b/backend/app/services/geo_assistant_service.py index c3fe4668..36a328e7 100644 --- a/backend/app/services/geo_assistant_service.py +++ b/backend/app/services/geo_assistant_service.py @@ -22,8 +22,11 @@ from app.schemas.assistant import ( AssistantTemporalSeries, ) from app.schemas.flood_hazard import FloodHazardSelectionRequest +from app.schemas.thematic_raster import ThematicRasterSelectionRequest from app.services.flood_hazard_acquisition_service import FloodHazardAcquisitionService from app.services.flood_hazard_analysis_service import FloodHazardAnalysisService +from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService +from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService from app.services.vector_feature_service import VectorFeatureService @@ -42,9 +45,17 @@ class GeoAssistantService: ) ESTIMATE_TOPIC_TERMS = { "population": ("bevolk", "inwoner"), + "space_occupation": ("ruimtebeslag",), + "open_space": ("open ruimte",), + "accessibility": ("bereikbaar", "knooppunt"), + "services": ("voorziening",), } ESTIMATE_TOPIC_LABELS = { "population": "bevolkingswaarden", + "space_occupation": "ruimtebeslagoppervlakten", + "open_space": "openruimte-oppervlakten", + "accessibility": "bereikbaarheidsscores", + "services": "voorzieningenscores", } @classmethod @@ -264,6 +275,19 @@ class GeoAssistantService: if dataset.dataset_type == "raster" and dataset.source_name == FloodHazardAcquisitionService.PROVIDER and (area is None or dataset.area_id is None or dataset.area_id == area.id) ] + thematic_candidates = [ + dataset + for dataset in datasets + if dataset.dataset_type == "raster" and dataset.source_name == ThematicRasterAcquisitionService.PROVIDER + and (area is None or dataset.area_id is None or dataset.area_id == area.id) + ] + thematic_by_product: dict[str, Dataset] = {} + for dataset in thematic_candidates: + product_key = str((dataset.source_metadata or {}).get("product_key") or "") + current = thematic_by_product.get(product_key) + if product_key and (current is None or (dataset.imported_at or datetime.min.replace(tzinfo=timezone.utc)) > (current.imported_at or datetime.min.replace(tzinfo=timezone.utc))): + thematic_by_product[product_key] = dataset + thematic_datasets = list(thematic_by_product.values()) warnings: list[str] = [] context_metrics: list[AssistantContextMetric] = [] source_dataset_ids: list[UUID] = [] @@ -323,6 +347,47 @@ class GeoAssistantService: } ) + for dataset in sorted(thematic_datasets, key=lambda item: str((item.source_metadata or {}).get("product_key") or item.name)): + try: + result = ThematicRasterAnalysisService.analyze( + db, + project_id, + dataset.id, + ThematicRasterSelectionRequest(bbox=bbox, area_id=area.id if area is not None else None), + settings=self.settings, + ) + except AppError as exc: + warnings.append(f"{dataset.name}: {exc.message}") + continue + serialized_metrics: list[dict[str, Any]] = [] + for metric in result["summary"]["metrics"]: + item = AssistantContextMetric( + theme=result["theme"], + label=str(metric["metric_label"]), + value=float(metric["metric_value"]), + unit=str(metric["metric_unit"]), + source=ThematicRasterAcquisitionService.ATTRIBUTION, + dataset_id=dataset.id, + observed_at=dataset.observed_at, + is_estimate=bool(metric.get("is_estimate", True)), + ) + context_metrics.append(item) + serialized_metrics.append(item.model_dump(mode="json")) + serialized_metrics[-1]["measurement_quality"] = "resolutiegebonden_bronmeting" + source_dataset_ids.append(dataset.id) + current_context.append( + { + "dataset_name": dataset.name, + "dataset_id": str(dataset.id), + "theme": result["theme"], + "source": ThematicRasterAcquisitionService.ATTRIBUTION, + "observed_at": dataset.observed_at.isoformat() if dataset.observed_at else None, + "metrics": serialized_metrics, + "unsupported_metrics": result["unsupported_metrics"], + "warning": result["limitation_message"], + } + ) + for dataset in sorted( flood_hazard_datasets, key=lambda item: str((item.source_metadata or {}).get("product_key") or item.name), @@ -425,6 +490,7 @@ class GeoAssistantService: "water_volume_available": False, "water_volume_reason": "Geen bathymetrie gekoppeld voor de permanente inhoud van waterlichamen.", "flood_hazard_scenarios_available": bool(flood_hazard_datasets), + "thematic_policy_rasters_available": bool(thematic_datasets), "flood_depth_area_integral_is_concurrent_volume": False, "object_counts_are_supporting_metrics": True, "causal_explanations_available": False, diff --git a/backend/app/services/thematic_raster_acquisition_service.py b/backend/app/services/thematic_raster_acquisition_service.py new file mode 100644 index 00000000..2d204083 --- /dev/null +++ b/backend/app/services/thematic_raster_acquisition_service.py @@ -0,0 +1,597 @@ +from __future__ import annotations + +import hashlib +import json +import math +import time +from dataclasses import dataclass +from datetime import UTC, datetime +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.thematic_raster import ( + ThematicRasterAcquireRequest, + ThematicRasterAcquisitionResult, + ThematicRasterProductRead, +) +from app.services.dataset_service import DatasetService + + +@dataclass(frozen=True) +class ThematicRasterProduct: + key: str + display_name: str + theme: str + metric_kind: str + coverage_id: str + native_resolution_m: float + source_value_unit: str + observation_year: int + source_version: str + catalog_url: str + legend_min_label: str + legend_max_label: str + limitation_message: str + + +class ThematicRasterAcquisitionService: + """Acquire bounded, allowlisted policy rasters from MercatorNet WCS.""" + + PROVIDER = "department_omgeving_thematic_raster" + SOURCE_CRS = "EPSG:31370" + WCS_VERSION = "1.0.0" + NODATA = -9999.0 + WCS_TILE_SIDE_M = 10_000.0 + WCS_REQUEST_INTERVAL_SECONDS = 0.5 + ATTRIBUTION = "Bron: Departement Omgeving, MercatorNet" + LICENSE_NOTE = "Publieke GDI-Vlaanderen bron; bronvermelding en productspecifieke gebruiksvoorwaarden blijven van toepassing." + + @staticmethod + def _products() -> dict[str, ThematicRasterProduct]: + products = ( + ThematicRasterProduct( + key="space_occupation_2025", + display_name="Ruimtebeslag Vlaanderen 2025", + theme="space_occupation", + metric_kind="binary_area", + coverage_id="lu:lu_ruibes_vlaa_2025_v3", + native_resolution_m=10.0, + source_value_unit="class_0_1", + observation_year=2025, + source_version="Toestand 2025 versie 3", + catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/ruimtebeslag-vlaanderen-toestand-2025", + legend_min_label="Geen ruimtebeslag", + legend_max_label="Ruimtebeslag", + limitation_message=( + "Binaire 10 m-kaart volgens de beleidsdefinitie van ruimtebeslag. Celgebaseerde oppervlakte is een " + "resolutiegebonden schatting en is niet gelijk aan uitsluitend bebouwde oppervlakte of verharding." + ), + ), + ThematicRasterProduct( + key="open_space_2022", + display_name="Open ruimte Vlaanderen 2022", + theme="open_space", + metric_kind="binary_area", + coverage_id="lu:lu_openruimte_vlaa_2022_v3", + native_resolution_m=10.0, + source_value_unit="class_0_1", + observation_year=2022, + source_version="Toestand 2022 versie 3", + catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/open-ruimte-vlaanderen-toestand-2022", + legend_min_label="Geen open ruimte", + legend_max_label="Open ruimte", + limitation_message=( + "Binaire 10 m-beleidskaart afgeleid uit landgebruik, ruimtebeslag en kernen. Open ruimte is niet " + "synoniem met natuur, bos, publieke toegankelijkheid of planologische bestemming." + ), + ), + ThematicRasterProduct( + key="population_density_2019", + display_name="Inwonersdichtheid per hectare 2019", + theme="population", + metric_kind="population_density", + coverage_id="ni:ni_inw_ha_vlaa_2019", + native_resolution_m=100.0, + source_value_unit="inhabitants_per_hectare", + observation_year=2019, + source_version="Toestand 2019", + catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/inwonersdichtheid-per-ha-vlaanderen-toestand-2019", + legend_min_label="0 inwoners/ha", + legend_max_label="Hogere dichtheid", + limitation_message=( + "Statistische 1 ha-rasterinschatting voor 2019, gecorrigeerd op statistische-sectorbasis. De som " + "binnen een getekende grens is een rasterraming en geen actuele registertelling." + ), + ), + ThematicRasterProduct( + key="node_value_2022", + display_name="Knooppuntwaarde collectief vervoer 2022", + theme="accessibility", + metric_kind="index_score", + coverage_id="lu:lu_knptw_ha_2022_v3", + native_resolution_m=100.0, + source_value_unit="source_index_score", + observation_year=2022, + source_version="Toestand 2022 versie 3", + catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/knooppuntwaarde-per-ha-toestand-2022", + legend_min_label="Lagere knooppuntwaarde", + legend_max_label="Hogere knooppuntwaarde", + limitation_message=( + "Bronindex per hectare op basis van collectief-vervoerknooppunten en afstandsverval. De score is " + "geen percentage, reistijd, dienstregeling van vandaag of garantie op bereikbaarheid." + ), + ), + ThematicRasterProduct( + key="service_level_2022", + display_name="Totaal voorzieningenniveau 2022", + theme="services", + metric_kind="normalized_score", + coverage_id="lu:lu_totvznv_ha_2022_v3", + native_resolution_m=100.0, + source_value_unit="score_0_1", + observation_year=2022, + source_version="Toestand 2022 versie 3", + catalog_url="https://www.vlaanderen.be/datavindplaats/catalogus/totaal-voorzieningenniveau-toestand-2022", + legend_min_label="Lager voorzieningenniveau", + legend_max_label="Hoger voorzieningenniveau", + limitation_message=( + "Genormaliseerde 0-1 nabijheidsscore voor basis-, regionale en metropolitane voorzieningen in " + "referentiejaar 2022. Dit is geen objecttelling, openingsurencontrole of actuele reistijd." + ), + ), + ) + return {product.key: product for product in products} + + @staticmethod + def list_products() -> list[dict[str, Any]]: + return [ + ThematicRasterProductRead( + key=product.key, + display_name=product.display_name, + theme=product.theme, + metric_kind=product.metric_kind, + coverage_id=product.coverage_id, + native_resolution_m=product.native_resolution_m, + source_crs=ThematicRasterAcquisitionService.SOURCE_CRS, + source_value_unit=product.source_value_unit, + observation_year=product.observation_year, + source_version=product.source_version, + catalog_url=product.catalog_url, + attribution=ThematicRasterAcquisitionService.ATTRIBUTION, + license_note=ThematicRasterAcquisitionService.LICENSE_NOTE, + legend_min_label=product.legend_min_label, + legend_max_label=product.legend_max_label, + limitation_message=product.limitation_message, + ).model_dump() + for product in ThematicRasterAcquisitionService._products().values() + ] + + @staticmethod + def _product(product_key: str) -> ThematicRasterProduct: + product = ThematicRasterAcquisitionService._products().get(product_key.strip().lower()) + if product is None: + raise AppError( + code="THEMATIC_RASTER_PRODUCT_NOT_SUPPORTED", + message="Select a product from the governed Flemish thematic raster registry", + details={"product_key": product_key}, + status_code=422, + ) + return product + + @staticmethod + def _prepared_request(payload: ThematicRasterAcquireRequest, settings: Settings) -> dict[str, Any]: + if not settings.thematic_raster_enabled: + raise AppError(code="THEMATIC_RASTER_NOT_CONFIGURED", message="Official thematic raster acquisition is disabled", status_code=503) + product = ThematicRasterAcquisitionService._product(payload.product_key) + values = (payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y) + if payload.bbox.crs.upper() != "EPSG:4326": + raise AppError(code="INVALID_BBOX_CRS", message="Thematic raster selection requires EPSG:4326", status_code=400) + if not all(math.isfinite(value) for value in values) or values[0] >= values[2] or values[1] >= values[3]: + raise AppError(code="INVALID_BBOX", message="Thematic raster selection must be a finite non-empty rectangle", status_code=400) + + transformer = Transformer.from_crs("EPSG:4326", ThematicRasterAcquisitionService.SOURCE_CRS, always_xy=True) + raw_bounds = transformer.transform_bounds(*values, densify_pts=21) + resolution = product.native_resolution_m + lambert_bounds = ( + math.floor(raw_bounds[0] / resolution) * resolution, + math.floor(raw_bounds[1] / resolution) * resolution, + math.ceil(raw_bounds[2] / resolution) * resolution, + math.ceil(raw_bounds[3] / resolution) * resolution, + ) + width_m = lambert_bounds[2] - lambert_bounds[0] + height_m = lambert_bounds[3] - lambert_bounds[1] + if width_m < settings.thematic_raster_min_side_m or height_m < settings.thematic_raster_min_side_m: + raise AppError( + code="THEMATIC_RASTER_SELECTION_TOO_SMALL", + message=f"Select an area of at least {settings.thematic_raster_min_side_m:g} by {settings.thematic_raster_min_side_m:g} metres", + status_code=422, + ) + if width_m > settings.thematic_raster_max_side_m or height_m > settings.thematic_raster_max_side_m: + raise AppError( + code="THEMATIC_RASTER_SELECTION_TOO_LARGE", + message=f"Select an area no larger than {settings.thematic_raster_max_side_m:g} by {settings.thematic_raster_max_side_m:g} metres", + details={"width_m": width_m, "height_m": height_m}, + status_code=422, + ) + width = max(1, round(width_m / resolution)) + height = max(1, round(height_m / resolution)) + if width * height > settings.thematic_raster_max_pixels: + raise AppError( + code="THEMATIC_RASTER_SELECTION_TOO_LARGE", + message="Thematic raster selection exceeds the configured cell limit", + details={"pixel_count": width * height, "max_pixels": settings.thematic_raster_max_pixels}, + status_code=422, + ) + request_identity = { + "provider": ThematicRasterAcquisitionService.PROVIDER, + "coverage_id": product.coverage_id, + "bbox_epsg4326": [round(float(value), 8) for value in values], + "bbox_epsg31370": [round(float(value), 3) for value in lambert_bounds], + "resolution_m": resolution, + "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() + return { + **request_identity, + "product": product, + "request_hash": request_hash, + "width": width, + "height": height, + } + + @staticmethod + def _wcs_request_url(settings: Settings, product: ThematicRasterProduct, bounds: tuple[float, float, float, float]) -> str: + query = { + "SERVICE": "WCS", + "VERSION": ThematicRasterAcquisitionService.WCS_VERSION, + "REQUEST": "GetCoverage", + "COVERAGE": product.coverage_id, + "CRS": ThematicRasterAcquisitionService.SOURCE_CRS, + "BBOX": ",".join(f"{value:.3f}" for value in bounds), + "RESX": f"{product.native_resolution_m:g}", + "RESY": f"{product.native_resolution_m:g}", + "FORMAT": "image/tiff", + "RESPONSE_CRS": ThematicRasterAcquisitionService.SOURCE_CRS, + } + return f"{settings.thematic_raster_wcs_url}?{urlencode(query)}" + + @staticmethod + def _tile_bounds(prepared: dict[str, Any]) -> list[tuple[float, float, float, float]]: + min_x, min_y, max_x, max_y = prepared["bbox_epsg31370"] + resolution = prepared["product"].native_resolution_m + side = max(resolution, math.floor(ThematicRasterAcquisitionService.WCS_TILE_SIDE_M / resolution) * resolution) + tiles: list[tuple[float, float, float, float]] = [] + y = min_y + while y < max_y: + tile_max_y = min(y + side, max_y) + x = min_x + while x < max_x: + tile_max_x = min(x + side, max_x) + tiles.append((x, y, tile_max_x, tile_max_y)) + x = tile_max_x + y = tile_max_y + return tiles + + @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="THEMATIC_RASTER_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422) + return intersection + + @staticmethod + def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]: + request = Request(request_url, headers={"Accept": "image/tiff,*/*", "User-Agent": "GeoIntel/0.1 bounded-thematic-raster"}) + max_bytes = settings.thematic_raster_max_response_mb * 1024 * 1024 + try: + with (opener or urlopen)(request, timeout=settings.thematic_raster_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="THEMATIC_RASTER_RESPONSE_TOO_LARGE", message="Official raster response exceeds the configured size limit", status_code=502) + content = response.read(max_bytes + 1) + except AppError: + raise + except HTTPError as exc: + preview = exc.read(300).decode("utf-8", errors="replace") + raise AppError( + code="THEMATIC_RASTER_PROVIDER_UNAVAILABLE", + message="The official MercatorNet WCS could not complete the bounded request", + details={"reason": str(exc), "provider_status_code": int(exc.code), "response_preview": preview}, + status_code=502, + ) from exc + except (URLError, TimeoutError, OSError) as exc: + raise AppError( + code="THEMATIC_RASTER_PROVIDER_UNAVAILABLE", + message="The official MercatorNet WCS could not complete the bounded request", + details={"reason": str(exc)}, + status_code=502, + ) from exc + if len(content) > max_bytes: + raise AppError(code="THEMATIC_RASTER_RESPONSE_TOO_LARGE", message="Official raster response exceeds the configured size limit", status_code=502) + if not content.startswith((b"II*\x00", b"MM\x00*")): + preview = content[:300].decode("utf-8", errors="replace") + raise AppError( + code="THEMATIC_RASTER_PROVIDER_INVALID_RESPONSE", + message="The official MercatorNet service did not return a GeoTIFF coverage", + details={"content_type": content_type, "response_preview": preview}, + status_code=502, + ) + return content, content_type + + @staticmethod + def _mosaic(coverages: list[bytes], product: ThematicRasterProduct) -> bytes: + if len(coverages) == 1: + return coverages[0] + try: + import rasterio + from rasterio.io import MemoryFile + from rasterio.merge import merge + except ImportError as exc: + raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio is required to assemble thematic raster tiles", status_code=503) from exc + memories = [MemoryFile(content) for content in coverages] + sources = [] + try: + sources = [memory.open() for memory in memories] + for source in sources: + if source.crs is None or source.crs.to_epsg() != 31370 or source.count != 1: + raise AppError(code="THEMATIC_RASTER_TILE_MISMATCH", message="Thematic raster tiles have incompatible CRS or bands", status_code=502) + if not all(math.isclose(abs(float(value)), product.native_resolution_m, abs_tol=0.05) for value in source.res): + raise AppError(code="THEMATIC_RASTER_TILE_MISMATCH", message="Thematic raster tile resolution differs from the registry", status_code=502) + mosaic, transform = merge(sources, res=(product.native_resolution_m, product.native_resolution_m), nodata=ThematicRasterAcquisitionService.NODATA, dtype="float32") + profile = sources[0].profile.copy() + profile.pop("blockxsize", None) + profile.pop("blockysize", None) + profile.update(driver="GTiff", width=mosaic.shape[2], height=mosaic.shape[1], count=1, dtype="float32", crs=ThematicRasterAcquisitionService.SOURCE_CRS, transform=transform, nodata=ThematicRasterAcquisitionService.NODATA, compress="deflate", predictor=3) + with MemoryFile() as output_memory: + with output_memory.open(**profile) as output: + output.write(mosaic) + return output_memory.read() + except AppError: + raise + except Exception as exc: + raise AppError(code="THEMATIC_RASTER_TILE_MOSAIC_FAILED", message="Thematic raster tiles could not be assembled", details={"reason": str(exc)}, status_code=502) from exc + finally: + for source in sources: + source.close() + for memory in memories: + memory.close() + + @staticmethod + def _fetch_coverage(prepared: dict[str, Any], settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, dict[str, Any]]: + product: ThematicRasterProduct = prepared["product"] + request_urls = [ThematicRasterAcquisitionService._wcs_request_url(settings, product, bounds) for bounds in ThematicRasterAcquisitionService._tile_bounds(prepared)] + coverages: list[bytes] = [] + digest = hashlib.sha256() + content_types: list[str] = [] + for index, request_url in enumerate(request_urls): + if index and opener is None: + time.sleep(ThematicRasterAcquisitionService.WCS_REQUEST_INTERVAL_SECONDS) + content, content_type = ThematicRasterAcquisitionService._fetch(request_url, settings, opener) + digest.update(len(content).to_bytes(8, "big")) + digest.update(content) + coverages.append(content) + content_types.append(content_type) + return ThematicRasterAcquisitionService._mosaic(coverages, product), { + "tile_count": len(request_urls), + "request_urls": request_urls, + "response_content_types": content_types, + "coverage_sha256": digest.hexdigest(), + } + + @staticmethod + def _validate_values(values, product: ThematicRasterProduct) -> None: + import numpy as np + + if values.size == 0: + raise AppError(code="THEMATIC_RASTER_NO_VALID_DATA", message="The official product contains no valid cells in this selection", status_code=422) + minimum = float(values.min()) + maximum = float(values.max()) + if minimum < 0: + raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Official thematic raster contains unexpected negative values", details={"minimum": minimum}, status_code=502) + if product.metric_kind == "binary_area" and not set(np.unique(values).tolist()).issubset({0.0, 1.0}): + raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Binary thematic raster contains classes outside 0 and 1", status_code=502) + if product.metric_kind == "normalized_score" and maximum > 1.0001: + raise AppError(code="THEMATIC_RASTER_INVALID_VALUES", message="Normalized thematic score falls outside the documented 0-1 range", details={"maximum": maximum}, 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 + 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 thematic raster validation", status_code=503) from exc + product: ThematicRasterProduct = prepared["product"] + 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="THEMATIC_RASTER_INVALID_CRS", message="Official thematic raster must use EPSG:31370", status_code=502) + if source.count != 1: + raise AppError(code="THEMATIC_RASTER_INVALID_BANDS", message="Official thematic raster must contain one band", status_code=502) + if not all(math.isclose(abs(float(value)), product.native_resolution_m, abs_tol=0.05) for value in source.res): + raise AppError(code="THEMATIC_RASTER_INVALID_RESOLUTION", message="Official thematic raster resolution differs from the registry", status_code=502) + transformer = Transformer.from_crs("EPSG:4326", ThematicRasterAcquisitionService.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") + raw = np.asarray(band.filled(np.nan), dtype="float32") + invalid = np.ma.getmaskarray(band) | ~np.isfinite(raw) + if source.nodata is not None: + invalid |= np.isclose(raw, float(source.nodata)) + normalized = np.ma.array(raw, mask=invalid) + values = normalized.compressed().astype("float64") + ThematicRasterAcquisitionService._validate_values(values, product) + profile = source.profile.copy() + profile.pop("blockxsize", None) + profile.pop("blockysize", None) + profile.update(driver="GTiff", width=normalized.shape[1], height=normalized.shape[0], count=1, dtype="float32", crs=ThematicRasterAcquisitionService.SOURCE_CRS, transform=transform, nodata=ThematicRasterAcquisitionService.NODATA, compress="deflate", predictor=3) + with MemoryFile() as output_memory: + with output_memory.open(**profile) as output: + output.write(normalized.filled(ThematicRasterAcquisitionService.NODATA), 1) + normalized_content = output_memory.read() + return normalized_content, { + "width": int(normalized.shape[1]), + "height": int(normalized.shape[0]), + "valid_pixel_count": int(values.size), + "nodata_value": ThematicRasterAcquisitionService.NODATA, + "resolution_m": product.native_resolution_m, + "minimum_value": float(values.min()), + "maximum_value": float(values.max()), + "p02_value": float(np.percentile(values, 2)), + "p98_value": float(np.percentile(values, 98)), + } + except AppError: + raise + except Exception as exc: + raise AppError(code="THEMATIC_RASTER_INVALID", message="The official thematic raster could not be validated", details={"reason": str(exc)}, status_code=502) from exc + + @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 == ThematicRasterAcquisitionService.PROVIDER, Dataset.status == "ready") + .order_by(Dataset.imported_at.desc()) + .first() + ) + return candidate if candidate and candidate.storage_path and Path(candidate.storage_path).is_file() else None + + @staticmethod + def acquire(db, project_id: UUID, payload: ThematicRasterAcquireRequest, *, settings: Settings | None = None, opener: Callable[..., Any] | None = None) -> dict[str, Any]: + resolved_settings = settings or get_settings() + prepared = ThematicRasterAcquisitionService._prepared_request(payload, resolved_settings) + product: ThematicRasterProduct = prepared["product"] + scope_geometry = ThematicRasterAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"]) + filename = f"thematic_{product.key}_{prepared['request_hash'][:12]}.tif" + if not payload.force_refresh: + cached = ThematicRasterAcquisitionService._cached_dataset(db, project_id, filename) + if cached is not None: + metadata = cached.source_metadata or {} + raster_metadata = cached.metadata_json or {} + return ThematicRasterAcquisitionResult( + output_dataset_id=cached.id, + reused=True, + provider=ThematicRasterAcquisitionService.PROVIDER, + product_key=product.key, + display_name=product.display_name, + theme=product.theme, + metric_kind=product.metric_kind, + coverage_id=product.coverage_id, + resolution_m=product.native_resolution_m, + width=int(raster_metadata.get("width", prepared["width"])), + height=int(raster_metadata.get("height", prepared["height"])), + valid_pixel_count=int(metadata.get("valid_pixel_count", 0)), + bbox_epsg4326=prepared["bbox_epsg4326"], + bbox_epsg31370=prepared["bbox_epsg31370"], + observation_year=product.observation_year, + source_value_unit=product.source_value_unit, + attribution=ThematicRasterAcquisitionService.ATTRIBUTION, + limitation_message=product.limitation_message, + ).model_dump(mode="json") + + coverage, transfer = ThematicRasterAcquisitionService._fetch_coverage(prepared, resolved_settings, opener) + normalized, validation = ThematicRasterAcquisitionService._normalize_raster(coverage, scope_geometry, prepared) + acquired_at = datetime.now(UTC) + observed_at = datetime(product.observation_year, 12, 31, 23, 59, 59, tzinfo=UTC) + scope_key = str(payload.area_id) if payload.area_id else prepared["request_hash"][:24] + dataset = DatasetService.import_raster_bytes( + db, + project_id=project_id, + area_id=payload.area_id, + filename=filename, + content=normalized, + source=f"Departement Omgeving MercatorNet WCS {product.coverage_id}", + source_name=ThematicRasterAcquisitionService.PROVIDER, + temporal_series_key=f"department-omgeving:thematic-raster:{product.key}:{scope_key}", + observed_at=observed_at, + valid_from=datetime(product.observation_year, 1, 1, tzinfo=UTC), + valid_to=observed_at, + temporal_granularity="year", + source_version=product.source_version, + source_metadata={ + "provider": ThematicRasterAcquisitionService.PROVIDER, + "service": "WCS", + "service_version": ThematicRasterAcquisitionService.WCS_VERSION, + "product_key": product.key, + "product_display_name": product.display_name, + "theme": product.theme, + "metric_kind": product.metric_kind, + "coverage_id": product.coverage_id, + "native_resolution_m": product.native_resolution_m, + "analysis_resolution_m": product.native_resolution_m, + "source_crs": ThematicRasterAcquisitionService.SOURCE_CRS, + "source_value_unit": product.source_value_unit, + "observation_year": product.observation_year, + "observation_date_precision": "year", + "valid_pixel_count": validation["valid_pixel_count"], + "minimum_value": validation["minimum_value"], + "maximum_value": validation["maximum_value"], + "render_min_value": validation["p02_value"], + "render_max_value": validation["p98_value"], + "bbox_epsg4326": prepared["bbox_epsg4326"], + "bbox_epsg31370": prepared["bbox_epsg31370"], + "catalog_url": product.catalog_url, + "attribution": ThematicRasterAcquisitionService.ATTRIBUTION, + "license_note": ThematicRasterAcquisitionService.LICENSE_NOTE, + "legend_min_label": product.legend_min_label, + "legend_max_label": product.legend_max_label, + "coverage_scope": "municipality" if payload.area_id else "bounded_selection", + }, + provenance_metadata={ + "acquisition": "explicit_bounded_tiled_wcs_coverage", + "acquired_at": acquired_at.isoformat(), + "request_hash": prepared["request_hash"], + "tile_count": transfer["tile_count"], + "tile_request_urls": transfer["request_urls"], + "response_content_types": transfer["response_content_types"], + "coverage_sha256": transfer["coverage_sha256"], + "normalized_sha256": hashlib.sha256(normalized).hexdigest(), + "bbox_epsg4326": prepared["bbox_epsg4326"], + "bbox_epsg31370": prepared["bbox_epsg31370"], + "clipped_to_area_id": str(payload.area_id) if payload.area_id else None, + "validation": validation, + "limitation_message": product.limitation_message, + }, + ) + return ThematicRasterAcquisitionResult( + output_dataset_id=dataset.id, + reused=False, + provider=ThematicRasterAcquisitionService.PROVIDER, + product_key=product.key, + display_name=product.display_name, + theme=product.theme, + metric_kind=product.metric_kind, + coverage_id=product.coverage_id, + resolution_m=product.native_resolution_m, + width=validation["width"], + height=validation["height"], + valid_pixel_count=validation["valid_pixel_count"], + bbox_epsg4326=prepared["bbox_epsg4326"], + bbox_epsg31370=prepared["bbox_epsg31370"], + observation_year=product.observation_year, + source_value_unit=product.source_value_unit, + attribution=ThematicRasterAcquisitionService.ATTRIBUTION, + limitation_message=product.limitation_message, + ).model_dump(mode="json") diff --git a/backend/app/services/thematic_raster_analysis_service.py b/backend/app/services/thematic_raster_analysis_service.py new file mode 100644 index 00000000..473d4a58 --- /dev/null +++ b/backend/app/services/thematic_raster_analysis_service.py @@ -0,0 +1,260 @@ +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.thematic_raster import ( + ThematicRasterMetric, + ThematicRasterSelectionRequest, + ThematicRasterSelectionResponse, + ThematicRasterSelectionSummary, +) +from app.services.thematic_raster_acquisition_service import ( + ThematicRasterAcquisitionService, + ThematicRasterProduct, +) + + +class ThematicRasterAnalysisService: + @staticmethod + def _load_dataset(db, project_id: UUID, dataset_id: UUID) -> Dataset: + dataset = db.get(Dataset, dataset_id) + if not dataset or dataset.project_id != project_id: + raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) + if dataset.dataset_type != "raster" or dataset.source_name != ThematicRasterAcquisitionService.PROVIDER: + raise AppError( + code="INVALID_THEMATIC_RASTER_DATASET", + message="Thematic analysis requires a governed Departement Omgeving raster dataset", + status_code=400, + ) + if dataset.status != "ready" or not dataset.storage_path or not Path(dataset.storage_path).is_file(): + raise AppError(code="DATASET_FILE_MISSING", message="Persisted thematic raster file is unavailable", status_code=404) + return dataset + + @staticmethod + def _product(dataset: Dataset) -> ThematicRasterProduct: + source_metadata = dataset.source_metadata or {} + product = ThematicRasterAcquisitionService._products().get(str(source_metadata.get("product_key") or "")) + if product is None or source_metadata.get("coverage_id") != product.coverage_id: + raise AppError(code="INVALID_THEMATIC_RASTER_METADATA", message="Thematic raster provenance is incomplete", status_code=409) + return product + + @staticmethod + def _selection_geometry(db, project_id: UUID, payload: ThematicRasterSelectionRequest): + selection = box(payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y) + if payload.area_id is None: + return selection + area = db.get(Area, payload.area_id) + if not area: + raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404) + if area.project_id != project_id: + raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400) + selection = selection.intersection(to_shape(area.geometry)) + if selection.is_empty or selection.area <= 0: + raise AppError(code="THEMATIC_RASTER_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422) + return selection + + @staticmethod + def _unsupported_metrics(product: ThematicRasterProduct) -> list[str]: + if product.metric_kind == "binary_area": + return ["object_count", "parcel_area", "current_land_use"] + if product.metric_kind == "population_density": + return ["current_population", "household_count", "address_level_population"] + if product.metric_kind == "index_score": + return ["travel_time_minutes", "current_timetable", "stop_count"] + return ["facility_count", "opening_hours", "current_service_availability"] + + @staticmethod + def analyze( + db, + project_id: UUID, + dataset_id: UUID, + payload: ThematicRasterSelectionRequest, + *, + settings: Settings | None = None, + ) -> dict: + resolved_settings = settings or get_settings() + dataset = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id) + product = ThematicRasterAnalysisService._product(dataset) + selection_4326 = ThematicRasterAnalysisService._selection_geometry(db, project_id, payload) + try: + import numpy as np + import rasterio + from rasterio.features import geometry_mask + from rasterio.mask import mask + except ImportError as exc: + raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for thematic raster analysis", status_code=503) from exc + + try: + with rasterio.open(dataset.storage_path) as source: + if source.crs is None or source.crs.to_epsg() != 31370: + raise AppError(code="INVALID_DATASET_CRS", message="Thematic raster CRS must be EPSG:31370", status_code=409) + transformer = Transformer.from_crs("EPSG:4326", source.crs, always_xy=True) + selection_metric = shapely_transform(transformer.transform, selection_4326) + analysis_geometry = selection_metric.intersection(box(*source.bounds)) + if analysis_geometry.is_empty or analysis_geometry.area <= 0: + raise AppError(code="THEMATIC_RASTER_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted thematic raster", status_code=422) + min_x, min_y, max_x, max_y = analysis_geometry.bounds + expected_cells = math.ceil((max_x - min_x) / abs(source.res[0])) * math.ceil((max_y - min_y) / abs(source.res[1])) + if expected_cells > resolved_settings.thematic_raster_max_pixels: + raise AppError( + code="THEMATIC_RASTER_SELECTION_TOO_LARGE", + message="Thematic raster analysis exceeds the configured cell limit", + details={"pixel_count": expected_cells, "max_pixels": resolved_settings.thematic_raster_max_pixels}, + status_code=422, + ) + clipped, clipped_transform = mask(source, [mapping(analysis_geometry)], crop=True, filled=False, indexes=[1]) + band = np.ma.asarray(clipped[0], dtype="float64") + raw = band.filled(np.nan) + selected = geometry_mask([mapping(analysis_geometry)], out_shape=band.shape, transform=clipped_transform, invert=True) + valid = selected & ~np.ma.getmaskarray(band) & np.isfinite(raw) + if source.nodata is not None: + valid &= ~np.isclose(raw, float(source.nodata)) + values = raw[valid] + ThematicRasterAcquisitionService._validate_values(values, product) + selected_cell_count = int(selected.sum()) + valid_cell_count = int(values.size) + resolution_x = abs(float(source.res[0])) + resolution_y = abs(float(source.res[1])) + cell_area_m2 = resolution_x * resolution_y + except AppError: + raise + except Exception as exc: + raise AppError( + code="THEMATIC_RASTER_ANALYSIS_FAILED", + message="The persisted thematic raster could not be analysed", + details={"reason": str(exc)}, + status_code=500, + ) from exc + + def metric(key: str, label: str, value: float, unit: str, method: str, *, estimate: bool = True) -> ThematicRasterMetric: + return ThematicRasterMetric( + metric_key=key, + metric_label=label, + metric_value=round(float(value), 4), + metric_unit=unit, + aggregation_method=method, + is_estimate=estimate, + ) + + if product.metric_kind == "binary_area": + positive_count = int(np.count_nonzero(values >= 0.5)) + positive_area_ha = positive_count * cell_area_m2 / 10_000.0 + positive_share = positive_count / max(1, valid_cell_count) * 100.0 + label = "Ruimtebeslag" if product.theme == "space_occupation" else "Open ruimte" + metrics = [ + metric(f"{product.theme}_area_ha", f"{label} in selectie", positive_area_ha, "ha", "positive_source_cells_times_cell_area"), + metric(f"{product.theme}_share_pct", f"Aandeel {label.lower()}", positive_share, "%", "positive_source_cells_divided_by_valid_selected_cells"), + metric("valid_raster_area_ha", "Rasteroppervlakte met bronwaarde", valid_cell_count * cell_area_m2 / 10_000.0, "ha", "valid_selected_cells_times_cell_area"), + ] + elif product.metric_kind == "population_density": + estimated_population = float(values.sum() * (cell_area_m2 / 10_000.0)) + metrics = [ + metric("estimated_inhabitants", "Geraamd aantal inwoners (2019)", estimated_population, "inwoners", "sum_density_times_selected_cell_area_hectares"), + metric("population_density_mean_per_ha", "Gemiddelde inwonersdichtheid", values.mean(), "inwoners/ha", "mean_valid_one_hectare_source_cells"), + metric("population_density_p90_per_ha", "90e percentiel inwonersdichtheid", np.percentile(values, 90), "inwoners/ha", "percentile_90_valid_source_cells"), + ] + else: + unit = "score" if product.metric_kind == "index_score" else "score (0-1)" + label = "Knooppuntwaarde" if product.metric_kind == "index_score" else "Voorzieningenniveau" + metrics = [ + metric(f"{product.theme}_mean", f"Gemiddelde {label.lower()}", values.mean(), unit, "mean_valid_source_cells"), + metric(f"{product.theme}_p10", f"10e percentiel {label.lower()}", np.percentile(values, 10), unit, "percentile_10_valid_source_cells"), + metric(f"{product.theme}_median", f"Mediaan {label.lower()}", np.percentile(values, 50), unit, "median_valid_source_cells"), + metric(f"{product.theme}_p90", f"90e percentiel {label.lower()}", np.percentile(values, 90), unit, "percentile_90_valid_source_cells"), + ] + + primary = metrics[0] + response = ThematicRasterSelectionResponse( + dataset_id=dataset.id, + product_key=product.key, + theme=product.theme, + metric_kind=product.metric_kind, + selection_bbox=payload.bbox, + selection_area_id=payload.area_id, + selected_cell_count=selected_cell_count, + valid_cell_count=valid_cell_count, + coverage_ratio=round(valid_cell_count / max(1, selected_cell_count), 6), + resolution_m=round(max(resolution_x, resolution_y), 4), + observation_year=product.observation_year, + summary=ThematicRasterSelectionSummary( + metric_label=primary.metric_label, + metric_value=primary.metric_value, + metric_unit=primary.metric_unit, + aggregation_method=primary.aggregation_method, + primary_metric_key=primary.metric_key, + metrics=metrics, + ), + unsupported_metrics=ThematicRasterAnalysisService._unsupported_metrics(product), + limitation_message=product.limitation_message, + generated_at=datetime.now(UTC).isoformat(), + ) + return response.model_dump(mode="json") + + @staticmethod + def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800) -> bytes: + dataset = ThematicRasterAnalysisService._load_dataset(db, project_id, dataset_id) + product = ThematicRasterAnalysisService._product(dataset) + try: + import numpy as np + import rasterio + from PIL import Image + from rasterio.enums import Resampling + except ImportError as exc: + raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio, numpy and Pillow are required for thematic raster rendering", status_code=503) from exc + palettes = { + "space_occupation": np.asarray([[251, 231, 211], [190, 62, 51]], dtype="float64"), + "open_space": np.asarray([[221, 238, 219], [38, 122, 70]], dtype="float64"), + "population": np.asarray([[238, 231, 246], [103, 58, 151]], dtype="float64"), + "accessibility": np.asarray([[233, 241, 244], [15, 118, 110]], dtype="float64"), + "services": np.asarray([[255, 244, 191], [182, 109, 22]], dtype="float64"), + } + try: + with rasterio.open(dataset.storage_path) as source: + scale = min(1.0, max_dimension / max(source.width, source.height)) + width = max(1, round(source.width * scale)) + height = max(1, round(source.height * scale)) + resampling = Resampling.nearest if product.metric_kind == "binary_area" else Resampling.bilinear + data = source.read(1, out_shape=(height, width), masked=True, resampling=resampling) + values = np.asarray(data.filled(np.nan), dtype="float64") + valid = np.isfinite(values) & ~np.ma.getmaskarray(data) + if source.nodata is not None: + valid &= ~np.isclose(values, float(source.nodata)) + if product.metric_kind == "binary_area": + valid &= values >= 0.5 + normalized = np.where(valid, 1.0, 0.0) + else: + source_metadata = dataset.source_metadata or {} + lower = float(source_metadata.get("render_min_value", np.nanpercentile(values[valid], 2) if valid.any() else 0.0)) + upper = float(source_metadata.get("render_max_value", np.nanpercentile(values[valid], 98) if valid.any() else 1.0)) + if upper <= lower: + upper = lower + 1.0 + normalized = np.clip((values - lower) / (upper - lower), 0.0, 1.0) + colors = palettes[product.theme] + rgba = np.zeros((height, width, 4), dtype="uint8") + for channel in range(3): + rgba[:, :, channel] = (colors[0, channel] + normalized * (colors[1, channel] - colors[0, channel])).astype("uint8") + rgba[:, :, 3] = np.where(valid, 205, 0).astype("uint8") + output = io.BytesIO() + Image.fromarray(rgba).save(output, format="PNG", optimize=True) + return output.getvalue() + except AppError: + raise + except Exception as exc: + raise AppError( + code="THEMATIC_RASTER_PREVIEW_FAILED", + message="The persisted thematic raster could not be rendered", + details={"reason": str(exc)}, + status_code=500, + ) from exc diff --git a/backend/app/services/vector_feature_service.py b/backend/app/services/vector_feature_service.py index 65297dd3..3a119ddb 100644 --- a/backend/app/services/vector_feature_service.py +++ b/backend/app/services/vector_feature_service.py @@ -29,6 +29,7 @@ FULL_AREA_CLIPPED_OPERATOR_TOOLS = { "provision_regional_bwk_natura2000.py", "provision_agricultural_parcel_history.py", "provision_buildings_addresses_register.py", + "provision_mol_soil_map.py", } SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS = { @@ -101,6 +102,16 @@ SEMANTIC_SELECTION_METRICS: dict[str, tuple[dict[str, Any], ...]] = { ), "nature_value": (), "agriculture": (), + "soil": ( + { + "metric_key": "soil_mapped_area", + "method": "intersection_area", + "label": "Bodemkaartoppervlakte", + "unit": "ha", + "geometry_dimension": 2, + "warning": "Historische bodemkartering op schaal 1:20.000; actuele lokale bodem- en drainagetoestand kan afwijken.", + }, + ), } SEMANTIC_COUNT_LABELS = { @@ -112,6 +123,7 @@ SEMANTIC_COUNT_LABELS = { "parcels": "Percelen", "nature_value": "BWK-kaartvlakken", "agriculture": "Landbouwgebruikspercelen", + "soil": "Bodemkaartvlakken", } # Sprint 205 initially normalized two official comma-separated ALZ group labels @@ -158,6 +170,8 @@ class VectorFeatureService: "landbouw": "agriculture", "landbouwgebruik": "agriculture", "building_registry": "buildings", + "soil_map": "soil", + "bodem": "soil", } for candidate in candidates: if not isinstance(candidate, str) or not candidate.strip(): diff --git a/backend/tests/test_sprint213_thematic_rasters.py b/backend/tests/test_sprint213_thematic_rasters.py new file mode 100644 index 00000000..019928b3 --- /dev/null +++ b/backend/tests/test_sprint213_thematic_rasters.py @@ -0,0 +1,339 @@ +from __future__ import annotations + +from pathlib import Path +from types import SimpleNamespace +from uuid import uuid4 + +import numpy as np +import pytest +import rasterio +from fastapi.testclient import TestClient +from pyproj import Transformer +from rasterio.io import MemoryFile +from rasterio.transform import from_origin + +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 Dataset, Job, Project +from app.schemas.thematic_raster import ThematicRasterAcquireRequest, ThematicRasterSelectionRequest +from app.schemas.assistant import AssistantQueryRequest +from app.services.geo_assistant_service import GeoAssistantService +from app.services.thematic_raster_acquisition_service import ThematicRasterAcquisitionService +from app.services.thematic_raster_analysis_service import ThematicRasterAnalysisService +from app.services.dataset_service import DatasetService + + +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 + + def all(self): + return self.result if isinstance(self.result, list) else [] + + +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): + self.content = content + self.headers = {"Content-Type": "image/tiff", "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 payload(product_key: str = "space_occupation_2025", *, side_m: float = 1000.0) -> ThematicRasterAcquireRequest: + transformer = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) + min_x, min_y = transformer.transform(200_000, 210_000) + max_x, max_y = transformer.transform(200_000 + side_m, 210_000 + side_m) + return ThematicRasterAcquireRequest( + bbox={"min_x": min_x, "min_y": min_y, "max_x": max_x, "max_y": max_y, "crs": "EPSG:4326"}, + product_key=product_key, + force_refresh=True, + ) + + +def raster_bytes(values: np.ndarray, resolution: float, *, nodata: float = -9999.0) -> bytes: + with MemoryFile() as memory: + with memory.open( + driver="GTiff", + width=values.shape[1], + height=values.shape[0], + count=1, + dtype=str(values.dtype), + crs="EPSG:31370", + transform=from_origin(200_000, 210_000 + values.shape[0] * resolution, resolution, resolution), + nodata=nodata, + ) as output: + output.write(values, 1) + return memory.read() + + +def test_registry_contains_five_governed_non_water_policy_products() -> None: + products = ThematicRasterAcquisitionService.list_products() + + assert [item["key"] for item in products] == [ + "space_occupation_2025", + "open_space_2022", + "population_density_2019", + "node_value_2022", + "service_level_2022", + ] + assert {item["theme"] for item in products} == {"space_occupation", "open_space", "population", "accessibility", "services"} + assert {item["native_resolution_m"] for item in products} == {10.0, 100.0} + assert all(item["coverage_id"].startswith(("lu:", "ni:")) for item in products) + assert all(item["source_crs"] == "EPSG:31370" for item in products) + assert all(item["attribution"] and item["license_note"] and item["limitation_message"] for item in products) + + +def test_request_is_bounded_allowlisted_and_uses_native_wcs_resolution() -> None: + settings = Settings(_env_file=None) + prepared = ThematicRasterAcquisitionService._prepared_request(payload("population_density_2019"), settings) + url = ThematicRasterAcquisitionService._wcs_request_url(settings, prepared["product"], tuple(prepared["bbox_epsg31370"])) + + assert "VERSION=1.0.0" in url + assert "COVERAGE=ni%3Ani_inw_ha_vlaa_2019" in url + assert "RESX=100" in url and "RESY=100" in url + assert prepared["width"] * prepared["height"] <= settings.thematic_raster_max_pixels + + with pytest.raises(AppError) as exc_info: + ThematicRasterAcquisitionService._prepared_request(payload("arbitrary_remote_layer"), settings) + assert exc_info.value.code == "THEMATIC_RASTER_PRODUCT_NOT_SUPPORTED" + + +def test_binary_and_normalized_products_fail_closed_on_invalid_values() -> None: + binary = ThematicRasterAcquisitionService._product("space_occupation_2025") + score = ThematicRasterAcquisitionService._product("service_level_2022") + + with pytest.raises(AppError, match="Binary"): + ThematicRasterAcquisitionService._validate_values(np.asarray([0.0, 2.0]), binary) + with pytest.raises(AppError, match="0-1"): + ThematicRasterAcquisitionService._validate_values(np.asarray([0.2, 1.2]), score) + + +def test_acquisition_clips_validates_and_delegates_persistence(monkeypatch) -> None: + project_id, output_dataset_id = uuid4(), uuid4() + db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")}) + content = raster_bytes(np.ones((100, 100), dtype="float32"), 10.0) + captured: dict = {} + + def fake_import(_db, **kwargs): + captured.update(kwargs) + return SimpleNamespace(id=output_dataset_id) + + monkeypatch.setattr(DatasetService, "import_raster_bytes", fake_import) + result = ThematicRasterAcquisitionService.acquire( + db, + project_id, + payload(side_m=1000.0), + settings=Settings(_env_file=None), + opener=lambda *_args, **_kwargs: FakeResponse(content), + ) + + assert result["output_dataset_id"] == str(output_dataset_id) + assert captured["source_name"] == ThematicRasterAcquisitionService.PROVIDER + assert captured["source_metadata"]["product_key"] == "space_occupation_2025" + assert captured["source_metadata"]["metric_kind"] == "binary_area" + assert captured["source_metadata"]["valid_pixel_count"] > 9_800 + assert captured["provenance_metadata"]["acquisition"] == "explicit_bounded_tiled_wcs_coverage" + assert len(captured["provenance_metadata"]["normalized_sha256"]) == 64 + + +def test_binary_area_analysis_returns_hectares_and_share(tmp_path) -> None: + project_id, dataset_id = uuid4(), uuid4() + values = np.zeros((10, 10), dtype="float32") + values[:, :5] = 1.0 + path = tmp_path / "space.tif" + path.write_bytes(raster_bytes(values, 10.0)) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="space.tif", + dataset_type="raster", + source="official", + source_name=ThematicRasterAcquisitionService.PROVIDER, + source_metadata={"product_key": "space_occupation_2025", "coverage_id": "lu:lu_ruibes_vlaa_2025_v3"}, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + result = ThematicRasterAnalysisService.analyze( + db, + project_id, + dataset_id, + ThematicRasterSelectionRequest(bbox=payload(side_m=100.0).bbox), + ) + metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]} + + assert result["valid_cell_count"] == 100 + assert metrics["space_occupation_area_ha"]["metric_value"] == pytest.approx(0.5) + assert metrics["space_occupation_share_pct"]["metric_value"] == pytest.approx(50.0) + assert "object_count" in result["unsupported_metrics"] + + +def test_population_analysis_sums_one_hectare_density_cells_without_claiming_current_counts(tmp_path) -> None: + project_id, dataset_id = uuid4(), uuid4() + values = np.asarray([[10.0, 20.0], [30.0, 40.0]], dtype="float32") + path = tmp_path / "population.tif" + path.write_bytes(raster_bytes(values, 100.0)) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="population.tif", + dataset_type="raster", + source="official", + source_name=ThematicRasterAcquisitionService.PROVIDER, + source_metadata={"product_key": "population_density_2019", "coverage_id": "ni:ni_inw_ha_vlaa_2019"}, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + result = ThematicRasterAnalysisService.analyze( + db, + project_id, + dataset_id, + ThematicRasterSelectionRequest(bbox=payload("population_density_2019", side_m=200.0).bbox), + ) + metrics = {item["metric_key"]: item for item in result["summary"]["metrics"]} + + assert metrics["estimated_inhabitants"]["metric_value"] == pytest.approx(100.0) + assert metrics["population_density_mean_per_ha"]["metric_value"] == pytest.approx(25.0) + assert metrics["estimated_inhabitants"]["is_estimate"] is True + assert "current_population" in result["unsupported_metrics"] + + +def test_assistant_context_receives_persisted_thematic_metrics(tmp_path) -> None: + project_id, dataset_id = uuid4(), uuid4() + values = np.asarray([[10.0, 20.0], [30.0, 40.0]], dtype="float32") + path = tmp_path / "assistant-population.tif" + path.write_bytes(raster_bytes(values, 100.0)) + project = Project(id=project_id, name="Mol", region="Mol") + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="population.tif", + dataset_type="raster", + source="official", + source_name=ThematicRasterAcquisitionService.PROVIDER, + source_metadata={"product_key": "population_density_2019", "coverage_id": "ni:ni_inw_ha_vlaa_2019"}, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Project, project_id): project, (Dataset, dataset_id): dataset}, query_result=[dataset]) + context, metrics, _series, dataset_ids, _warnings, _scope = GeoAssistantService(Settings(_env_file=None))._build_context( + db, + project_id=project_id, + payload=AssistantQueryRequest(question="Hoeveel inwoners?", bbox=payload("population_density_2019", side_m=200.0).bbox), + ) + + assert any(metric.theme == "population" and metric.label.startswith("Geraamd aantal") for metric in metrics) + assert dataset_id in dataset_ids + assert context["rules"]["thematic_policy_rasters_available"] is True + + +def test_index_renderer_returns_browser_png(tmp_path) -> None: + project_id, dataset_id = uuid4(), uuid4() + values = np.linspace(0.1, 4.0, 100, dtype="float32").reshape((10, 10)) + path = tmp_path / "node.tif" + path.write_bytes(raster_bytes(values, 100.0)) + dataset = Dataset( + id=dataset_id, + project_id=project_id, + name="node.tif", + dataset_type="raster", + source="official", + source_name=ThematicRasterAcquisitionService.PROVIDER, + source_metadata={ + "product_key": "node_value_2022", + "coverage_id": "lu:lu_knptw_ha_2022_v3", + "render_min_value": 0.1, + "render_max_value": 4.0, + }, + status="ready", + storage_path=str(path), + ) + db = FakeSession({(Dataset, dataset_id): dataset}) + + assert ThematicRasterAnalysisService.render_png(db, project_id, dataset_id).startswith(b"\x89PNG\r\n\x1a\n") + + +def test_api_uses_canonical_envelopes(monkeypatch) -> None: + project_id, dataset_id = uuid4(), uuid4() + db = FakeSession({(Project, project_id): Project(id=project_id, name="Mol")}) + monkeypatch.setattr( + ThematicRasterAcquisitionService, + "acquire", + lambda *_args, **_kwargs: {"output_dataset_id": str(dataset_id), "provider": ThematicRasterAcquisitionService.PROVIDER}, + ) + monkeypatch.setattr( + ThematicRasterAnalysisService, + "analyze", + lambda *_args, **_kwargs: {"dataset_id": str(dataset_id), "theme": "population", "summary": {"metric_value": 10.0}}, + ) + app.dependency_overrides[get_db] = lambda: db + try: + client = TestClient(app) + products = client.get(f"/api/v1/projects/{project_id}/datasets/thematic-raster/products") + acquisition = client.post( + f"/api/v1/projects/{project_id}/datasets/thematic-raster/acquire", + json=payload().model_dump(mode="json"), + ) + selection = client.post( + f"/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select", + json={"bbox": payload().bbox.model_dump()}, + ) + finally: + app.dependency_overrides.clear() + + assert products.status_code == 200 and set(products.json()) == {"data"} + assert products.json()["data"]["total"] == 5 + assert acquisition.status_code == 200 and set(acquisition.json()) == {"data"} + assert acquisition.json()["data"]["job_type"] == "raster.thematic.acquire" + assert selection.status_code == 200 and selection.json()["data"]["theme"] == "population" + assert any(isinstance(item, Job) for item in db.added) diff --git a/backend/tests/test_sprint214_dov_soil_map.py b/backend/tests/test_sprint214_dov_soil_map.py new file mode 100644 index 00000000..ca52ebc3 --- /dev/null +++ b/backend/tests/test_sprint214_dov_soil_map.py @@ -0,0 +1,183 @@ +from __future__ import annotations + +import importlib.util +from pathlib import Path +from uuid import uuid4 + +import pytest +from shapely.geometry import box, mapping, shape +from shapely.ops import transform as transform_geometry + +from app.models import Dataset +from app.services.vector_feature_service import VectorFeatureService + + +ROOT = Path(__file__).resolve().parents[2] + + +def load_operator(): + path = ROOT / "scripts" / "provision_mol_soil_map.py" + spec = importlib.util.spec_from_file_location("dov_soil_map_operator", path) + assert spec and spec.loader + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module + + +class FakeResponse: + def __init__(self, payload: dict, url: str): + self._payload = payload + self.url = url + self.content = b'{"type":"FeatureCollection"}' + + def raise_for_status(self) -> None: + return None + + def json(self) -> dict: + return self._payload + + +class FakeSession: + def __init__(self, pages: list[dict]): + self.pages = pages + self.calls: list[dict] = [] + + def get(self, _url: str, *, params: dict, timeout: int): + self.calls.append({"params": dict(params), "timeout": timeout}) + return FakeResponse(self.pages[len(self.calls) - 1], f"https://example.test/page/{len(self.calls)}") + + +def soil_feature(module, feature_id: str = "bodemtypes.1") -> dict: + return { + "type": "Feature", + "id": feature_id, + "geometry": mapping(box(5.0, 51.0, 5.02, 51.02)), + "properties": { + "gid": 1, + "id_kaartvlak": 10, + "Bodemtype": "Zeg", + "Unibodemtype": "Zeg", + "Bodemserie": "Zeg", + "Beknopte_omschrijving_bodemserie": "Natte zandbodem", + "Gegeneraliseerde_legende": "Nat zand", + "Textuurklasse_code": "Z", + "Textuurklasse": "zand", + "Drainageklasse_code": "e", + "Drainageklasse": "nat", + "Profielontwikkelingsgroep_code": "g", + "Profielontwikkelingsgroep": "humus B horizont", + "Eenduidige_legende_titel": "bodemserie Zeg", + }, + } + + +def test_wfs_pagination_is_bounded_complete_and_deterministic() -> None: + module = load_operator() + feature = soil_feature(module) + pages = [ + { + "type": "FeatureCollection", + "numberMatched": 3, + "numberReturned": 2, + "features": [feature, {**feature, "id": "bodemtypes.2"}], + }, + { + "type": "FeatureCollection", + "numberMatched": 3, + "numberReturned": 1, + "features": [{**feature, "id": "bodemtypes.3"}], + }, + ] + session = FakeSession(pages) + + result = list( + module.iter_wfs_pages( + session, + (196000.0, 205000.0, 211000.0, 224000.0), + page_limit=2, + timeout=30, + ) + ) + + assert len(result) == 2 + assert [call["params"]["startIndex"] for call in session.calls] == ["0", "2"] + assert all(call["params"]["typeNames"] == "bodemkaart:bodemtypes" for call in session.calls) + assert all(call["params"]["bbox"].endswith("EPSG:31370") for call in session.calls) + assert all(call["params"]["sortBy"] == "gid" for call in session.calls) + + +def test_soil_feature_is_exactly_clipped_and_keeps_governed_properties() -> None: + module = load_operator() + boundary_wgs84 = box(5.005, 51.005, 5.015, 51.015) + boundary_lambert72 = transform_geometry(module.TO_LAMBERT72.transform, boundary_wgs84) + + normalized, was_clipped = module.normalize_feature(soil_feature(module), boundary_lambert72) + + assert normalized is not None and was_clipped is True + persisted_geometry = shape(normalized["geometry"]) + assert persisted_geometry.within(boundary_wgs84.buffer(1e-7)) + properties = normalized["properties"] + assert properties["source_name"] == "dov_soil_map" + assert properties["soil_texture_class"] == "zand" + assert properties["soil_drainage_class"] == "nat" + assert properties["survey_period"] == "1949-1971" + assert properties["clipped_area_ha"] > 0 + assert "may differ today" in properties["historical_drainage_limitation"] + + +def test_soil_map_uses_existing_semantic_selection_architecture() -> None: + dataset = Dataset( + id=uuid4(), + project_id=uuid4(), + name="dov_soil_map_mol.geojson", + dataset_type="vector", + source="operator_official_import", + source_name="dov_soil_map", + reference_layer_name="soil", + source_metadata={ + "theme": "soil", + "selection_aggregation": { + "method": "intersection_area", + "label": "Bodemkaartoppervlakte", + "unit": "ha", + }, + }, + status="ready", + ) + + assert VectorFeatureService._dataset_theme(dataset) == "soil" + assert VectorFeatureService.supports_selection_summary(dataset) is True + assert VectorFeatureService.can_use_full_area_fast_path(dataset, None) is False + + +def test_soil_operator_contract_has_no_direct_persistence_and_is_packaged() -> None: + operator = (ROOT / "scripts" / "provision_mol_soil_map.py").read_text(encoding="utf-8") + dockerfile = (ROOT / "deploy" / "unraid" / "Dockerfile.all-in-one").read_text(encoding="utf-8") + readiness = (ROOT / "scripts" / "run_readiness_check.sh").read_text(encoding="utf-8") + map_workspace = (ROOT / "frontend" / "src" / "components" / "map" / "MapWorkspace.tsx").read_text(encoding="utf-8") + + assert "/datasets/upload" in operator + assert "vector_features" in operator + assert "does not write directly" in " ".join(operator.split()) + assert "SessionLocal" not in operator and "INSERT INTO" not in operator + assert "COPY scripts/provision_mol_soil_map.py" in dockerfile + assert "py_compile scripts/provision_mol_soil_map.py" in readiness + assert "id: 'soil'" in map_workspace + assert "dataset.source_name === 'dov_soil_map'" in map_workspace + + +def test_incomplete_wfs_pagination_fails_closed() -> None: + module = load_operator() + session = FakeSession( + [ + { + "type": "FeatureCollection", + "numberMatched": 2, + "numberReturned": 0, + "features": [], + } + ] + ) + + with pytest.raises(RuntimeError, match="returned 0 of 2"): + list(module.iter_wfs_pages(session, (0.0, 0.0, 1.0, 1.0), page_limit=100, timeout=30)) diff --git a/deploy/unraid/Dockerfile.all-in-one b/deploy/unraid/Dockerfile.all-in-one index 82647336..0b190c90 100644 --- a/deploy/unraid/Dockerfile.all-in-one +++ b/deploy/unraid/Dockerfile.all-in-one @@ -77,6 +77,8 @@ COPY scripts/provision_mol_context_layers.py /app/scripts/provision_mol_context_ COPY scripts/provision_mol_dhmv.py /app/scripts/provision_mol_dhmv.py COPY scripts/provision_mol_flood_hazards.py /app/scripts/provision_mol_flood_hazards.py COPY scripts/provision_regional_flood_hazards.py /app/scripts/provision_regional_flood_hazards.py +COPY scripts/provision_thematic_rasters.py /app/scripts/provision_thematic_rasters.py +COPY scripts/provision_mol_soil_map.py /app/scripts/provision_mol_soil_map.py COPY scripts/provision_mol_population_history.py /app/scripts/provision_mol_population_history.py COPY scripts/provision_mol_historical_landuse.py /app/scripts/provision_mol_historical_landuse.py COPY scripts/provision_regional_historical_landuse.py /app/scripts/provision_regional_historical_landuse.py diff --git a/deploy/unraid/geointel-unraid-template.xml b/deploy/unraid/geointel-unraid-template.xml index 5fb0da1c..62168241 100644 --- a/deploy/unraid/geointel-unraid-template.xml +++ b/deploy/unraid/geointel-unraid-template.xml @@ -43,6 +43,10 @@ https://geoservice.waterinfo.be/OGRK/wcs 5.0 12000000 + true + https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs + 20000 + 12000000 true http://host.docker.internal:11434 qwen3.5:9b diff --git a/deploy/unraid/run-dockerman-container.sh b/deploy/unraid/run-dockerman-container.sh index 6551b747..f70695b8 100644 --- a/deploy/unraid/run-dockerman-container.sh +++ b/deploy/unraid/run-dockerman-container.sh @@ -43,6 +43,13 @@ FLOOD_HAZARD_MAX_SIDE_M="${FLOOD_HAZARD_MAX_SIDE_M:-20000}" FLOOD_HAZARD_MAX_PIXELS="${FLOOD_HAZARD_MAX_PIXELS:-12000000}" FLOOD_HAZARD_TIMEOUT_SECONDS="${FLOOD_HAZARD_TIMEOUT_SECONDS:-300}" FLOOD_HAZARD_MAX_RESPONSE_MB="${FLOOD_HAZARD_MAX_RESPONSE_MB:-160}" +THEMATIC_RASTER_ENABLED="${THEMATIC_RASTER_ENABLED:-true}" +THEMATIC_RASTER_WCS_URL="${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs}" +THEMATIC_RASTER_MIN_SIDE_M="${THEMATIC_RASTER_MIN_SIDE_M:-100}" +THEMATIC_RASTER_MAX_SIDE_M="${THEMATIC_RASTER_MAX_SIDE_M:-20000}" +THEMATIC_RASTER_MAX_PIXELS="${THEMATIC_RASTER_MAX_PIXELS:-12000000}" +THEMATIC_RASTER_TIMEOUT_SECONDS="${THEMATIC_RASTER_TIMEOUT_SECONDS:-300}" +THEMATIC_RASTER_MAX_RESPONSE_MB="${THEMATIC_RASTER_MAX_RESPONSE_MB:-160}" YOLO_ENABLED="${YOLO_ENABLED:-false}" YOLO_MODELS_DIR="${YOLO_MODELS_DIR:-/app/models}" YOLO_MODEL_PATH="${YOLO_MODEL_PATH:-}" @@ -135,6 +142,13 @@ docker run -d \ -e FLOOD_HAZARD_MAX_PIXELS="$FLOOD_HAZARD_MAX_PIXELS" \ -e FLOOD_HAZARD_TIMEOUT_SECONDS="$FLOOD_HAZARD_TIMEOUT_SECONDS" \ -e FLOOD_HAZARD_MAX_RESPONSE_MB="$FLOOD_HAZARD_MAX_RESPONSE_MB" \ + -e THEMATIC_RASTER_ENABLED="$THEMATIC_RASTER_ENABLED" \ + -e THEMATIC_RASTER_WCS_URL="$THEMATIC_RASTER_WCS_URL" \ + -e THEMATIC_RASTER_MIN_SIDE_M="$THEMATIC_RASTER_MIN_SIDE_M" \ + -e THEMATIC_RASTER_MAX_SIDE_M="$THEMATIC_RASTER_MAX_SIDE_M" \ + -e THEMATIC_RASTER_MAX_PIXELS="$THEMATIC_RASTER_MAX_PIXELS" \ + -e THEMATIC_RASTER_TIMEOUT_SECONDS="$THEMATIC_RASTER_TIMEOUT_SECONDS" \ + -e THEMATIC_RASTER_MAX_RESPONSE_MB="$THEMATIC_RASTER_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 6dbc5c51..342eff25 100644 --- a/docker-compose.unraid.yml +++ b/docker-compose.unraid.yml @@ -41,6 +41,13 @@ services: FLOOD_HAZARD_MAX_PIXELS: ${FLOOD_HAZARD_MAX_PIXELS:-12000000} FLOOD_HAZARD_TIMEOUT_SECONDS: ${FLOOD_HAZARD_TIMEOUT_SECONDS:-300} FLOOD_HAZARD_MAX_RESPONSE_MB: ${FLOOD_HAZARD_MAX_RESPONSE_MB:-160} + THEMATIC_RASTER_ENABLED: ${THEMATIC_RASTER_ENABLED:-true} + THEMATIC_RASTER_WCS_URL: ${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs} + THEMATIC_RASTER_MIN_SIDE_M: ${THEMATIC_RASTER_MIN_SIDE_M:-100} + THEMATIC_RASTER_MAX_SIDE_M: ${THEMATIC_RASTER_MAX_SIDE_M:-20000} + THEMATIC_RASTER_MAX_PIXELS: ${THEMATIC_RASTER_MAX_PIXELS:-12000000} + THEMATIC_RASTER_TIMEOUT_SECONDS: ${THEMATIC_RASTER_TIMEOUT_SECONDS:-300} + THEMATIC_RASTER_MAX_RESPONSE_MB: ${THEMATIC_RASTER_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 604105ee..64879ad8 100644 --- a/docker-compose.yml +++ b/docker-compose.yml @@ -46,6 +46,13 @@ services: FLOOD_HAZARD_MAX_PIXELS: ${FLOOD_HAZARD_MAX_PIXELS:-12000000} FLOOD_HAZARD_TIMEOUT_SECONDS: ${FLOOD_HAZARD_TIMEOUT_SECONDS:-300} FLOOD_HAZARD_MAX_RESPONSE_MB: ${FLOOD_HAZARD_MAX_RESPONSE_MB:-160} + THEMATIC_RASTER_ENABLED: ${THEMATIC_RASTER_ENABLED:-true} + THEMATIC_RASTER_WCS_URL: ${THEMATIC_RASTER_WCS_URL:-https://www.mercator.vlaanderen.be/raadpleegdienstenmercatorpubliek/wcs} + THEMATIC_RASTER_MIN_SIDE_M: ${THEMATIC_RASTER_MIN_SIDE_M:-100} + THEMATIC_RASTER_MAX_SIDE_M: ${THEMATIC_RASTER_MAX_SIDE_M:-20000} + THEMATIC_RASTER_MAX_PIXELS: ${THEMATIC_RASTER_MAX_PIXELS:-12000000} + THEMATIC_RASTER_TIMEOUT_SECONDS: ${THEMATIC_RASTER_TIMEOUT_SECONDS:-300} + THEMATIC_RASTER_MAX_RESPONSE_MB: ${THEMATIC_RASTER_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 7f1ea966..9f4f6856 100644 --- a/docs/API_CONTRACTS.md +++ b/docs/API_CONTRACTS.md @@ -322,6 +322,53 @@ Returns a constrained transparent PNG for a persisted governed VMM flood-depth Dataset. It never accepts an arbitrary path or coverage id and is used by the existing MapLibre image-overlay path. +### GET `/api/v1/projects/{project_id}/datasets/thematic-raster/products` + +Returns the fixed MercatorNet registry for `space_occupation_2025`, +`open_space_2022`, `population_density_2019`, `node_value_2022` and +`service_level_2022`. Every item includes the governed WCS coverage id, native +resolution, source unit, observation year, legend, attribution and limitation. + +### POST `/api/v1/projects/{project_id}/datasets/thematic-raster/acquire` + +Acquires one allowlisted official coverage behind the synchronous Job +abstraction. The request accepts only an EPSG:4326 bbox, optional project Area, +one registry product key and an explicit refresh flag: + +```json +{ + "bbox": {"min_x": 5.03, "min_y": 51.15, "max_x": 5.25, "max_y": 51.33, "crs": "EPSG:4326"}, + "area_id": "optional-project-area-uuid", + "product_key": "population_density_2019", + "force_refresh": false +} +``` + +The backend uses native 10 m or 100 m resolution, splits requests into bounded +WCS 1.0 tiles, validates EPSG:31370 and documented source values, masks the +exact Area and persists an ordinary raster Dataset and DatasetVersion. It does +not accept arbitrary URLs, coverage ids, resolutions or expressions. + +### POST `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select` + +Returns source-correct metrics for a bbox and optional exact Area mask: + +- occupied/open hectares, share and valid raster area for binary products; +- estimated inhabitants plus mean/P90 inhabitants per hectare for the 2019 + population raster; +- mean, P10, median and P90 source score for node value and service level. + +The response names estimate status, aggregation method, source unit, +observation year, attribution, unsupported metrics and product limitation. +Current register population, live public-transport availability and causal +interpretations are not produced. + +### GET `/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/image` + +Returns a constrained transparent PNG generated from the persisted governed +raster. It accepts neither an arbitrary file path nor a provider URL and feeds +the existing MapLibre image-overlay path. + ### GET `/api/v1/projects/{project_id}/datasets` List datasets. diff --git a/docs/CODEX_EXECUTION_LOG.md b/docs/CODEX_EXECUTION_LOG.md index d7db7451..73ae4997 100644 --- a/docs/CODEX_EXECUTION_LOG.md +++ b/docs/CODEX_EXECUTION_LOG.md @@ -1,3 +1,29 @@ +## Sprint 213-214 Cross-domain thematic rasters and DOV soil map (2026-07-16) + +Changed: +- Added a governed five-product MercatorNet thematic-raster registry covering + space, people, accessibility and services with fixed coverage ids, native + resolutions, units, years, legends and limitations. +- Added safe WCS tiling, exact Area masking, canonical raster persistence, + MapLibre rendering, selection metrics and persisted-metric Ollama context. +- Added an explicit DOV soil-map operator for Mol with deterministic WFS 2.0 + pagination, exact EPSG:31370 clipping, raw checksums and canonical vector + upload. No source writes directly to PostGIS. +- Added end-user themes for space occupation, open space, population, + accessibility, services and soil. Current-only official states no longer + make Evolution appear available without at least two observations. + +Validated during implementation: +- Focused thematic and soil suites passed with 14 tests. +- Frontend TypeScript typecheck passed after map and API integration. +- Live source contracts were checked against MercatorNet WCS and the DOV + production WFS; live runtime provisioning follows after deployment. + +Next: +- Run the complete readiness gate, deploy Tower, provision all six Mol layers, + verify live PostGIS metrics and inspect the map at desktop, widescreen and + mobile widths before starting regional rollout. + ## Sprint 195 Guided raster-to-detection workflow (2026-07-14) Changed: diff --git a/docs/DATAVINDPLAATS_SOURCE_ROADMAP.md b/docs/DATAVINDPLAATS_SOURCE_ROADMAP.md index 322ab36f..5a1f0c4e 100644 --- a/docs/DATAVINDPLAATS_SOURCE_ROADMAP.md +++ b/docs/DATAVINDPLAATS_SOURCE_ROADMAP.md @@ -28,8 +28,10 @@ CRS/unit validation, persisted provenance and selection-metric tests pass. - Historical context: official 1778, 1873 and 1969 mapped land-use classes and modern land-use editions where persisted. -This is a strong base, but population/services, soil classes and modeled -accessibility are still genuine coverage gaps. +The cross-domain Wave 1 implementation now adds population/services, modeled +accessibility, space occupation, open space and the DOV soil classes to Mol. +Regional completeness beyond Mol remains a controlled operator rollout rather +than an automatic startup fetch. ## Domain A - Space, Buildings And Economy @@ -142,6 +144,13 @@ DatasetService/VectorFeatureService path. This immediately gives every drawn rectangle a balanced profile across space, soil, people, services and accessibility. +Implementation status on 2026-07-16: complete for Mol. The five raster +products use one fixed MercatorNet registry and the soil map uses the official +`bodemkaart:bodemtypes` DOV WFS. All outputs are persisted through existing +Dataset services, carry source year/period and limitations, and are selectable +on the map. Regional raster provisioning is supported per persisted +municipality; regional soil partitioning remains the next scale-out step. + ### Wave 2 - Economic and mobility objects Add governed vector operators for business parks, Hoppinpoints and cycle diff --git a/docs/DATA_SOURCES.md b/docs/DATA_SOURCES.md index fb574461..58f85328 100644 --- a/docs/DATA_SOURCES.md +++ b/docs/DATA_SOURCES.md @@ -26,6 +26,43 @@ aanvraag. GeoIntel verzint geen historische pixelopnamedatum. Bronnen: Dit document verzamelt concrete databronnen voor GeoIntel Kempen. +## Cross-domain official area profile + +GeoIntel uses one governed MercatorNet WCS registry for five non-water policy +rasters. The operator is explicit and bounded; no source fetch occurs during +startup or directly from the browser. + +| Product | Coverage | Resolution | Selection output | +| --- | --- | ---: | --- | +| Ruimtebeslag 2025 | `lu:lu_ruibes_vlaa_2025_v3` | 10 m | occupied hectares and share | +| Open ruimte 2022 | `lu:lu_openruimte_vlaa_2022_v3` | 10 m | open-space hectares and share | +| Inwonersdichtheid 2019 | `ni:ni_inw_ha_vlaa_2019` | 100 m | estimated inhabitants and inhabitants/ha | +| Knooppuntwaarde 2022 | `lu:lu_knptw_ha_2022_v3` | 100 m | source-score distribution | +| Voorzieningenniveau 2022 | `lu:lu_totvznv_ha_2022_v3` | 100 m | normalized source-score distribution | + +All coverages are stored in EPSG:31370 with their native cells, exact Area +mask, request/checksum provenance and reference year. Space occupation is not +the same as buildings or paving. Open space is not automatically nature or +publicly accessible land. Population is a 2019 raster estimate, not a current +register count. Accessibility and service scores are not live travel times or +object counts. + +## DOV digital soil map + +- Service: `https://www.dov.vlaanderen.be/geoserver/wfs` +- Layer: `bodemkaart:bodemtypes` +- Source CRS and metric clipping: EPSG:31370 +- Persisted geometry: EPSG:4326 through DatasetService/VectorFeatureService +- Scale: 1:20,000 +- Survey evidence: field data collected between 1949 and 1971 + +The Mol operator follows all WFS pages, keeps checksummed raw responses and +persists soil type, series, generalized legend, texture, drainage, profile and +substrate fields. Selection output reports intersected mapped hectares and +fixed sand/anthropogenic context classes. The drainage field remains historical +baseline evidence; GeoIntel does not claim it describes current parcel +drainage or replace a site investigation. + ## GRB — Basiskaart Vlaanderen - Naam: Basiskaart Vlaanderen / GRB diff --git a/docs/STORAGE_ARCHITECTURE.md b/docs/STORAGE_ARCHITECTURE.md index d1cb2358..643bf25e 100644 --- a/docs/STORAGE_ARCHITECTURE.md +++ b/docs/STORAGE_ARCHITECTURE.md @@ -205,6 +205,29 @@ The smoke never passes `--apply`. It fails if the cleanup summary is not a dry-run, if any export/file deletion is reported, or if the dry-run candidate fields are missing. +## Governed cross-domain raster and soil evidence + +MercatorNet thematic products use the ordinary raster Dataset and immutable +DatasetVersion paths under `STORAGE_ROOT`. `source_metadata` records product, +coverage, native resolution/unit, render range, reference year and request +bounds. `provenance_metadata` records every tiled WCS URL, transfer checksum, +normalized checksum, exact Area id and validation result. No new storage table +or provider-side path is introduced. + +The DOV Mol soil operator retains evidence under: + +```text +storage/operator-evidence/dov-soil-map/mol/ + raw/dov_soil_map_page_*.json + dov_soil_map_mol.geojson + dov_soil_map_mol.manifest.json +``` + +Only the normalized, exactly clipped GeoJSON enters DatasetService and +`vector_features`. The manifest binds the persisted artifact, Mol boundary, +source page URLs/checksums, feature completeness, class-area summaries and +historical limitations. Raw source responses remain operator evidence. + ## Model storage Model artifacts live under: diff --git a/docs/TODO.md b/docs/TODO.md index 35d81e4d..9fb9b5cb 100644 --- a/docs/TODO.md +++ b/docs/TODO.md @@ -30,6 +30,10 @@ - [x] Add annual agricultural-use parcels through an explicit provider/operator contract. - [x] Extend the official 1778/1873/1969 historical buildings, water and roads series from Mol to the approved regional scope with partitioned source audits. - [x] Connect a drawn rectangle to bounded official orthophoto acquisition, local configured-YOLO detection and persisted GRB QA. +- [x] Add governed official space-occupation, open-space, population-density, accessibility and service-level rasters with semantic map metrics. +- [x] Add the DOV digital soil map for Mol with exact clipping, source attributes and historical survey limitations. +- [ ] Execute the governed thematic-raster operator for every persisted Kempen municipality after the Mol live gate passes. +- [ ] Generalize the DOV soil-map operator to all 28 approved Kempen municipality partitions with one regional snapshot manifest. ## Governed source expansion backlog diff --git a/frontend/README.md b/frontend/README.md index 3a8d4da7..512a06cc 100644 --- a/frontend/README.md +++ b/frontend/README.md @@ -555,6 +555,14 @@ provider technology: space/buildings, nature/agriculture, soil/relief, mobility/accessibility, population/services and climate/living environment. The central definitions live in `src/lib/sourcePortfolio.ts`. +The current Map explorer also recognizes the governed `Ruimtebeslag`, `Open +ruimte`, `Bevolking`, `Bereikbaarheid`, `Voorzieningen` and `Bodem` Datasets. +The first five render persisted raster PNGs with product-specific legends and +return cell-based semantic metrics; `Bodem` renders the persisted DOV polygons +and exposes soil attributes on feature selection. A single official snapshot +never activates Evolution by itself. Evolution is enabled only when at least +two comparable observations exist. + The domain cards count only matching `ready` Datasets as operational. Audited official sources that have not passed acquisition, persistence and metric validation remain inside the collapsed follow-up list with an explicit diff --git a/frontend/src/App.tsx b/frontend/src/App.tsx index 3ea4b771..049011fa 100644 --- a/frontend/src/App.tsx +++ b/frontend/src/App.tsx @@ -191,7 +191,10 @@ function App(): JSX.Element { const floodHazards = datasets.filter( (dataset) => dataset.dataset_type === 'raster' && dataset.source_name === 'vmm_flood_hazard' && dataset.status === 'ready', ) - return [...vectors, ...terrain, ...floodHazards] + const thematicRasters = datasets.filter( + (dataset) => dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster' && dataset.status === 'ready', + ) + return [...vectors, ...terrain, ...floodHazards, ...thematicRasters] }, [datasets], ) diff --git a/frontend/src/components/map/MapWorkspace.tsx b/frontend/src/components/map/MapWorkspace.tsx index 95266b09..48153947 100644 --- a/frontend/src/components/map/MapWorkspace.tsx +++ b/frontend/src/components/map/MapWorkspace.tsx @@ -8,6 +8,7 @@ import { getDatasetDisplayName, getDatasetSourceDisplayName } from '../../lib/da import { TemporalTrendChart } from './TemporalTrendChart' import { terrainImageUrl } from '../../lib/terrainImage' import { floodHazardImageUrl } from '../../lib/floodHazardImage' +import { thematicRasterImageUrl } from '../../lib/thematicRaster' const DEFAULT_SELECTED_FEATURE_FILENAME = 'selected-feature.geojson' const DEFAULT_AREA_SELECTION_FILENAME = 'area-selection.geojson' @@ -15,7 +16,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' | 'flood_hazard' | 'elevation' | 'roads' | 'parcels' +type DataThemeId = 'buildings' | 'space_occupation' | 'open_space' | 'population' | 'forest' | 'nature_value' | 'agriculture' | 'soil' | 'water' | 'flood_hazard' | 'elevation' | 'accessibility' | 'services' | 'roads' | 'parcels' interface DataTheme { id: DataThemeId @@ -39,6 +40,20 @@ const DATA_THEMES: DataTheme[] = [ description: 'Gebouwen en gebouwcontouren uit GRB of een andere persistente bron.', tokens: ['buildings', 'building', 'gebouwen', 'gebouw', 'bebouwing', 'gbg'], }, + { + id: 'space_occupation', + label: 'Ruimtebeslag', + shortLabel: 'Ruimtebeslag', + description: 'Officiële 10 m-beleidskaart van ruimte ingenomen door wonen, economie, infrastructuur en recreatie.', + tokens: ['space_occupation', 'ruimtebeslag', 'ruibes'], + }, + { + id: 'open_space', + label: 'Open ruimte', + shortLabel: 'Open ruimte', + description: 'Officiële 10 m-beleidskaart van open ruimte buiten kernen en ruimtebeslag.', + tokens: ['open_space', 'open ruimte', 'openruimte'], + }, { id: 'population', label: 'Bevolking', @@ -67,6 +82,13 @@ const DATA_THEMES: DataTheme[] = [ description: 'Jaarlijkse officiële landbouwgebruikspercelen en hoofdteeltgroepen.', tokens: ['agriculture', 'agricultural', 'landbouw', 'landbouwgebruik', 'agpa'], }, + { + id: 'soil', + label: 'Bodem', + shortLabel: 'Bodemkaart', + description: 'Historische DOV-bodemkartering met bodemtype, textuur en drainageklasse voor Mol.', + tokens: ['soil', 'bodem', 'bodemkaart', 'bodemtype', 'dov_soil_map'], + }, { id: 'water', label: 'Water', @@ -88,6 +110,20 @@ const DATA_THEMES: DataTheme[] = [ description: 'Maaiveld- of oppervlaktehoogte, reliëf en helling uit DHMV II.', tokens: ['dhmv', 'elevation', 'height', 'hoogte', 'terrain', 'surface', 'dtm', 'dsm', 'reliëf'], }, + { + id: 'accessibility', + label: 'Bereikbaarheid', + shortLabel: 'Knooppuntwaarde', + description: 'Knooppuntwaarde van collectief vervoer per hectare voor referentiejaar 2022.', + tokens: ['accessibility', 'bereikbaarheid', 'knooppuntwaarde', 'knptw'], + }, + { + id: 'services', + label: 'Voorzieningen', + shortLabel: 'Voorzieningenniveau', + description: 'Genormaliseerde nabijheid van basis-, regionale en metropolitane voorzieningen in 2022.', + tokens: ['services', 'voorzieningen', 'voorzieningenniveau', 'totvznv'], + }, { id: 'roads', label: 'Wegen', @@ -106,13 +142,18 @@ const DATA_THEMES: DataTheme[] = [ const DATA_THEME_MAP_STYLES: Record = { buildings: { fill: '#d45f3d', line: '#9f3e24' }, + space_occupation: { fill: '#be3e33', line: '#8f2c24' }, + open_space: { fill: '#267a46', line: '#175c32' }, population: { fill: '#7559a6', line: '#5b3f88' }, forest: { fill: '#347950', line: '#225f3b' }, nature_value: { fill: '#9a4f64', line: '#74364a' }, agriculture: { fill: '#7b8f32', line: '#53671d' }, + soil: { fill: '#9a7040', line: '#6f4c27' }, water: { fill: '#2676a8', line: '#155b85' }, flood_hazard: { fill: '#1597c2', line: '#075985' }, elevation: { fill: '#a57a4b', line: '#315f59' }, + accessibility: { fill: '#0f766e', line: '#115e59' }, + services: { fill: '#b66d16', line: '#854d0e' }, roads: { fill: '#6b7280', line: '#4b5563' }, parcels: { fill: '#a7792f', line: '#7d571f' }, } @@ -126,6 +167,12 @@ function datasetAvailabilityLabel(dataset: DatasetCreateResponse): string { const resolution = Number(dataset.source_metadata?.['analysis_resolution_m']) return `${Number.isFinite(resolution) ? `${resolution.toLocaleString('nl-BE')} m` : 'Raster'} overstromingsscenario` } + if (dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster') { + const resolution = Number(dataset.source_metadata?.['analysis_resolution_m']) + const year = Number(dataset.source_metadata?.['observation_year']) + const resolutionLabel = Number.isFinite(resolution) ? `${resolution.toLocaleString('nl-BE')} m` : 'Raster' + return `${resolutionLabel} officiële bron${Number.isFinite(year) ? ` · ${year}` : ''}` + } return `${(dataset.feature_count ?? dataset.vector_summary?.feature_count ?? 0).toLocaleString('nl-BE')} objecten beschikbaar` } @@ -156,6 +203,12 @@ function datasetMatchesTheme(dataset: DatasetCreateResponse, theme: DataTheme): if (dataset.source_name === 'digitaal_vlaanderen_dhmv') { return theme.id === 'elevation' } + if (dataset.source_name === 'department_omgeving_thematic_raster') { + return dataset.source_metadata?.['theme'] === theme.id + } + if (dataset.source_name === 'dov_soil_map') { + return theme.id === 'soil' + } const searchText = datasetSearchText(dataset) return theme.tokens.some((token) => searchText.includes(token)) } @@ -184,6 +237,7 @@ function pickThemeDataset( (dataset.source_name === 'department_omgeving_land_use' ? 90_000 : 0) + (dataset.source_name === 'inbo_bwk_natura2000' ? 95_000 : 0) + (dataset.source_name === 'agentschap_landbouw_zeevisserij_agricultural_parcels' ? 98_000 : 0) + + (dataset.source_name === 'department_omgeving_thematic_raster' ? 5_000_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_name === 'vmm_flood_hazard' ? 5_000_000 : 0) + @@ -763,7 +817,20 @@ export function MapWorkspace({ opacity: 0.82, } : null - const activeImageOverlay = floodHazardImageOverlay ?? terrainImageOverlay ?? orthophotoImageOverlay + const thematicRasterBounds = activeThemeDataset?.source_name === 'department_omgeving_thematic_raster' + ? activeThemeDataset.source_metadata?.['bbox_epsg4326'] + : null + const thematicRasterImageOverlay = activeThemeDataset?.source_name === 'department_omgeving_thematic_raster' && selectedProjectId && Array.isArray(thematicRasterBounds) && thematicRasterBounds.length === 4 + ? { + url: thematicRasterImageUrl(selectedProjectId, activeThemeDataset.id), + bbox: thematicRasterBounds.map(Number) as [number, number, number, number], + label: getDatasetDisplayName(activeThemeDataset), + opacity: 0.78, + } + : null + const thematicLegendMin = String(activeThemeDataset?.source_metadata?.['legend_min_label'] ?? 'Lagere waarde') + const thematicLegendMax = String(activeThemeDataset?.source_metadata?.['legend_max_label'] ?? 'Hogere waarde') + const activeImageOverlay = thematicRasterImageOverlay ?? floodHazardImageOverlay ?? 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 @@ -775,7 +842,9 @@ export function MapWorkspace({ [availableMapDatasets], ) const activeTemporalSeriesGroups = themeTemporalSeriesMap[activeTheme.id] - const availableEvolutionThemes = DATA_THEMES.filter((theme) => themeTemporalSeriesMap[theme.id].length > 0) + const availableEvolutionThemes = DATA_THEMES.filter((theme) => + themeTemporalSeriesMap[theme.id].some((group) => group.items.length >= 2), + ) const activeTemporalSeriesGroup = activeTemporalSeriesGroups.find((group) => group.key === selectedTemporalSeriesKey) ?? activeTemporalSeriesGroups[0] const activeTemporalSeries = activeTemporalSeriesGroup?.items ?? EMPTY_TEMPORAL_SERIES @@ -815,13 +884,28 @@ export function MapWorkspace({ (metric) => metric.metric_key !== activeSelectionResult?.summary?.primary_metric_key, ) const terrainReliefMetric = activeSupportingMetrics.find((metric) => metric.metric_key === 'relief_m') + const populationDensityMetric = activeSupportingMetrics.find((metric) => metric.metric_key === 'population_density_mean_per_ha') + const scoreMedianMetric = activeSupportingMetrics.find((metric) => metric.metric_key.endsWith('_median')) const activeSecondaryMetric = activeMetricUnit === 'm TAW' ? terrainReliefMetric ? `${terrainReliefMetric.metric_value.toLocaleString('nl-BE', { maximumFractionDigits: 2 })} m reliëf` : null - : selectedAreaSquareMetres && selectedAreaSquareMetres > 0 - ? activeMetricUnit === 'ha' + : activeMetricUnit === 'inwoners' + ? populationDensityMetric ? selectionMetricLabel(populationDensityMetric) : null + : activeMetricUnit.startsWith('score') + ? scoreMedianMetric ? selectionMetricLabel(scoreMedianMetric) : 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 + const activeSecondaryLabel = activeMetricUnit === 'ha' + ? 'Aandeel selectie' + : activeMetricUnit === 'm TAW' + ? 'Reliëf' + : activeMetricUnit === 'inwoners' + ? 'Gemiddelde dichtheid' + : activeMetricUnit.startsWith('score') + ? 'Mediaan' + : 'Dichtheid' const selectedResultProperties = useMemo(() => { const keys = new Map>() for (const feature of activeSelectionResult?.geojson.features ?? []) { @@ -1185,7 +1269,7 @@ export function MapWorkspace({ const dataset = themeDatasetMap[theme.id] const temporalGroups = themeTemporalSeriesMap[theme.id] const temporalGroup = temporalGroups[0] - const evolutionAvailable = temporalGroups.length > 0 + const evolutionAvailable = temporalGroups.some((group) => group.items.length >= 2) const available = Boolean(dataset) && (analysisMode === 'current' || evolutionAvailable) const active = activeThemeId === theme.id const firstObservation = temporalGroup?.items[0]?.observed_at @@ -1382,7 +1466,12 @@ export function MapWorkspace({ />
Werkgebied - {activeImageOverlay ? {activeImageOverlay.label} : null} + {thematicRasterImageOverlay ? ( + + + {thematicLegendMin} → {thematicLegendMax} + + ) : activeImageOverlay ? {activeImageOverlay.label} : null} {analysisOverlayActive ? ( <> AI-kandidaten @@ -1585,7 +1674,7 @@ export function MapWorkspace({ {activeSelectionResult ? resultMetricLabel(activeSelectionResult) : 'Geen resultaat'}
- {activeMetricUnit === 'ha' ? 'Aandeel selectie' : activeMetricUnit === 'm TAW' ? 'Reliëf' : 'Dichtheid'} + {activeSecondaryLabel} {activeSecondaryMetric ?? (selectedDensity === null ? 'n.v.t.' : `${selectedDensity.toLocaleString('nl-BE', { maximumFractionDigits: 1 })} / km2`)}
@@ -1658,7 +1747,7 @@ export function MapWorkspace({ {analysisMode === 'current' ? (
- +
) : null} diff --git a/frontend/src/hooks/useMapSelectionExtract.ts b/frontend/src/hooks/useMapSelectionExtract.ts index e6930313..1afbb3f4 100644 --- a/frontend/src/hooks/useMapSelectionExtract.ts +++ b/frontend/src/hooks/useMapSelectionExtract.ts @@ -4,6 +4,7 @@ import { formatError } from '../lib/formatError' import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types' import { terrainSelectionToMapSelection } from '../lib/terrainSelection' import { floodHazardSelectionToMapSelection } from '../lib/floodHazardSelection' +import { thematicRasterSelectionToMapSelection } from '../lib/thematicRaster' interface MapSelectionExtractOptions { selectedProjectId: string | null @@ -42,8 +43,9 @@ export function useMapSelectionExtract({ } const terrainDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'digitaal_vlaanderen_dhmv' const floodHazardDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'vmm_flood_hazard' - if (!isVectorDatasetType(selectedDataset.dataset_type) && !terrainDataset && !floodHazardDataset) { - setMapSelectionError('Gebiedsanalyse ondersteunt een vectorlaag, DHMV-hoogtemodel of beheerd VMM-overstromingsscenario.') + const thematicRasterDataset = selectedDataset.dataset_type === 'raster' && selectedDataset.source_name === 'department_omgeving_thematic_raster' + if (!isVectorDatasetType(selectedDataset.dataset_type) && !terrainDataset && !floodHazardDataset && !thematicRasterDataset) { + setMapSelectionError('Gebiedsanalyse ondersteunt een vectorlaag of een beheerd thematisch raster.') return null } @@ -63,6 +65,11 @@ export function useMapSelectionExtract({ bbox: { ...bbox, crs: 'EPSG:4326' }, area_id: areaId, })) + : thematicRasterDataset + ? thematicRasterSelectionToMapSelection(await datasetsApi.selectThematicRaster(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, diff --git a/frontend/src/hooks/useMapThemeSelectionInsights.ts b/frontend/src/hooks/useMapThemeSelectionInsights.ts index 0a5cea43..150eb817 100644 --- a/frontend/src/hooks/useMapThemeSelectionInsights.ts +++ b/frontend/src/hooks/useMapThemeSelectionInsights.ts @@ -4,6 +4,7 @@ import { datasetsApi } from '../services/api/datasets' import type { DatasetCreateResponse, VectorSelectionBBox, VectorSelectionResponse } from '../types' import { terrainSelectionToMapSelection } from '../lib/terrainSelection' import { floodHazardSelectionToMapSelection } from '../lib/floodHazardSelection' +import { thematicRasterSelectionToMapSelection } from '../lib/thematicRaster' export interface MapThemeQuery { themeId: TThemeId @@ -66,6 +67,11 @@ export function useMapThemeSelectionInsights( bbox, area_id: areaId, })) + : dataset.dataset_type === 'raster' && dataset.source_name === 'department_omgeving_thematic_raster' + ? thematicRasterSelectionToMapSelection(await datasetsApi.selectThematicRaster(selectedProjectId, dataset.id, { + bbox, + area_id: areaId, + })) : await datasetsApi.selectVectorFeatures(selectedProjectId, dataset.id, { bbox, area_id: areaId, diff --git a/frontend/src/lib/datasetDisplay.ts b/frontend/src/lib/datasetDisplay.ts index bb89ad63..04a610cc 100644 --- a/frontend/src/lib/datasetDisplay.ts +++ b/frontend/src/lib/datasetDisplay.ts @@ -9,8 +9,13 @@ const DATASET_LABEL_BY_LAYER: Record = { forest: 'Bos en groen', nature_value: 'Natuurwaarde', agriculture: 'Landbouwgebruikspercelen', + soil: 'Digitale bodemkaart', flood_hazard: 'Overstromingsgevaar', elevation: 'Hoogte en reliëf', + space_occupation: 'Ruimtebeslag', + open_space: 'Open ruimte', + accessibility: 'Knooppuntwaarde', + services: 'Voorzieningenniveau', building_registry: 'Gebouwenregister', regional_boundary: 'Grens vervoerregio Kempen', municipality_boundaries: 'Gemeentegrenzen Kempen', @@ -29,6 +34,8 @@ const DATASET_SOURCE_LABELS: Record = { vrbg: 'Digitaal Vlaanderen', waterinfo: 'Waterinfo Vlaanderen', vmm_flood_hazard: 'Vlaamse Milieumaatschappij', + department_omgeving_thematic_raster: 'Departement Omgeving', + dov_soil_map: 'Databank Ondergrond Vlaanderen', } export function getDatasetSourceDisplayName(dataset: DatasetCreateResponse): string { @@ -45,6 +52,10 @@ export function getDatasetDisplayName(dataset: DatasetCreateResponse): string { const productName = dataset.source_metadata?.['product_display_name'] return typeof productName === 'string' && productName.trim() ? productName : 'VMM-overstromingsscenario' } + if (dataset.source_name === 'department_omgeving_thematic_raster') { + const productName = dataset.source_metadata?.['product_display_name'] + return typeof productName === 'string' && productName.trim() ? productName : 'Officieel Vlaams themaraster' + } const layer = (dataset.reference_layer_name ?? dataset.source_metadata?.layer_name ?? dataset.source_metadata?.layer_type ?? '') .toString() .toLowerCase() diff --git a/frontend/src/lib/sourcePortfolio.ts b/frontend/src/lib/sourcePortfolio.ts index db41ef45..3569442c 100644 --- a/frontend/src/lib/sourcePortfolio.ts +++ b/frontend/src/lib/sourcePortfolio.ts @@ -35,6 +35,9 @@ const sourceNameIs = (dataset: DatasetCreateResponse, sourceName: string): boole const referenceLayerIs = (dataset: DatasetCreateResponse, layerName: string): boolean => String(dataset.reference_layer_name ?? dataset.source_metadata?.['theme'] ?? '').toLowerCase() === layerName +const thematicProductIs = (dataset: DatasetCreateResponse, productKey: string): boolean => + sourceNameIs(dataset, 'department_omgeving_thematic_raster') && dataset.source_metadata?.['product_key'] === productKey + export const SOURCE_DOMAINS: SourceDomainDefinition[] = [ { key: 'space', @@ -127,7 +130,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ metricExamples: 'hectare ruimtebeslag, aandeel, vergelijking met open ruimte', priority: 'next', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/ruimtebeslag-vlaanderen-toestand-2025', - matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_space_occupation'), + matches: (dataset) => thematicProductIs(dataset, 'space_occupation_2025'), }, { key: 'settlement_typology', @@ -137,7 +140,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ coverage: 'Toestand 2022', value: 'Begrijpelijke morfologie voor kern, lint en verspreide bebouwing.', metricExamples: 'oppervlakte en aandeel per morfologisch type', - priority: 'planned', + priority: 'next', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/kernen-linten-verspreide-bebouwing-in-vlaanderen-kernen-toestand-2022', matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_settlement_typology'), }, @@ -187,7 +190,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ metricExamples: 'hectare open ruimte, aandeel en fragmentatie', priority: 'planned', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/open-ruimte-vlaanderen-toestand-2022', - matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_open_space'), + matches: (dataset) => thematicProductIs(dataset, 'open_space_2022'), }, { key: 'dhmv', @@ -244,10 +247,10 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ owner: 'Departement Omgeving', coverage: '1 ha raster, toestand 2022', value: 'Modelmatige bereikbaarheid via collectief vervoer.', - metricExamples: 'gemiddelde score en aandeel per bereikbaarheidsklasse', + metricExamples: 'gemiddelde, mediaan en percentielen van de bronindex', priority: 'next', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/knooppuntwaarde-per-ha-toestand-2022', - matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_node_value'), + matches: (dataset) => thematicProductIs(dataset, 'node_value_2022'), }, { key: 'hoppin', @@ -295,7 +298,7 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ metricExamples: 'geschat aantal inwoners en inwoners per hectare', priority: 'next', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/inwonersdichtheid-per-ha-vlaanderen-toestand-2019', - matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_population_density'), + matches: (dataset) => thematicProductIs(dataset, 'population_density_2019'), }, { key: 'service_level', @@ -304,10 +307,10 @@ export const OFFICIAL_SOURCE_PORTFOLIO: OfficialSourceDefinition[] = [ owner: 'Departement Omgeving', coverage: '1 ha raster, toestand 2022', value: 'Nabijheid van dagelijkse voorzieningen in één brongetrouwe score.', - metricExamples: 'gemiddelde score en aandeel per voorzieningsklasse', + metricExamples: 'gemiddelde, mediaan en percentielen van de 0-1-score', priority: 'next', url: 'https://www.vlaanderen.be/datavindplaats/catalogus/totaal-voorzieningenniveau-toestand-2022', - matches: (dataset) => sourceNameIs(dataset, 'department_omgeving_service_level'), + matches: (dataset) => thematicProductIs(dataset, 'service_level_2022'), }, { key: 'municipality_indicators', diff --git a/frontend/src/lib/thematicRaster.ts b/frontend/src/lib/thematicRaster.ts new file mode 100644 index 00000000..fe1af233 --- /dev/null +++ b/frontend/src/lib/thematicRaster.ts @@ -0,0 +1,27 @@ +import type { ThematicRasterSelectionResponse, VectorSelectionResponse } from '../types' + +export function thematicRasterImageUrl(projectId: string, datasetId: string): string { + return `/api/v1/projects/${projectId}/datasets/${datasetId}/raster/thematic/image` +} + +export function thematicRasterSelectionToMapSelection(result: ThematicRasterSelectionResponse): VectorSelectionResponse { + return { + selection_bbox: result.selection_bbox, + selection_area_id: result.selection_area_id, + feature_count: result.valid_cell_count, + total_feature_count: result.valid_cell_count, + limit: 0, + truncated: false, + geojson: { type: 'FeatureCollection', features: [] }, + summary: { + ...result.summary, + feature_count: result.valid_cell_count, + is_estimate: true, + warning: result.limitation_message, + metrics: result.summary.metrics.map((metric) => ({ + ...metric, + is_estimate: metric.is_estimate ?? true, + })), + }, + } +} diff --git a/frontend/src/services/api/datasets.ts b/frontend/src/services/api/datasets.ts index da030fba..96c54c04 100644 --- a/frontend/src/services/api/datasets.ts +++ b/frontend/src/services/api/datasets.ts @@ -24,6 +24,9 @@ import type { FloodHazardSelectionResponse, DhmvProductRead, TerrainSelectionResponse, + ThematicRasterAcquireRequest, + ThematicRasterProductRead, + ThematicRasterSelectionResponse, } from '../../types' const DATASET_PAGE_SIZE = 200 @@ -141,6 +144,16 @@ export const datasetsApi = { payload: { bbox: VectorSelectionRequest['bbox']; area_id?: string }, ): Promise => apiPost(`/api/v1/projects/${projectId}/datasets/${datasetId}/raster/flood-hazard/select`, payload), + acquireThematicRaster: (projectId: string, payload: ThematicRasterAcquireRequest): Promise => + apiPost(`/api/v1/projects/${projectId}/datasets/thematic-raster/acquire`, payload), + listThematicRasterProducts: (projectId: string): Promise<{ items: ThematicRasterProductRead[]; total: number }> => + apiGet<{ items: ThematicRasterProductRead[]; total: number }>(`/api/v1/projects/${projectId}/datasets/thematic-raster/products`), + selectThematicRaster: ( + projectId: string, + datasetId: string, + payload: { bbox: VectorSelectionRequest['bbox']; area_id?: string }, + ): Promise => + apiPost(`/api/v1/projects/${projectId}/datasets/${datasetId}/raster/thematic/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 91fffe11..6526984a 100644 --- a/frontend/src/styles/app.css +++ b/frontend/src/styles/app.css @@ -5751,13 +5751,18 @@ section { } .geo-theme-symbol-buildings { background: #d45f3d; } +.geo-theme-symbol-space_occupation { background: #be3e33; } +.geo-theme-symbol-open_space { background: #267a46; } .geo-theme-symbol-population { background: #7559a6; } .geo-theme-symbol-forest { background: #347950; } .geo-theme-symbol-nature_value { background: #9a4f64; } .geo-theme-symbol-agriculture { background: #7b8f32; } +.geo-theme-symbol-soil { background: #9a7040; } .geo-theme-symbol-water { background: #2676a8; } .geo-theme-symbol-flood_hazard { background: #1597c2; } .geo-theme-symbol-elevation { background: #a57a4b; } +.geo-theme-symbol-accessibility { background: #0f766e; } +.geo-theme-symbol-services { background: #b66d16; } .geo-theme-symbol-roads { background: #6b7280; } .geo-theme-symbol-parcels { background: #a7792f; } @@ -5932,6 +5937,16 @@ section { background: rgba(117, 89, 166, 0.24); } +.geo-map-legend .geo-legend-layer-space_occupation { + border-color: #8f2c24; + background: rgba(190, 62, 51, 0.28); +} + +.geo-map-legend .geo-legend-layer-open_space { + border-color: #175c32; + background: rgba(38, 122, 70, 0.26); +} + .geo-map-legend .geo-legend-layer-forest { border-color: #225f3b; background: rgba(52, 121, 80, 0.24); @@ -5942,6 +5957,11 @@ section { background: rgba(38, 118, 168, 0.24); } +.geo-map-legend .geo-legend-layer-soil { + border-color: #6f4c27; + background: rgba(154, 112, 64, 0.24); +} + .geo-map-legend .geo-legend-layer-flood_hazard { border-color: #075985; background: rgba(21, 151, 194, 0.28); @@ -5952,6 +5972,16 @@ section { background: rgba(165, 122, 75, 0.26); } +.geo-map-legend .geo-legend-layer-accessibility { + border-color: #115e59; + background: rgba(15, 118, 110, 0.26); +} + +.geo-map-legend .geo-legend-layer-services { + border-color: #854d0e; + background: rgba(182, 109, 22, 0.26); +} + .geo-map-legend .geo-legend-layer-roads { border-color: #4b5563; background: rgba(107, 114, 128, 0.24); @@ -5972,6 +6002,32 @@ section { background: linear-gradient(135deg, #7a9b68 0 33%, #d1b37a 33% 66%, #8eb5cb 66%); } +.geo-map-legend .geo-legend-thematic { + display: inline-flex; + align-items: center; + gap: 0.42rem; +} + +.geo-map-legend .geo-legend-thematic small { + color: #42504c; + font-size: 0.66rem; + font-weight: 700; +} + +.geo-map-legend .geo-legend-ramp { + display: block; + width: 3.4rem; + height: 0.62rem; + border: 1px solid rgba(23, 39, 34, 0.25); + border-radius: 2px; +} + +.geo-legend-ramp-space_occupation { background: linear-gradient(90deg, #fbe7d3, #be3e33); } +.geo-legend-ramp-open_space { background: linear-gradient(90deg, #ddeedb, #267a46); } +.geo-legend-ramp-population { background: linear-gradient(90deg, #eee7f6, #673a97); } +.geo-legend-ramp-accessibility { background: linear-gradient(90deg, #e9f1f4, #0f766e); } +.geo-legend-ramp-services { background: linear-gradient(90deg, #fff4bf, #b66d16); } + .geo-map-legend .geo-legend-added { border-color: #15803d; background: rgba(22, 163, 74, 0.2); diff --git a/frontend/src/types.ts b/frontend/src/types.ts index 71536518..db6d3499 100644 --- a/frontend/src/types.ts +++ b/frontend/src/types.ts @@ -434,6 +434,57 @@ export interface FloodHazardSelectionResponse { generated_at: string } +export interface ThematicRasterAcquireRequest { + bbox: VectorSelectionBBox + area_id?: string | null + product_key: string + force_refresh?: boolean +} + +export interface ThematicRasterProductRead { + key: string + display_name: string + theme: 'space_occupation' | 'open_space' | 'population' | 'accessibility' | 'services' + metric_kind: 'binary_area' | 'population_density' | 'index_score' | 'normalized_score' + coverage_id: string + native_resolution_m: number + source_crs: 'EPSG:31370' + source_value_unit: string + observation_year: number + source_version: string + catalog_url: string + attribution: string + license_note: string + legend_min_label: string + legend_max_label: string + limitation_message: string +} + +export interface ThematicRasterSelectionResponse { + dataset_id: string + product_key: string + theme: ThematicRasterProductRead['theme'] + metric_kind: ThematicRasterProductRead['metric_kind'] + selection_bbox: VectorSelectionBBox + selection_area_id?: string | null + selected_cell_count: number + valid_cell_count: number + coverage_ratio: number + resolution_m: number + observation_year: number + summary: { + metric_label: string + metric_value: number + metric_unit: string + aggregation_method: string + primary_metric_key: string + metrics: Array + } + 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 b907c718..50534b6b 100644 --- a/scripts/README.md +++ b/scripts/README.md @@ -1613,6 +1613,33 @@ take a long time because every VMM WCS tile is bounded, rate-limited and validated. This is expected operator work; the app never fetches these rasters on page load or map click. +## Cross-domain Mol profile + +Load the five official policy rasters for the exact Mol municipality Area and +immediately verify each persisted selection result: + +```bash +docker exec geointel python /app/scripts/provision_thematic_rasters.py +``` + +Plan the later complete Kempen rollout without source fetches or writes: + +```bash +docker exec geointel python /app/scripts/provision_thematic_rasters.py \ + --project-name "Kempen Regional Workbench" --all-municipalities --dry-run +``` + +Load the official DOV soil map for Mol: + +```bash +docker exec geointel python /app/scripts/provision_mol_soil_map.py +``` + +`--fetch-only` builds the soil artifact and manifest without API import; +`--force` is the only way to bypass an existing ready soil Dataset. Both +operators use canonical APIs and persistent operator-evidence storage. They do +not run on application startup. + ## Tower deployment Push the local branch to Gitea, then rebuild the Unraid/Tower Docker runtime: diff --git a/scripts/provision_mol_soil_map.py b/scripts/provision_mol_soil_map.py new file mode 100644 index 00000000..a738920d --- /dev/null +++ b/scripts/provision_mol_soil_map.py @@ -0,0 +1,661 @@ +"""Provision the official DOV digital soil map for the municipality of Mol. + +The operator follows every bounded WFS page, retains checksummed source +responses, clips soil polygons to the persisted Mol Area in EPSG:31370 and +imports the result through GeoIntel's canonical dataset upload route. It does +not write directly to vector_features and it does not treat the historical +1949-1971 field survey as a current drainage observation. +""" + +from __future__ import annotations + +import argparse +import hashlib +import json +import os +import sys +from collections import defaultdict +from datetime import datetime, timezone +from pathlib import Path +from typing import Any, Iterable + +import requests +from pyproj import Transformer +from requests.adapters import HTTPAdapter +from shapely.geometry import MultiPolygon, Polygon, mapping, shape +from shapely.ops import transform as transform_geometry +from shapely.ops import unary_union +from shapely.validation import make_valid +from urllib3.util.retry import Retry + + +WFS_URL = "https://www.dov.vlaanderen.be/geoserver/wfs" +CATALOG_URL = ( + "https://www.vlaanderen.be/datavindplaats/catalogus/" + "digitale-bodemkaart-van-het-vlaams-gewest-bodemtypes" +) +TYPE_NAME = "bodemkaart:bodemtypes" +SOURCE_NAME = "dov_soil_map" +SOURCE_VERSION = "Digitale uitgave juni 2017" +SURVEY_PERIOD = "1949-1971" +OBSERVED_AT = "1971-12-31T23:59:59Z" +VALID_FROM = "1949-01-01T00:00:00Z" +VALID_TO = OBSERVED_AT +ATTRIBUTION = "Databank Ondergrond Vlaanderen - Digitale bodemkaart: bodemtypes" +DEFAULT_API_URL = "http://127.0.0.1:8000" +DEFAULT_OUTPUT_DIR = "/app/storage/operator-evidence/dov-soil-map/mol" +DEFAULT_PROJECT_NAME = "Mol Municipality Workbench" +DEFAULT_AREA_FRAGMENT = "Gemeente Mol" +DATASET_FILENAME = "dov_soil_map_mol.geojson" +MANIFEST_FILENAME = "dov_soil_map_mol.manifest.json" +SCHEMA_VERSION = 1 + +TO_LAMBERT72 = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True) +TO_WGS84 = Transformer.from_crs("EPSG:31370", "EPSG:4326", always_xy=True) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Provision the official DOV soil map 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-fragment", default=DEFAULT_AREA_FRAGMENT) + parser.add_argument( + "--output-dir", + type=Path, + default=Path(os.environ.get("GEOINTEL_SOIL_MAP_OUTPUT_DIR", DEFAULT_OUTPUT_DIR)), + ) + parser.add_argument("--page-limit", type=int, default=500) + parser.add_argument("--max-features", type=int, default=20_000) + parser.add_argument("--request-timeout", type=int, default=180) + parser.add_argument("--import-timeout", type=int, default=1800) + parser.add_argument("--force", action="store_true") + parser.add_argument("--fetch-only", action="store_true") + return parser.parse_args() + + +def utc_now() -> str: + return datetime.now(timezone.utc).isoformat() + + +def sha256_bytes(value: bytes) -> str: + return hashlib.sha256(value).hexdigest() + + +def sha256_file(path: Path) -> str: + digest = hashlib.sha256() + with path.open("rb") as handle: + for chunk in iter(lambda: handle.read(1024 * 1024), b""): + digest.update(chunk) + return digest.hexdigest() + + +def write_bytes_atomic(path: Path, value: bytes) -> None: + path.parent.mkdir(parents=True, exist_ok=True) + temporary = path.with_suffix(path.suffix + ".tmp") + temporary.write_bytes(value) + temporary.replace(path) + + +def write_json_atomic(path: Path, value: Any, *, pretty: bool = False) -> None: + encoded = json.dumps( + value, + ensure_ascii=False, + indent=2 if pretty else None, + separators=None if pretty else (",", ":"), + ).encode("utf-8") + write_bytes_atomic(path, encoded) + + +def source_session() -> requests.Session: + retry = Retry( + total=5, + connect=5, + read=5, + status=5, + backoff_factor=1.0, + status_forcelist=(429, 500, 502, 503, 504), + allowed_methods=frozenset({"GET"}), + raise_on_status=True, + ) + session = requests.Session() + session.headers.update({"User-Agent": "GeoIntel-DOV-Soil-Mol-Operator/1.0"}) + adapter = HTTPAdapter(max_retries=retry) + session.mount("https://", adapter) + session.mount("http://", adapter) + return session + + +def response_data(response: requests.Response) -> Any: + try: + payload = response.json() + except ValueError as exc: + raise RuntimeError(f"GeoIntel API returned non-JSON ({response.status_code}): {response.text[:300]}") from exc + if not response.ok: + raise RuntimeError( + f"GeoIntel API failed ({response.status_code}): {json.dumps(payload, ensure_ascii=False)[:800]}" + ) + if not isinstance(payload, dict) or "data" not in payload: + raise RuntimeError("GeoIntel API response does not use the canonical data envelope") + return payload["data"] + + +def paginated_items(session: requests.Session, url: str, *, timeout: int) -> list[dict[str, Any]]: + items: list[dict[str, Any]] = [] + offset = 0 + total: int | None = None + while total is None or offset < total: + page = response_data(session.get(url, params={"limit": 200, "offset": offset}, timeout=timeout)) + page_items = list(page.get("items") or []) + page_total = int(page.get("total") or 0) + if total is None: + total = page_total + elif page_total != total: + raise RuntimeError("GeoIntel pagination total changed while locating the Mol workspace") + items.extend(page_items) + if not page_items: + break + offset += len(page_items) + if total is not None and len(items) != total: + raise RuntimeError(f"GeoIntel pagination returned {len(items)} of {total} items") + return items + + +def polygonal_geometry(geometry): + if geometry is None or geometry.is_empty: + return None + if not geometry.is_valid: + geometry = make_valid(geometry) + polygons: list[Polygon] = [] + + def collect(candidate) -> None: + if candidate is None or candidate.is_empty: + return + if isinstance(candidate, Polygon): + polygons.append(candidate) + elif isinstance(candidate, MultiPolygon): + polygons.extend(part for part in candidate.geoms if not part.is_empty) + elif hasattr(candidate, "geoms"): + for part in candidate.geoms: + collect(part) + + collect(geometry) + if not polygons: + return None + result = unary_union(polygons) + if not result.is_valid: + result = make_valid(result) + return result if not result.is_empty and result.is_valid else None + + +def locate_workspace( + session: requests.Session, + base_url: str, + project_name: str, + area_fragment: str, + timeout: int, +) -> tuple[str, str, Any, list[dict[str, Any]]]: + projects = paginated_items(session, f"{base_url}/api/v1/projects", timeout=timeout) + project = next((item for item in projects if item.get("name") == project_name), None) + if not project: + raise RuntimeError(f"Project {project_name!r} is missing") + project_id = str(project["id"]) + areas = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", timeout=timeout) + area = next( + (item for item in areas if area_fragment.casefold() in str(item.get("name") or "").casefold()), + None, + ) + if not area or not area.get("geometry"): + raise RuntimeError(f"Persisted Mol Area containing {area_fragment!r} is missing") + boundary_wgs84 = polygonal_geometry(shape(area["geometry"])) + if boundary_wgs84 is None: + raise RuntimeError("Persisted Mol Area is not valid polygonal geometry") + datasets = paginated_items(session, f"{base_url}/api/v1/projects/{project_id}/datasets", timeout=timeout) + return project_id, str(area["id"]), boundary_wgs84, datasets + + +def iter_wfs_pages( + session: requests.Session, + bbox_lambert72: tuple[float, float, float, float], + *, + page_limit: int, + timeout: int, +) -> Iterable[tuple[dict[str, Any], str, bytes]]: + start_index = 0 + expected_total: int | None = None + while True: + params = { + "service": "WFS", + "version": "2.0.0", + "request": "GetFeature", + "typeNames": TYPE_NAME, + "srsName": "EPSG:4326", + "bbox": ",".join(f"{value:.3f}" for value in bbox_lambert72) + ",EPSG:31370", + "count": str(page_limit), + "startIndex": str(start_index), + "sortBy": "gid", + "outputFormat": "application/json", + } + response = session.get(WFS_URL, params=params, timeout=timeout) + response.raise_for_status() + payload = response.json() + if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection": + raise RuntimeError("DOV WFS returned an invalid FeatureCollection") + features = list(payload.get("features") or []) + matched = int(payload.get("numberMatched") or payload.get("totalFeatures") or 0) + if expected_total is None: + expected_total = matched + elif matched != expected_total: + raise RuntimeError("DOV WFS numberMatched changed during pagination") + yield payload, response.url, response.content + returned = int(payload.get("numberReturned") or len(features)) + if returned != len(features): + raise RuntimeError("DOV WFS numberReturned does not match its feature payload") + start_index += returned + if returned == 0 or start_index >= expected_total: + if start_index != expected_total: + raise RuntimeError(f"DOV WFS returned {start_index} of {expected_total} matched features") + break + + +def normalize_feature(feature: dict[str, Any], boundary_lambert72) -> tuple[dict[str, Any] | None, bool]: + geometry_payload = feature.get("geometry") + if not geometry_payload: + return None, False + source_wgs84 = polygonal_geometry(shape(geometry_payload)) + if source_wgs84 is None: + return None, False + source_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, source_wgs84)) + if source_lambert72 is None or not source_lambert72.intersects(boundary_lambert72): + return None, False + was_clipped = not source_lambert72.within(boundary_lambert72) + clipped_lambert72 = polygonal_geometry(source_lambert72.intersection(boundary_lambert72)) + if clipped_lambert72 is None or clipped_lambert72.area <= 0: + return None, was_clipped + clipped_wgs84 = polygonal_geometry(transform_geometry(TO_WGS84.transform, clipped_lambert72)) + if clipped_wgs84 is None: + return None, was_clipped + + raw = dict(feature.get("properties") or {}) + gid = raw.get("gid") + map_polygon_id = raw.get("id_kaartvlak") + source_id = str(feature.get("id") or f"{TYPE_NAME}:{gid or map_polygon_id}") + properties = { + "source_name": SOURCE_NAME, + "source_collection": TYPE_NAME, + "source_feature_id": source_id, + "source_gid": gid, + "source_map_polygon_id": map_polygon_id, + "reference_layer_name": "soil", + "theme": "soil", + "authority_level": "authoritative_historical_baseline", + "coverage_scope": "municipality", + "municipality": "Mol", + "nis_code": "13025", + "source_version": SOURCE_VERSION, + "survey_period": SURVEY_PERIOD, + "soil_type_code": raw.get("Bodemtype"), + "unified_soil_type_code": raw.get("Unibodemtype"), + "soil_series_code": raw.get("Bodemserie"), + "soil_series_description": raw.get("Beknopte_omschrijving_bodemserie"), + "soil_generalized_legend": raw.get("Gegeneraliseerde_legende"), + "soil_texture_class_code": raw.get("Textuurklasse_code"), + "soil_texture_class": raw.get("Textuurklasse"), + "soil_drainage_class_code": raw.get("Drainageklasse_code"), + "soil_drainage_class": raw.get("Drainageklasse"), + "soil_profile_group_code": raw.get("Profielontwikkelingsgroep_code"), + "soil_profile_group": raw.get("Profielontwikkelingsgroep"), + "soil_substrate_code": raw.get("Substraat_code"), + "soil_substrate": raw.get("Substraat_Vlaanderen") or raw.get("Substraat_legende"), + "soil_region": raw.get("Streek"), + "classification_type": raw.get("Type_classificatie"), + "soil_map_title": raw.get("Eenduidige_legende_titel"), + "clipped_area_ha": round(float(clipped_lambert72.area) / 10_000.0, 8), + "attribution": ATTRIBUTION, + "historical_drainage_limitation": ( + "Drainage class derives from field data collected between 1949 and 1971 and may differ today." + ), + } + return { + "type": "Feature", + "id": source_id, + "geometry": mapping(clipped_wgs84), + "properties": properties, + }, was_clipped + + +def prepare_artifact( + session: requests.Session, + boundary_wgs84, + output_dir: Path, + *, + page_limit: int, + max_features: int, + timeout: int, +) -> tuple[Path, Path, dict[str, Any]]: + output_dir.mkdir(parents=True, exist_ok=True) + raw_dir = output_dir / "raw" + raw_dir.mkdir(parents=True, exist_ok=True) + boundary_lambert72 = polygonal_geometry(transform_geometry(TO_LAMBERT72.transform, boundary_wgs84)) + if boundary_lambert72 is None: + raise RuntimeError("Mol boundary could not be transformed to EPSG:31370") + + retained: list[dict[str, Any]] = [] + raw_pages: list[dict[str, Any]] = [] + source_urls: list[str] = [] + seen_ids: set[str] = set() + raw_feature_count = 0 + duplicate_count = 0 + rejected_count = 0 + clipped_count = 0 + area_by_legend: dict[str, float] = defaultdict(float) + area_by_texture: dict[str, float] = defaultdict(float) + area_by_drainage: dict[str, float] = defaultdict(float) + + for page_number, (payload, source_url, raw_bytes) in enumerate( + iter_wfs_pages( + session, + boundary_lambert72.bounds, + page_limit=page_limit, + timeout=timeout, + ), + start=1, + ): + page_path = raw_dir / f"dov_soil_map_page_{page_number:05d}.json" + write_bytes_atomic(page_path, raw_bytes) + features = list(payload.get("features") or []) + raw_feature_count += len(features) + if raw_feature_count > max_features: + raise RuntimeError( + f"DOV WFS exceeded the {max_features} feature safety limit; refusing a truncated import" + ) + raw_pages.append( + { + "path": str(page_path.relative_to(output_dir)), + "sha256": sha256_bytes(raw_bytes), + "size_bytes": len(raw_bytes), + "feature_count": len(features), + "source_url": source_url, + } + ) + source_urls.append(source_url) + for feature in features: + raw = dict(feature.get("properties") or {}) + source_id = str(feature.get("id") or f"{TYPE_NAME}:{raw.get('gid')}") + if source_id in seen_ids: + duplicate_count += 1 + continue + seen_ids.add(source_id) + normalized, was_clipped = normalize_feature(feature, boundary_lambert72) + if normalized is None: + rejected_count += 1 + continue + if was_clipped: + clipped_count += 1 + retained.append(normalized) + properties = normalized["properties"] + area = float(properties["clipped_area_ha"]) + area_by_legend[str(properties.get("soil_generalized_legend") or "Onbekend")] += area + area_by_texture[str(properties.get("soil_texture_class") or "Onbekend")] += area + area_by_drainage[str(properties.get("soil_drainage_class") or "Onbekend")] += area + + if not retained: + raise RuntimeError("DOV WFS returned no valid soil polygons inside the persisted Mol Area") + generated_at = utc_now() + artifact = { + "type": "FeatureCollection", + "name": "Digitale bodemkaart - Gemeente Mol", + "features": retained, + "source": ATTRIBUTION, + "source_version": SOURCE_VERSION, + "survey_period": SURVEY_PERIOD, + "catalog_url": CATALOG_URL, + "generated_at": generated_at, + } + artifact_path = output_dir / DATASET_FILENAME + write_json_atomic(artifact_path, artifact) + manifest = { + "schema_version": SCHEMA_VERSION, + "source_version": SOURCE_VERSION, + "survey_period": SURVEY_PERIOD, + "source_type_name": TYPE_NAME, + "wfs_url": WFS_URL, + "catalog_url": CATALOG_URL, + "attribution": ATTRIBUTION, + "generated_at": generated_at, + "crs_source_service": "EPSG:31370", + "crs_response_and_persisted": "EPSG:4326", + "crs_clip_and_area_measurement": "EPSG:31370", + "boundary_sha256": sha256_bytes(json.dumps(mapping(boundary_wgs84), sort_keys=True).encode("utf-8")), + "boundary_bbox_wgs84": list(boundary_wgs84.bounds), + "boundary_bbox_epsg31370": list(boundary_lambert72.bounds), + "page_limit": page_limit, + "page_count": len(raw_pages), + "raw_source_feature_count": raw_feature_count, + "feature_count": len(retained), + "duplicate_count": duplicate_count, + "rejected_or_outside_count": rejected_count, + "clipped_feature_count": clipped_count, + "reference_truncated": False, + "raw_pages": raw_pages, + "source_urls": source_urls, + "area_by_generalized_legend_ha": {key: round(value, 6) for key, value in sorted(area_by_legend.items())}, + "area_by_texture_ha": {key: round(value, 6) for key, value in sorted(area_by_texture.items())}, + "area_by_drainage_ha": {key: round(value, 6) for key, value in sorted(area_by_drainage.items())}, + "artifact_path": str(artifact_path), + "artifact_sha256": sha256_file(artifact_path), + "artifact_size_bytes": artifact_path.stat().st_size, + "limitations": [ + "The map is based on field data collected between 1949 and 1971.", + "Current drainage, land use and local soil disturbance may differ from the mapped class.", + "The 1:20,000 source is contextual evidence and not a parcel-scale soil investigation.", + ], + } + manifest_path = output_dir / MANIFEST_FILENAME + write_json_atomic(manifest_path, manifest, pretty=True) + return artifact_path, manifest_path, manifest + + +def selection_metrics() -> list[dict[str, Any]]: + return [ + { + "metric_key": "soil_dry_sand_area", + "method": "intersection_area", + "label": "Gekarteerd als droog zand", + "unit": "ha", + "geometry_dimension": 2, + "filter_property": "soil_generalized_legend", + "filter_values": ["Droog zand", "Zeer droog zand"], + }, + { + "metric_key": "soil_moist_sand_area", + "method": "intersection_area", + "label": "Gekarteerd als vochtig zand", + "unit": "ha", + "geometry_dimension": 2, + "filter_property": "soil_generalized_legend", + "filter_values": ["Vochtig zand"], + }, + { + "metric_key": "soil_wet_sand_area", + "method": "intersection_area", + "label": "Gekarteerd als nat zand", + "unit": "ha", + "geometry_dimension": 2, + "filter_property": "soil_generalized_legend", + "filter_values": ["Nat zand", "Zeer nat zand"], + }, + { + "metric_key": "soil_anthropogenic_area", + "method": "intersection_area", + "label": "Antropogene bodemklasse", + "unit": "ha", + "geometry_dimension": 2, + "filter_property": "soil_generalized_legend", + "filter_values": ["Antropogeen"], + }, + ] + + +def upload_artifact( + session: requests.Session, + *, + base_url: str, + project_id: str, + area_id: str, + artifact_path: Path, + manifest_path: Path, + manifest: dict[str, Any], + timeout: int, +) -> dict[str, Any]: + limitation = ( + "Historische bodemkartering op schaal 1:20.000 op basis van veldwerk 1949-1971; " + "de huidige drainage en lokale bodemtoestand kunnen afwijken." + ) + source_metadata = { + "provider": SOURCE_NAME, + "theme": "soil", + "layer_type": "soil", + "source_collection": TYPE_NAME, + "source_crs": "EPSG:31370", + "persisted_crs": "EPSG:4326", + "authority_level": "authoritative_historical_baseline", + "coverage_scope": "municipality", + "municipality": "Mol", + "nis_code": "13025", + "feature_count": manifest["feature_count"], + "geometry_clipped_to_area": True, + "semantic_metrics": False, + "survey_period": SURVEY_PERIOD, + "source_scale": "1:20,000", + "attribution": ATTRIBUTION, + "catalog_url": CATALOG_URL, + "license_note": "DOV standard attribution and public GDI reuse conditions apply.", + "limitation_message": limitation, + "selection_aggregation": { + "metric_key": "soil_mapped_area", + "method": "intersection_area", + "label": "Bodemkaartoppervlakte", + "unit": "ha", + "geometry_dimension": 2, + "warning": limitation, + }, + "selection_metrics": selection_metrics(), + } + provenance_metadata = { + "operator_tool": "provision_mol_soil_map.py", + "operator_explicit_fetch": True, + "geometry_clipped_to_area": True, + "source_type_name": TYPE_NAME, + "wfs_url": WFS_URL, + "catalog_url": CATALOG_URL, + "manifest_path": str(manifest_path), + "artifact_sha256": manifest["artifact_sha256"], + "raw_page_checksums": {page["path"]: page["sha256"] for page in manifest["raw_pages"]}, + "source_urls": manifest["source_urls"], + "reference_truncated": False, + "generated_at": manifest["generated_at"], + "limitations": manifest["limitations"], + } + with artifact_path.open("rb") as handle: + response = session.post( + f"{base_url}/api/v1/projects/{project_id}/datasets/upload", + data={ + "dataset_type": "vector", + "source": "operator_official_import", + "dataset_role": "reference", + "source_name": SOURCE_NAME, + "reference_layer_name": "soil", + "source_metadata_json": json.dumps(source_metadata, ensure_ascii=False), + "provenance_metadata_json": json.dumps(provenance_metadata, ensure_ascii=False), + "area_id": area_id, + "temporal_series_key": "dov:digital-soil-map:mol", + "observed_at": OBSERVED_AT, + "valid_from": VALID_FROM, + "valid_to": VALID_TO, + "temporal_granularity": "period", + "source_version": SOURCE_VERSION, + }, + files={"file": (artifact_path.name, handle, "application/geo+json")}, + timeout=timeout, + ) + return response_data(response) + + +def main() -> int: + args = parse_args() + if args.page_limit < 1 or args.page_limit > 2000 or args.max_features < args.page_limit: + print(json.dumps({"status": "error", "message": "Invalid page or feature safety limits"}), file=sys.stderr) + return 2 + base_url = args.base_url.rstrip("/") + api_session = requests.Session() + try: + project_id, area_id, boundary, datasets = locate_workspace( + api_session, + base_url, + args.project_name, + args.area_fragment, + args.request_timeout, + ) + existing = next( + ( + item + for item in datasets + if item.get("source_name") == SOURCE_NAME + and str(item.get("area_id") or "") == area_id + and item.get("status") == "ready" + ), + None, + ) + if existing and not args.force: + result = { + "status": "reused", + "project_id": project_id, + "area_id": area_id, + "dataset_id": existing["id"], + "feature_count": existing.get("feature_count"), + } + else: + artifact_path, manifest_path, manifest = prepare_artifact( + source_session(), + boundary, + args.output_dir, + page_limit=args.page_limit, + max_features=args.max_features, + timeout=args.request_timeout, + ) + if args.fetch_only: + result = { + "status": "prepared", + "project_id": project_id, + "area_id": area_id, + "artifact_path": str(artifact_path), + "feature_count": manifest["feature_count"], + } + else: + dataset = upload_artifact( + api_session, + base_url=base_url, + project_id=project_id, + area_id=area_id, + artifact_path=artifact_path, + manifest_path=manifest_path, + manifest=manifest, + timeout=args.import_timeout, + ) + result = { + "status": "imported", + "project_id": project_id, + "area_id": area_id, + "dataset_id": dataset["id"], + "feature_count": dataset.get("feature_count") or manifest["feature_count"], + "artifact_path": str(artifact_path), + } + print(json.dumps(result, ensure_ascii=False, indent=2)) + return 0 + except (OSError, RuntimeError, requests.RequestException, ValueError) as exc: + print(json.dumps({"status": "error", "message": str(exc)}, ensure_ascii=False), file=sys.stderr) + return 1 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/provision_thematic_rasters.py b/scripts/provision_thematic_rasters.py new file mode 100644 index 00000000..0781c3a5 --- /dev/null +++ b/scripts/provision_thematic_rasters.py @@ -0,0 +1,182 @@ +"""Provision governed Flemish thematic rasters through the GeoIntel API. + +The safe default loads all five products for the official Mol municipality +Area. Use --all-municipalities with an explicitly named regional project to +load every persisted municipality Area. The operator never writes to PostGIS +or storage directly and never accepts an arbitrary external service URL. +""" + +from __future__ import annotations + +import argparse +import json +import os +from typing import Any + +import requests + + +DEFAULT_API_URL = "http://127.0.0.1:8000" +DEFAULT_PROJECT_NAME = "Mol Municipality Workbench" +DEFAULT_AREA_FRAGMENT = "Gemeente Mol" +DEFAULT_PRODUCTS = ( + "space_occupation_2025", + "open_space_2022", + "population_density_2019", + "node_value_2022", + "service_level_2022", +) + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description="Provision governed Flemish thematic raster products.") + 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", default=DEFAULT_AREA_FRAGMENT, help="Case-insensitive Area name fragment.") + parser.add_argument("--products", default=",".join(DEFAULT_PRODUCTS)) + parser.add_argument("--all-municipalities", action="store_true", help="Process every Area whose name starts with 'Gemeente '.") + parser.add_argument("--force-refresh", action="store_true") + parser.add_argument("--timeout", type=int, default=900) + parser.add_argument("--dry-run", action="store_true") + return parser.parse_args() + + +def unwrap(response: requests.Response) -> Any: + try: + payload = response.json() + except ValueError as exc: + raise RuntimeError(f"GeoIntel returned non-JSON HTTP {response.status_code}: {response.text[:300]}") from exc + if not response.ok: + error = payload.get("error") if isinstance(payload, dict) else None + message = error.get("message") if isinstance(error, dict) else response.text[:300] + raise RuntimeError(f"GeoIntel HTTP {response.status_code}: {message}") + return payload.get("data") if isinstance(payload, dict) and "data" in payload else payload + + +def paged_items(session: requests.Session, url: str, timeout: int) -> list[dict[str, Any]]: + items: list[dict[str, Any]] = [] + offset = 0 + while True: + separator = "&" if "?" in url else "?" + page = unwrap(session.get(f"{url}{separator}limit=200&offset={offset}", timeout=timeout)) + rows = list(page.get("items") or []) + items.extend(rows) + total = int(page.get("total") or 0) + if not rows or len(items) >= total: + return items + offset += len(rows) + + +def geometry_bbox(geometry: dict[str, Any]) -> dict[str, Any]: + points: list[tuple[float, float]] = [] + + def visit(value: Any) -> None: + if isinstance(value, list) and len(value) >= 2 and all(isinstance(item, (int, float)) for item in value[:2]): + points.append((float(value[0]), float(value[1]))) + return + if isinstance(value, list): + for item in value: + visit(item) + + visit(geometry.get("coordinates")) + if not points: + raise RuntimeError("Persisted Area geometry contains no coordinates") + return { + "min_x": min(point[0] for point in points), + "min_y": min(point[1] for point in points), + "max_x": max(point[0] for point in points), + "max_y": max(point[1] for point in points), + "crs": "EPSG:4326", + } + + +def find_project(projects: list[dict[str, Any]], name: str) -> dict[str, Any]: + matches = [project for project in projects if str(project.get("name", "")).casefold() == name.casefold()] + if len(matches) != 1: + raise RuntimeError(f"Expected exactly one project named {name!r}, found {len(matches)}") + return matches[0] + + +def select_areas(areas: list[dict[str, Any]], fragment: str, all_municipalities: bool) -> list[dict[str, Any]]: + if all_municipalities: + selected = [area for area in areas if str(area.get("name", "")).casefold().startswith("gemeente ")] + else: + selected = [area for area in areas if fragment.casefold() in str(area.get("name", "")).casefold()] + if not selected: + raise RuntimeError("No persisted Area matches the requested scope") + selected.sort(key=lambda item: str(item.get("name", "")).casefold()) + return selected + + +def main() -> int: + args = parse_args() + base_url = args.base_url.rstrip("/") + requested_products = [value.strip() for value in args.products.split(",") if value.strip()] + if not requested_products: + raise RuntimeError("Select at least one thematic raster product") + session = requests.Session() + session.headers.update({"User-Agent": "GeoIntel-Thematic-Raster-Operator/1.0"}) + + projects = paged_items(session, f"{base_url}/api/v1/projects", args.timeout) + project = find_project(projects, args.project_name) + project_id = str(project["id"]) + areas = paged_items(session, f"{base_url}/api/v1/projects/{project_id}/areas", args.timeout) + selected_areas = select_areas(areas, args.area, args.all_municipalities) + registry = unwrap(session.get(f"{base_url}/api/v1/projects/{project_id}/datasets/thematic-raster/products", timeout=args.timeout)) + products = {str(item["key"]): item for item in registry.get("items") or []} + unknown = sorted(set(requested_products) - set(products)) + if unknown: + raise RuntimeError(f"Products are not present in the canonical registry: {', '.join(unknown)}") + + print(json.dumps({ + "status": "planned" if args.dry_run else "running", + "project_id": project_id, + "project_name": project["name"], + "area_count": len(selected_areas), + "products": requested_products, + }, ensure_ascii=False)) + if args.dry_run: + for area in selected_areas: + print(json.dumps({"area_id": area["id"], "area_name": area["name"], "bbox": geometry_bbox(area["geometry"])}, ensure_ascii=False)) + return 0 + + results: list[dict[str, Any]] = [] + for area in selected_areas: + bbox = geometry_bbox(area["geometry"]) + for product_key in requested_products: + acquisition = unwrap(session.post( + f"{base_url}/api/v1/projects/{project_id}/datasets/thematic-raster/acquire", + json={ + "bbox": bbox, + "area_id": area["id"], + "product_key": product_key, + "force_refresh": args.force_refresh, + }, + timeout=args.timeout, + )) + if acquisition.get("status") != "success" or not acquisition.get("output_dataset_id"): + raise RuntimeError(f"Acquisition failed for {area['name']} / {product_key}: {acquisition}") + dataset_id = str(acquisition["output_dataset_id"]) + analysis = unwrap(session.post( + f"{base_url}/api/v1/projects/{project_id}/datasets/{dataset_id}/raster/thematic/select", + json={"bbox": bbox, "area_id": area["id"]}, + timeout=args.timeout, + )) + result = { + "area_id": area["id"], + "area_name": area["name"], + "product_key": product_key, + "dataset_id": dataset_id, + "reused": bool((acquisition.get("result_json") or {}).get("reused")), + "metric": analysis.get("summary"), + "coverage_ratio": analysis.get("coverage_ratio"), + } + results.append(result) + print(json.dumps(result, ensure_ascii=False)) + + print(json.dumps({"status": "complete", "dataset_count": len(results), "project_id": project_id}, ensure_ascii=False)) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/scripts/run_readiness_check.sh b/scripts/run_readiness_check.sh index 6a9a2574..77e2e15c 100755 --- a/scripts/run_readiness_check.sh +++ b/scripts/run_readiness_check.sh @@ -56,6 +56,8 @@ ${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_mol_flood_hazards.py ${PYTHON_BIN} -m py_compile scripts/provision_regional_flood_hazards.py +${PYTHON_BIN} -m py_compile scripts/provision_thematic_rasters.py +${PYTHON_BIN} -m py_compile scripts/provision_mol_soil_map.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