geointel/backend/app/services/thematic_raster_acquisition_service.py

from __future__ import annotations

import hashlib
import json
import math
import time
from dataclasses import dataclass
from datetime import UTC, datetime
from http.client import HTTPException
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.services.outbound_request_guard import guarded_opener
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
    included_source_values: tuple[int, ...] = ()


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
    WCS_FETCH_ATTEMPTS = 3
    WCS_RETRY_DELAY_SECONDS = 1.0
    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="forest_land_use_2025",
                display_name="Bos volgens Landgebruik Vlaanderen 2025",
                theme="forest",
                metric_kind="binary_area",
                coverage_id="lu:lu_landgebruik_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/landgebruik-vlaanderen-toestand-2025",
                legend_min_label="Geen bosklasse",
                legend_max_label="Bos",
                limitation_message=(
                    "10 m-afleiding van bronklasse 12 (bos) uit Landgebruik Vlaanderen 2025. De oppervlakte is "
                    "resolutiegebonden en vormt geen juridische bosgrens, boomtelling, kroonbedekking of houtvolume."
                ),
                included_source_values=(12,),
            ),
            ThematicRasterProduct(
                key="agricultural_land_use_2025",
                display_name="Akker en landbouwgrasland 2025",
                theme="agriculture",
                metric_kind="binary_area",
                coverage_id="lu:lu_landgebruik_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/landgebruik-vlaanderen-toestand-2025",
                legend_min_label="Ander landgebruik",
                legend_max_label="Akker of landbouwgrasland",
                limitation_message=(
                    "10 m-afleiding van bronklassen 13 (akker) en 14 (grasland in landbouwgebruik). Dit is werkelijk "
                    "landgebruik en geen ALZ-perceelaangifte, eigendomsgrens, teeltregister of juridische bestemming."
                ),
                included_source_values=(13, 14),
            ),
            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,
                included_source_values=list(product.included_source_values),
                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 _coverage_scope(db, area_id: UUID | None) -> str:
        if area_id is None:
            return "bounded_selection"
        area = db.get(Area, area_id)
        if area and str(area.name).casefold().startswith("gemeente "):
            return "municipality"
        return "bounded_selection"

    @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
        for attempt in range(1, ThematicRasterAcquisitionService.WCS_FETCH_ATTEMPTS + 1):
            try:
                with (opener or guarded_opener(request_url))(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)
                break
            except AppError:
                raise
            except HTTPError as exc:
                preview = exc.read(300).decode("utf-8", errors="replace")
                if attempt < ThematicRasterAcquisitionService.WCS_FETCH_ATTEMPTS and int(exc.code) in {429, 500, 502, 503, 504}:
                    time.sleep(ThematicRasterAcquisitionService.WCS_RETRY_DELAY_SECONDS * attempt)
                    continue
                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, "attempts": attempt},
                    status_code=502,
                ) from exc
            except (URLError, TimeoutError, OSError, HTTPException) as exc:
                if attempt < ThematicRasterAcquisitionService.WCS_FETCH_ATTEMPTS:
                    time.sleep(ThematicRasterAcquisitionService.WCS_RETRY_DELAY_SECONDS * attempt)
                    continue
                raise AppError(
                    code="THEMATIC_RASTER_PROVIDER_UNAVAILABLE",
                    message="The official MercatorNet WCS could not complete the bounded request",
                    details={"reason": str(exc), "attempts": attempt},
                    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))
                source_values = np.ma.array(raw, mask=invalid).compressed().astype("float64")
                if product.included_source_values:
                    rounded = np.rint(source_values)
                    if not np.allclose(source_values, rounded, atol=0.0001):
                        raise AppError(
                            code="THEMATIC_RASTER_INVALID_VALUES",
                            message="Categorical land-use coverage contains non-integer source classes",
                            status_code=502,
                        )
                    if source_values.size and (
                        float(source_values.min()) < 0
                        or float(source_values.max()) > 255
                    ):
                        raise AppError(
                            code="THEMATIC_RASTER_INVALID_VALUES",
                            message="Categorical land-use coverage contains source classes outside the governed range",
                            status_code=502,
                        )
                    source_classes = np.where(invalid, 0, np.rint(raw)).astype("int16")
                    binary = np.isin(source_classes, product.included_source_values).astype("float32")
                    normalized = np.ma.array(binary, mask=invalid)
                else:
                    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)),
                    "included_source_values": list(product.included_source_values),
                    "source_minimum_value": float(source_values.min()),
                    "source_maximum_value": float(source_values.max()),
                }
        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,
                "included_source_values": list(product.included_source_values),
                "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": ThematicRasterAcquisitionService._coverage_scope(db, payload.area_id),
            },
            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")