from __future__ import annotations from dataclasses import dataclass from datetime import UTC, datetime import hashlib import io import json import math from pathlib import Path from typing import Any 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, ThematicRasterMetric, ThematicRasterProductRead, ThematicRasterSelectionRequest, ThematicRasterSelectionResponse, ThematicRasterSelectionSummary, WalousAcquisitionResult, ) from app.services.dataset_service import DatasetService @dataclass(frozen=True) class WalousProduct: key: str display_name: str observation_year: int source_filename: str source_version: str catalog_url: str download_url: str source_sha256_filename: str attribution: str accuracy_label: str raw_class_crosswalk: dict[int, int] | None comparability_note: str observation_start: datetime observation_end: datetime class WalousLandCoverService: PROVIDER = "spw_walous_land_cover" SOURCE_CRS = "EPSG:3812" SOURCE_RESOLUTION_M = 1.0 SOURCE_VALUE_UNIT = "walous_class_code" THEME = "land_cover_use" METRIC_KIND = "categorical_area" NODATA = 255 ATTRIBUTION = "Service public de Wallonie (SPW), Aerospacelab S.A." LICENSE_NOTE = ( "CC BY 4.0; cite the official SPW WALOUS edition and identify modifications." ) LIMITATION = ( "GeoIntel analyseert een nearest-neighbour afgeleide van het officiele 1 m WALOUS-raster op de " "geconfigureerde analyseresolutie. Oppervlakten zijn celgebaseerde schattingen; de kaart is landbedekking, " "geen juridisch landgebruik, eigendom, boomtelling of actuele terreinwaarneming." ) # WALOUS has 11 semantic classes, but its official raster codes are not a # continuous 1..11 range. Codes 80 and 90 distinguish low woody cover. CLASS_LABELS = { 1: "Kunstmatige bodembedekking", 2: "Kunstmatige constructies boven maaiveld", 3: "Spoorweg", 4: "Kale bodem", 5: "Oppervlaktewater", 6: "Jaarlijks wisselende kruidlaag", 7: "Jaarronde kruidlaag", 8: "Naaldbomen hoger dan 3 m", 9: "Loofbomen hoger dan 3 m", 80: "Naaldbomen tot 3 m", 90: "Loofbomen tot 3 m", } CLASS_COLORS = { 1: (155, 155, 155), 2: (183, 72, 67), 3: (68, 68, 68), 4: (194, 165, 119), 5: (44, 129, 185), 6: (236, 202, 73), 7: (161, 201, 78), 8: (28, 89, 51), 9: (52, 132, 72), 80: (78, 125, 70), 90: (107, 164, 87), } # The original 2018 product retains stacked two-digit codes. The official # "Classe vue" legend resolves those codes to the visible top class. The # only 2018-only visible class, greenhouses (62), is explicitly normalized # to artificial constructions so the stable 11-class series can be used. WALOUS_2018_CLASS_CROSSWALK = { 0: NODATA, 1: 1, 11: 1, 15: 1, 18: 1, 19: 1, 31: 1, 51: 1, 71: 1, 81: 1, 91: 1, 2: 2, 28: 2, 29: 2, 62: 2, 3: 3, 38: 3, 39: 3, 73: 3, 83: 3, 93: 3, 4: 4, 5: 5, 55: 5, 58: 5, 59: 5, 75: 5, 85: 5, 95: 5, 6: 6, 7: 7, 8: 8, 9: 9, 80: 80, 90: 90, } @staticmethod def _products() -> dict[str, WalousProduct]: products = ( WalousProduct( key="walous_land_cover_2018", display_name="WALOUS landbedekking 2018", observation_year=2018, source_filename="walous_land_cover_2018_3812.tif", source_version="WALOUS_OCS__2018", catalog_url="https://geoportail.wallonie.be/catalogue/a0ad23a1-1845-4bd5-8c2f-0f62d3f1ec75.html", download_url=( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "a0ad23a1-1845-4bd5-8c2f-0f62d3f1ec75/WALOUS_OCS__2018_GEOTIFF_3812.zip" ), source_sha256_filename="walous_land_cover_2018_3812.sha256", attribution="Service public de Wallonie (SPW), UCLouvain, ULB, ISSeP", accuracy_label="Officiele globale nauwkeurigheid 91,5%", raw_class_crosswalk=WalousLandCoverService.WALOUS_2018_CLASS_CROSSWALK, comparability_note=( "De 2018-editie gebruikt een eerdere, deels handmatig geconsolideerde methode. GeoIntel past de " "officiele 'Classe vue'-crosswalk toe en groepeert de 2018-only serreklasse bij constructies; " "trends blijven methodologisch begrensde schattingen." ), observation_start=datetime(2018, 1, 1, tzinfo=UTC), observation_end=datetime(2018, 12, 31, 23, 59, 59, tzinfo=UTC), ), WalousProduct( key="walous_land_cover_2020", display_name="WALOUS landbedekking 2020", observation_year=2020, source_filename="walous_land_cover_2020_3812.tif", source_version="WAL_OCS_IA__2020", catalog_url="https://geoportail.wallonie.be/catalogue/47b348f1-6e7a-4baa-963c-0232a43c0cff.html", download_url=( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "47b348f1-6e7a-4baa-963c-0232a43c0cff/WAL_OCS_IA__2020_GEOTIFF_3812.zip" ), source_sha256_filename="walous_land_cover_2020_3812.sha256", attribution=WalousLandCoverService.ATTRIBUTION, accuracy_label="Officiele globale nauwkeurigheid 83,30%", raw_class_crosswalk=None, comparability_note="", observation_start=datetime(2020, 4, 1, tzinfo=UTC), observation_end=datetime(2020, 4, 24, 23, 59, 59, tzinfo=UTC), ), WalousProduct( key="walous_land_cover_2023", display_name="WALOUS landbedekking 2023", observation_year=2023, source_filename="walous_land_cover_2023_3812.tif", source_version="WAL_OCS_IA__2023", catalog_url="https://geoportail.wallonie.be/catalogue/4e780ba1-463c-478e-95df-d2f1963a150d.html", download_url=( "https://geoservices.wallonie.be/geotraitement/spwdatadownload/results/" "4e780ba1-463c-478e-95df-d2f1963a150d/WAL_OCS_IA__2023_GEOTIFF_3812.zip" ), source_sha256_filename="walous_land_cover_2023_3812.sha256", attribution=WalousLandCoverService.ATTRIBUTION, accuracy_label="Officiele globale nauwkeurigheid 87,10%", raw_class_crosswalk=None, comparability_note="", observation_start=datetime(2023, 5, 27, tzinfo=UTC), observation_end=datetime(2023, 6, 25, 23, 59, 59, tzinfo=UTC), ), ) return {product.key: product for product in products} @staticmethod def _source_path(settings: Settings, product: WalousProduct) -> Path: return Path(settings.walous_source_dir) / product.source_filename @staticmethod def list_products(*, settings: Settings | None = None) -> list[dict[str, Any]]: resolved = settings or get_settings() result: list[dict[str, Any]] = [] for product in WalousLandCoverService._products().values(): configured = ( resolved.walous_enabled and WalousLandCoverService._source_path(resolved, product).is_file() ) result.append( ThematicRasterProductRead( key=product.key, display_name=product.display_name, theme=WalousLandCoverService.THEME, metric_kind=WalousLandCoverService.METRIC_KIND, coverage_id=product.source_version, native_resolution_m=WalousLandCoverService.SOURCE_RESOLUTION_M, analysis_resolution_m=resolved.walous_analysis_resolution_m, source_crs=WalousLandCoverService.SOURCE_CRS, source_value_unit=WalousLandCoverService.SOURCE_VALUE_UNIT, observation_year=product.observation_year, source_version=product.source_version, catalog_url=product.catalog_url, attribution=product.attribution, license_note=WalousLandCoverService.LICENSE_NOTE, legend_min_label="WALOUS klasse 1 (kunstmatige bodem)", legend_max_label="WALOUS klasse 90 (loofbomen tot 3 m)", included_source_values=list(WalousLandCoverService.CLASS_LABELS), limitation_message=" ".join( part for part in ( WalousLandCoverService.LIMITATION, f"{product.accuracy_label}.", product.comparability_note, ) if part ), coverage_zones=["wallonia"], configured=configured, status="configured" if configured else "source_not_provisioned", ).model_dump() ) return result @staticmethod def _product(product_key: str) -> WalousProduct: product = WalousLandCoverService._products().get(product_key.strip().lower()) if product is None: raise AppError( code="WALOUS_PRODUCT_NOT_SUPPORTED", message="Select a product from the governed WALOUS registry", details={"product_key": product_key}, status_code=422, ) return product @staticmethod def _scope_geometry(db, project_id: UUID, payload: ThematicRasterAcquireRequest): if not db.get(Project, project_id): raise AppError( code="PROJECT_NOT_FOUND", message="Project not found", status_code=404 ) if payload.bbox.crs.upper() != "EPSG:4326": raise AppError( code="INVALID_BBOX_CRS", message="WALOUS acquisition requires EPSG:4326", status_code=400, ) values = [ payload.bbox.min_x, payload.bbox.min_y, payload.bbox.max_x, payload.bbox.max_y, ] if ( not all(math.isfinite(value) for value in values) or values[0] >= values[2] or values[1] >= values[3] ): raise AppError( code="INVALID_BBOX", message="WALOUS selection must be a finite non-empty rectangle", status_code=400, ) selection = box(*values) if payload.area_id is None: return selection, values area = db.get(Area, payload.area_id) if area is None or area.project_id != project_id: raise AppError( code="AREA_NOT_FOUND", message="Area not found", status_code=404 ) selection = selection.intersection(to_shape(area.geometry)) if selection.is_empty or selection.area <= 0: raise AppError( code="WALOUS_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422, ) return selection, values @staticmethod def _read_source_window( source_path: Path, scope_4326, settings: Settings, product: WalousProduct, ) -> tuple[bytes, dict[str, Any]]: try: import numpy as np import rasterio from rasterio.enums import Resampling from rasterio.features import geometry_mask from rasterio.io import MemoryFile from rasterio.transform import from_bounds from rasterio.windows import from_bounds as window_from_bounds except ImportError as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio and numpy are required for WALOUS", status_code=503, ) from exc resolution = float(settings.walous_analysis_resolution_m) transformer = Transformer.from_crs( "EPSG:4326", WalousLandCoverService.SOURCE_CRS, always_xy=True ) scope_metric = shapely_transform(transformer.transform, scope_4326) try: with rasterio.open(source_path) as source: if ( source.crs is None or source.crs.to_epsg() != 3812 or source.count != 1 ): raise AppError( code="WALOUS_SOURCE_INVALID", message="WALOUS source must be a one-band EPSG:3812 raster", status_code=409, ) if not all( math.isclose(abs(float(value)), 1.0, abs_tol=0.05) for value in source.res ): raise AppError( code="WALOUS_SOURCE_INVALID", message="WALOUS source must retain the official 1 m resolution", status_code=409, ) clipped_geometry = scope_metric.intersection(box(*source.bounds)) if clipped_geometry.is_empty or clipped_geometry.area <= 0: raise AppError( code="WALOUS_SELECTION_OUTSIDE_COVERAGE", message="Selection does not overlap WALOUS coverage", status_code=422, ) min_x, min_y, max_x, max_y = clipped_geometry.bounds bounds = ( math.floor(min_x / resolution) * resolution, math.floor(min_y / resolution) * resolution, math.ceil(max_x / resolution) * resolution, math.ceil(max_y / resolution) * resolution, ) width_m, height_m = bounds[2] - bounds[0], bounds[3] - bounds[1] if ( width_m > settings.walous_max_side_m or height_m > settings.walous_max_side_m ): raise AppError( code="WALOUS_SELECTION_TOO_LARGE", message=f"Select no more than {settings.walous_max_side_m:g} by {settings.walous_max_side_m:g} metres", details={"width_m": width_m, "height_m": height_m}, status_code=422, ) width, height = ( max(1, round(width_m / resolution)), max(1, round(height_m / resolution)), ) if width * height > settings.walous_max_pixels: raise AppError( code="WALOUS_SELECTION_TOO_LARGE", message="WALOUS selection exceeds the configured cell limit", details={ "pixel_count": width * height, "max_pixels": settings.walous_max_pixels, }, status_code=422, ) window = window_from_bounds(*bounds, transform=source.transform) band = source.read( 1, window=window, out_shape=(height, width), masked=True, resampling=Resampling.nearest, ) output_transform = from_bounds(*bounds, width, height) outside_scope = geometry_mask( [mapping(clipped_geometry)], out_shape=(height, width), transform=output_transform, invert=False, ) # The official 2023 GeoTIFF is signed int8 while GDAL exposes # its nodata sentinel as 255. Filling before widening would # therefore reject the sentinel as out of range for int8. raw = np.asarray(np.ma.getdata(band), dtype="uint8") invalid = np.ma.getmaskarray(band) | outside_scope if source.nodata is not None: invalid |= np.isclose(raw.astype("float64"), float(source.nodata)) raw[invalid] = WalousLandCoverService.NODATA source_valid = raw[raw != WalousLandCoverService.NODATA] if source_valid.size == 0: raise AppError( code="WALOUS_NO_VALID_DATA", message="WALOUS contains no valid cells in this selection", status_code=422, ) source_classes = set(np.unique(source_valid).astype(int).tolist()) governed_source_classes = set( product.raw_class_crosswalk or WalousLandCoverService.CLASS_LABELS ) unexpected = sorted(source_classes - governed_source_classes) if unexpected: raise AppError( code="WALOUS_SOURCE_INVALID_VALUES", message="WALOUS contains classes outside the governed 11-class code set", details={"unexpected_classes": unexpected}, status_code=409, ) if product.raw_class_crosswalk: normalized = np.full( raw.shape, WalousLandCoverService.NODATA, dtype="uint8" ) for ( source_value, normalized_value, ) in product.raw_class_crosswalk.items(): normalized[(raw == source_value) & ~invalid] = normalized_value raw = normalized valid = raw[raw != WalousLandCoverService.NODATA] classes = set(np.unique(valid).astype(int).tolist()) profile = { "driver": "GTiff", "width": width, "height": height, "count": 1, "dtype": "uint8", "crs": WalousLandCoverService.SOURCE_CRS, "transform": output_transform, "nodata": WalousLandCoverService.NODATA, "compress": "deflate", "predictor": 2, } with MemoryFile() as memory: with memory.open(**profile) as output: output.write(raw, 1) content = memory.read() return content, { "width": width, "height": height, "valid_pixel_count": int(valid.size), "classes_present": sorted(classes), "source_classes_present": sorted(source_classes), "class_crosswalk": product.raw_class_crosswalk, "bbox_epsg3812": list(bounds), "source_width": int(source.width), "source_height": int(source.height), "source_nodata": None if source.nodata is None else float(source.nodata), "source_resolution_m": 1.0, "analysis_resolution_m": resolution, } except AppError: raise except Exception as exc: raise AppError( code="WALOUS_SOURCE_READ_FAILED", message="The provisioned WALOUS source could not be read", details={"reason": str(exc)}, status_code=500, ) 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 == WalousLandCoverService.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, ) -> dict[str, Any]: resolved = settings or get_settings() if not resolved.walous_enabled: raise AppError( code="WALOUS_NOT_CONFIGURED", message="WALOUS bounded analysis is disabled", status_code=503, ) product = WalousLandCoverService._product(payload.product_key) source_path = WalousLandCoverService._source_path(resolved, product) if not source_path.is_file(): raise AppError( code="WALOUS_SOURCE_NOT_PROVISIONED", message="The official WALOUS source archive has not been provisioned on this runtime", details={ "expected_path": str(source_path), "operator_command": "python scripts/provision_walous_sources.py --years 2018 2020 2023", }, status_code=503, ) scope, bbox_4326 = WalousLandCoverService._scope_geometry( db, project_id, payload ) identity = { "product_key": product.key, "bbox_epsg4326": [round(float(value), 8) for value in bbox_4326], "area_id": str(payload.area_id) if payload.area_id else None, "analysis_resolution_m": resolved.walous_analysis_resolution_m, } request_hash = hashlib.sha256( json.dumps(identity, sort_keys=True).encode() ).hexdigest() filename = f"walous_{product.observation_year}_{request_hash[:12]}_3812.tif" if not payload.force_refresh: cached = WalousLandCoverService._cached_dataset(db, project_id, filename) if cached is not None: metadata = cached.source_metadata or {} return WalousAcquisitionResult( output_dataset_id=cached.id, reused=True, provider=WalousLandCoverService.PROVIDER, product_key=product.key, display_name=product.display_name, theme=WalousLandCoverService.THEME, metric_kind=WalousLandCoverService.METRIC_KIND, resolution_m=float( metadata.get( "analysis_resolution_m", resolved.walous_analysis_resolution_m, ) ), width=int((cached.metadata_json or {}).get("width", 0)), height=int((cached.metadata_json or {}).get("height", 0)), valid_pixel_count=int(metadata.get("valid_pixel_count", 0)), bbox_epsg4326=bbox_4326, bbox_epsg3812=list(metadata.get("bbox_epsg3812") or []), observation_year=product.observation_year, source_value_unit=WalousLandCoverService.SOURCE_VALUE_UNIT, attribution=product.attribution, limitation_message=" ".join( part for part in ( WalousLandCoverService.LIMITATION, f"{product.accuracy_label}.", product.comparability_note, ) if part ), ).model_dump(mode="json") content, validation = WalousLandCoverService._read_source_window( source_path, scope, resolved, product ) source_sha256_path = source_path.with_name(product.source_sha256_filename) source_sha256 = ( source_sha256_path.read_text(encoding="ascii").strip().split()[0] if source_sha256_path.is_file() else None ) acquired_at = datetime.now(UTC) observed_at = product.observation_end spatial_series_hash = hashlib.sha256( json.dumps( { "bbox": identity["bbox_epsg4326"], "area_id": identity["area_id"], "resolution": identity["analysis_resolution_m"], }, sort_keys=True, ).encode() ).hexdigest()[:24] dataset = DatasetService.import_raster_bytes( db, project_id=project_id, area_id=payload.area_id, filename=filename, content=content, source=f"SPW WALOUS {product.source_version} operator-provisioned GeoTIFF", source_name=WalousLandCoverService.PROVIDER, temporal_series_key=f"spw:walous:land-cover:{spatial_series_hash}", observed_at=observed_at, valid_from=product.observation_start, valid_to=product.observation_end, temporal_granularity="year", source_version=product.source_version, source_metadata={ "provider": WalousLandCoverService.PROVIDER, "service": "official_predefined_dataset_atom", "product_key": product.key, "product_display_name": product.display_name, "theme": WalousLandCoverService.THEME, "metric_kind": WalousLandCoverService.METRIC_KIND, "source_crs": WalousLandCoverService.SOURCE_CRS, "source_resolution_m": WalousLandCoverService.SOURCE_RESOLUTION_M, "analysis_resolution_m": validation["analysis_resolution_m"], "source_value_unit": WalousLandCoverService.SOURCE_VALUE_UNIT, "class_labels": WalousLandCoverService.CLASS_LABELS, "observation_year": product.observation_year, "observation_start": product.observation_start.isoformat(), "observation_end": product.observation_end.isoformat(), "valid_pixel_count": validation["valid_pixel_count"], "classes_present": validation["classes_present"], "source_classes_present": validation["source_classes_present"], "class_crosswalk": validation["class_crosswalk"], "bbox_epsg4326": bbox_4326, "bbox_epsg3812": validation["bbox_epsg3812"], "coverage_zones": ["wallonia"], "catalog_url": product.catalog_url, "download_url": product.download_url, "attribution": product.attribution, "license_note": WalousLandCoverService.LICENSE_NOTE, "limitation_message": " ".join( part for part in ( WalousLandCoverService.LIMITATION, f"{product.accuracy_label}.", product.comparability_note, ) if part ), }, provenance_metadata={ "acquisition": "operator_provisioned_official_archive_bounded_window", "acquired_at": acquired_at.isoformat(), "request_hash": request_hash, "source_filename": product.source_filename, "source_sha256": source_sha256, "derived_sha256": hashlib.sha256(content).hexdigest(), "resampling": "nearest", "validation": validation, }, ) return WalousAcquisitionResult( output_dataset_id=dataset.id, reused=False, provider=WalousLandCoverService.PROVIDER, product_key=product.key, display_name=product.display_name, theme=WalousLandCoverService.THEME, metric_kind=WalousLandCoverService.METRIC_KIND, resolution_m=validation["analysis_resolution_m"], width=validation["width"], height=validation["height"], valid_pixel_count=validation["valid_pixel_count"], bbox_epsg4326=bbox_4326, bbox_epsg3812=validation["bbox_epsg3812"], observation_year=product.observation_year, source_value_unit=WalousLandCoverService.SOURCE_VALUE_UNIT, attribution=product.attribution, limitation_message=" ".join( part for part in ( WalousLandCoverService.LIMITATION, f"{product.accuracy_label}.", product.comparability_note, ) if part ), ).model_dump(mode="json") @staticmethod def _load_dataset( db, project_id: UUID, dataset_id: UUID ) -> tuple[Dataset, WalousProduct]: 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 != WalousLandCoverService.PROVIDER ): raise AppError( code="INVALID_WALOUS_DATASET", message="WALOUS analysis requires a governed WALOUS raster", 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 WALOUS raster is unavailable", status_code=404, ) product = WalousLandCoverService._product( str((dataset.source_metadata or {}).get("product_key") or "") ) return dataset, product @staticmethod def _analysis_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 area is None or area.project_id != project_id: raise AppError( code="AREA_NOT_FOUND", message="Area not found", status_code=404 ) selection = selection.intersection(to_shape(area.geometry)) if selection.is_empty or selection.area <= 0: raise AppError( code="WALOUS_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422, ) return selection @staticmethod def analyze( db, project_id: UUID, dataset_id: UUID, payload: ThematicRasterSelectionRequest ) -> dict[str, Any]: dataset, product = WalousLandCoverService._load_dataset( db, project_id, dataset_id ) selection_4326 = WalousLandCoverService._analysis_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 WALOUS analysis", status_code=503, ) from exc try: with rasterio.open(dataset.storage_path) as source: transformer = Transformer.from_crs( "EPSG:4326", source.crs, always_xy=True ) selection_metric = shapely_transform( transformer.transform, selection_4326 ) geometry = selection_metric.intersection(box(*source.bounds)) if geometry.is_empty or geometry.area <= 0: raise AppError( code="WALOUS_SELECTION_OUTSIDE_DATASET", message="Selection does not overlap the persisted WALOUS raster", status_code=422, ) clipped, transform = mask( source, [mapping(geometry)], crop=True, filled=False, indexes=[1] ) band = np.ma.asarray(clipped[0]) raw = np.asarray(np.ma.getdata(band), dtype="uint8") selected = geometry_mask( [mapping(geometry)], out_shape=raw.shape, transform=transform, invert=True, ) valid = ( selected & ~np.ma.getmaskarray(band) & (raw != WalousLandCoverService.NODATA) ) values = raw[valid] selected_count = int(selected.sum()) valid_count = int(values.size) if not valid_count: raise AppError( code="WALOUS_NO_VALID_DATA", message="WALOUS contains no valid cells in this selection", status_code=422, ) cell_area_m2 = abs(float(source.res[0]) * float(source.res[1])) except AppError: raise except Exception as exc: raise AppError( code="WALOUS_ANALYSIS_FAILED", message="The persisted WALOUS raster could not be analysed", details={"reason": str(exc)}, status_code=500, ) from exc def area_for(classes: set[int]) -> float: return float( np.count_nonzero(np.isin(values, list(classes))) * cell_area_m2 / 10_000.0 ) metric_specs = [ ( "land_cover_observed_area_ha", "Gekarteerde landbedekking", set(WalousLandCoverService.CLASS_LABELS), ), ("forest_cover_area_ha", "Boom- en bosbedekking", {8, 9, 80, 90}), ("surface_water_area_ha", "Oppervlaktewater", {5}), ( "artificial_cover_area_ha", "Kunstmatige bedekking en constructies", {1, 2, 3}, ), ("annual_herbaceous_cover_area_ha", "Jaarlijks wisselende kruidlaag", {6}), ("permanent_herbaceous_cover_area_ha", "Jaarronde kruidlaag", {7}), ("bare_soil_area_ha", "Kale bodem", {4}), ] metrics = [ ThematicRasterMetric( metric_key=key, metric_label=label, metric_value=round(area_for(classes), 4), metric_unit="ha", aggregation_method="nearest_resampled_cells_times_cell_area", is_estimate=True, ) for key, label, classes in metric_specs ] primary = metrics[0] return ThematicRasterSelectionResponse( dataset_id=dataset.id, product_key=product.key, theme=WalousLandCoverService.THEME, metric_kind=WalousLandCoverService.METRIC_KIND, selection_bbox=payload.bbox, selection_area_id=payload.area_id, selected_cell_count=selected_count, valid_cell_count=valid_count, coverage_ratio=round(valid_count / max(1, selected_count), 6), resolution_m=round(math.sqrt(cell_area_m2), 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=[ "legal_land_use", "ownership", "tree_count", "timber_volume", "water_volume", ], limitation_message=" ".join( part for part in ( WalousLandCoverService.LIMITATION, f"{product.accuracy_label}.", product.comparability_note, ) if part ), generated_at=datetime.now(UTC).isoformat(), ).model_dump(mode="json") @staticmethod def render_png( db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1800 ) -> bytes: dataset, _product = WalousLandCoverService._load_dataset( db, project_id, dataset_id ) try: import numpy as np import rasterio from PIL import Image from rasterio.enums import Resampling except ImportError as exc: raise AppError( code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio, numpy and Pillow are required for WALOUS rendering", status_code=503, ) from exc with rasterio.open(dataset.storage_path) as source: scale = min(1.0, max_dimension / max(source.width, source.height)) width, height = ( max(1, round(source.width * scale)), max(1, round(source.height * scale)), ) values = source.read( 1, out_shape=(height, width), masked=True, resampling=Resampling.nearest ) raw = np.asarray(np.ma.getdata(values), dtype="uint8") rgba = np.zeros((height, width, 4), dtype="uint8") for value, color in WalousLandCoverService.CLASS_COLORS.items(): selected = raw == value rgba[:, :, 0][selected] = color[0] rgba[:, :, 1][selected] = color[1] rgba[:, :, 2][selected] = color[2] rgba[:, :, 3][selected] = 205 output = io.BytesIO() Image.fromarray(rgba).save(output, format="PNG", optimize=True) return output.getvalue()