684 lines
37 KiB
Python
684 lines
37 KiB
Python
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.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 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)
|
|
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")
|