945 lines
39 KiB
Python
945 lines
39 KiB
Python
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()
|