feat: complete governed Walloon coverage sources
This commit is contained in:
@@ -0,0 +1,551 @@
|
||||
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
|
||||
accuracy_label: str
|
||||
|
||||
|
||||
class WalousLandCoverService:
|
||||
PROVIDER = "spw_walous_land_cover"
|
||||
SOURCE_CRS = "EPSG:3812"
|
||||
SOURCE_RESOLUTION_M = 1.0
|
||||
SOURCE_VALUE_UNIT = "class_1_11"
|
||||
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."
|
||||
)
|
||||
CLASS_LABELS = {
|
||||
1: "Jaarlijks wisselende kruidlaag",
|
||||
2: "Jaarronde kruidlaag",
|
||||
3: "Naaldbomen hoger dan 3 m",
|
||||
4: "Loofbomen hoger dan 3 m",
|
||||
5: "Naaldbomen tot 3 m",
|
||||
6: "Loofbomen tot 3 m",
|
||||
7: "Kale bodem",
|
||||
8: "Oppervlaktewater",
|
||||
9: "Kunstmatige bodembedekking",
|
||||
10: "Spoorweg",
|
||||
11: "Kunstmatige constructies boven maaiveld",
|
||||
}
|
||||
CLASS_COLORS = {
|
||||
1: (236, 202, 73),
|
||||
2: (161, 201, 78),
|
||||
3: (28, 89, 51),
|
||||
4: (52, 132, 72),
|
||||
5: (78, 125, 70),
|
||||
6: (107, 164, 87),
|
||||
7: (194, 165, 119),
|
||||
8: (44, 129, 185),
|
||||
9: (155, 155, 155),
|
||||
10: (68, 68, 68),
|
||||
11: (183, 72, 67),
|
||||
}
|
||||
|
||||
@staticmethod
|
||||
def _products() -> dict[str, WalousProduct]:
|
||||
products = (
|
||||
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",
|
||||
accuracy_label="Officiele globale nauwkeurigheid 83,30%",
|
||||
),
|
||||
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",
|
||||
accuracy_label="Officiele globale nauwkeurigheid 87,10%",
|
||||
),
|
||||
)
|
||||
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=WalousLandCoverService.ATTRIBUTION,
|
||||
license_note=WalousLandCoverService.LICENSE_NOTE,
|
||||
legend_min_label="WALOUS klasse 1",
|
||||
legend_max_label="WALOUS klasse 11",
|
||||
included_source_values=list(WalousLandCoverService.CLASS_LABELS),
|
||||
limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
|
||||
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,
|
||||
) -> 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)
|
||||
raw = np.asarray(band.filled(WalousLandCoverService.NODATA), 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
|
||||
valid = raw[raw != WalousLandCoverService.NODATA]
|
||||
if valid.size == 0:
|
||||
raise AppError(code="WALOUS_NO_VALID_DATA", message="WALOUS contains no valid cells in this selection", status_code=422)
|
||||
classes = set(np.unique(valid).astype(int).tolist())
|
||||
unexpected = sorted(classes - set(WalousLandCoverService.CLASS_LABELS))
|
||||
if unexpected:
|
||||
raise AppError(code="WALOUS_SOURCE_INVALID_VALUES", message="WALOUS contains classes outside the governed 1-11 legend", details={"unexpected_classes": unexpected}, status_code=409)
|
||||
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),
|
||||
"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 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=WalousLandCoverService.ATTRIBUTION,
|
||||
limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
|
||||
).model_dump(mode="json")
|
||||
|
||||
content, validation = WalousLandCoverService._read_source_window(source_path, scope, resolved)
|
||||
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 = datetime(product.observation_year, 12, 31, 23, 59, 59, tzinfo=UTC)
|
||||
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=datetime(product.observation_year, 1, 1, tzinfo=UTC),
|
||||
valid_to=observed_at,
|
||||
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,
|
||||
"valid_pixel_count": validation["valid_pixel_count"],
|
||||
"classes_present": validation["classes_present"],
|
||||
"bbox_epsg4326": bbox_4326,
|
||||
"bbox_epsg3812": validation["bbox_epsg3812"],
|
||||
"coverage_zones": ["wallonia"],
|
||||
"catalog_url": product.catalog_url,
|
||||
"download_url": product.download_url,
|
||||
"attribution": WalousLandCoverService.ATTRIBUTION,
|
||||
"license_note": WalousLandCoverService.LICENSE_NOTE,
|
||||
"limitation_message": f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
|
||||
},
|
||||
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=WalousLandCoverService.ATTRIBUTION,
|
||||
limitation_message=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
|
||||
).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(band.filled(WalousLandCoverService.NODATA), 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", {3, 4, 5, 6}),
|
||||
("surface_water_area_ha", "Oppervlaktewater", {8}),
|
||||
("artificial_cover_area_ha", "Kunstmatige bedekking en constructies", {9, 10, 11}),
|
||||
("annual_herbaceous_cover_area_ha", "Jaarlijks wisselende kruidlaag", {1}),
|
||||
("permanent_herbaceous_cover_area_ha", "Jaarronde kruidlaag", {2}),
|
||||
("bare_soil_area_ha", "Kale bodem", {7}),
|
||||
]
|
||||
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=f"{WalousLandCoverService.LIMITATION} {product.accuracy_label}.",
|
||||
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(values.filled(WalousLandCoverService.NODATA), 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()
|
||||
Reference in New Issue
Block a user