730 lines
34 KiB
Python
730 lines
34 KiB
Python
from __future__ import annotations
|
|
|
|
import hashlib
|
|
import json
|
|
import math
|
|
import time
|
|
from dataclasses import dataclass
|
|
from datetime import UTC, datetime
|
|
from email.parser import BytesParser
|
|
from email.policy import default
|
|
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
|
|
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.dhmv import DhmvAcquireRequest, DhmvAcquisitionResult, DhmvProductRead
|
|
from app.services.dataset_service import DatasetService
|
|
|
|
|
|
@dataclass(frozen=True)
|
|
class DhmvProduct:
|
|
key: str
|
|
display_name: str
|
|
surface_model: str
|
|
coverage_id: str
|
|
native_resolution_m: float
|
|
catalog_url: str
|
|
limitation_message: str
|
|
|
|
|
|
class DhmvAcquisitionService:
|
|
PROVIDER = "digitaal_vlaanderen_dhmv"
|
|
SOURCE_CRS = "EPSG:31370"
|
|
VERTICAL_REFERENCE = "TAW (Tweede Algemene Waterpassing)"
|
|
ACQUISITION_PERIOD = "2013-2015"
|
|
SOURCE_VERSION = "DHMV II 2014.01"
|
|
NODATA = -9999.0
|
|
ATTRIBUTION = "Bron: Digitaal Vlaanderen, Digitaal Hoogtemodel Vlaanderen II"
|
|
LICENSE_NOTE = "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen."
|
|
WCS_TILE_SIDE_M = 10_000.0
|
|
WCS_REQUEST_INTERVAL_SECONDS = 2.0
|
|
WCS_RETRY_DELAY_SECONDS = 4.0
|
|
WCS_TRANSIENT_STATUS_CODES = frozenset({400, 429, 502, 503, 504})
|
|
WCS_EDGE_RESOLUTION_REL_TOLERANCE = 0.05
|
|
WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M = 0.25
|
|
DTM_CATALOG_URL = (
|
|
"https://www.vlaanderen.be/datavindplaats/catalogus/"
|
|
"digitaal-hoogtemodel-vlaanderen-ii-dtm-raster-1-m"
|
|
)
|
|
DSM_CATALOG_URL = (
|
|
"https://www.vlaanderen.be/datavindplaats/catalogus/"
|
|
"digitaal-hoogtemodel-vlaanderen-ii-dsm-raster-1-m"
|
|
)
|
|
|
|
@staticmethod
|
|
def _products() -> dict[str, DhmvProduct]:
|
|
products = (
|
|
DhmvProduct(
|
|
key="dtm_1m",
|
|
display_name="DHMV II terreinmodel (DTM)",
|
|
surface_model="terrain",
|
|
coverage_id="DHMVII_DTM_1m",
|
|
native_resolution_m=1.0,
|
|
catalog_url=DhmvAcquisitionService.DTM_CATALOG_URL,
|
|
limitation_message=(
|
|
"Maaiveldhoogte uit de opnameperiode 2013-2015. Gebouwen en andere objecten zijn verwijderd. "
|
|
"Afstroming is een afgeleide interpretatie; dit product bevat geen waterdiepte."
|
|
),
|
|
),
|
|
DhmvProduct(
|
|
key="dsm_1m",
|
|
display_name="DHMV II oppervlaktemodel (DSM)",
|
|
surface_model="surface",
|
|
coverage_id="DHMVII_DSM_1m",
|
|
native_resolution_m=1.0,
|
|
catalog_url=DhmvAcquisitionService.DSM_CATALOG_URL,
|
|
limitation_message=(
|
|
"Oppervlaktehoogte uit de opnameperiode 2013-2015, inclusief gebouwen en vegetatie. "
|
|
"Dit is geen maaiveldmodel, waterdiepte of rechtstreeks gebouwhoogteproduct."
|
|
),
|
|
),
|
|
)
|
|
return {product.key: product for product in products}
|
|
|
|
@staticmethod
|
|
def list_products() -> list[dict[str, Any]]:
|
|
return [
|
|
DhmvProductRead(
|
|
key=product.key,
|
|
display_name=product.display_name,
|
|
surface_model=product.surface_model,
|
|
coverage_id=product.coverage_id,
|
|
native_resolution_m=product.native_resolution_m,
|
|
source_crs=DhmvAcquisitionService.SOURCE_CRS,
|
|
vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
|
|
acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
|
|
catalog_url=product.catalog_url,
|
|
attribution=DhmvAcquisitionService.ATTRIBUTION,
|
|
limitation_message=product.limitation_message,
|
|
).model_dump()
|
|
for product in DhmvAcquisitionService._products().values()
|
|
]
|
|
|
|
@staticmethod
|
|
def _product(product_key: str) -> DhmvProduct:
|
|
product = DhmvAcquisitionService._products().get(product_key.strip().lower())
|
|
if product is None:
|
|
raise AppError(
|
|
code="DHMV_PRODUCT_NOT_SUPPORTED",
|
|
message="Select DTM or DSM from the governed DHMV II product registry",
|
|
details={"product_key": product_key},
|
|
status_code=422,
|
|
)
|
|
return product
|
|
|
|
@staticmethod
|
|
def _prepared_request(payload: DhmvAcquireRequest, settings: Settings) -> dict[str, Any]:
|
|
if not settings.dhmv_enabled:
|
|
raise AppError(code="DHMV_NOT_CONFIGURED", message="DHMV acquisition is disabled", status_code=503)
|
|
product = DhmvAcquisitionService._product(payload.product_key)
|
|
resolution_m = float(payload.resolution_m or settings.dhmv_resolution_m)
|
|
if resolution_m < product.native_resolution_m or resolution_m > 10.0:
|
|
raise AppError(
|
|
code="DHMV_RESOLUTION_NOT_SUPPORTED",
|
|
message="DHMV analysis resolution must be between the native 1 metre and 10 metres",
|
|
status_code=422,
|
|
)
|
|
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 payload.bbox.min_x >= payload.bbox.max_x or payload.bbox.min_y >= payload.bbox.max_y:
|
|
raise AppError(code="INVALID_BBOX", message="DHMV selection must be a finite non-empty rectangle", status_code=400)
|
|
|
|
transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.SOURCE_CRS, always_xy=True)
|
|
lambert_bounds = transformer.transform_bounds(*values, densify_pts=21)
|
|
width_m = float(lambert_bounds[2] - lambert_bounds[0])
|
|
height_m = float(lambert_bounds[3] - lambert_bounds[1])
|
|
if width_m < settings.dhmv_min_side_m or height_m < settings.dhmv_min_side_m:
|
|
raise AppError(
|
|
code="DHMV_SELECTION_TOO_SMALL",
|
|
message=f"Select an area of at least {settings.dhmv_min_side_m:g} by {settings.dhmv_min_side_m:g} metres",
|
|
status_code=422,
|
|
)
|
|
if width_m > settings.dhmv_max_side_m or height_m > settings.dhmv_max_side_m:
|
|
raise AppError(
|
|
code="DHMV_SELECTION_TOO_LARGE",
|
|
message=f"Select an area no larger than {settings.dhmv_max_side_m:g} by {settings.dhmv_max_side_m:g} metres",
|
|
details={"width_m": width_m, "height_m": height_m},
|
|
status_code=422,
|
|
)
|
|
width = max(1, math.ceil(width_m / resolution_m))
|
|
height = max(1, math.ceil(height_m / resolution_m))
|
|
if width * height > settings.dhmv_max_pixels:
|
|
raise AppError(
|
|
code="DHMV_SELECTION_TOO_LARGE",
|
|
message="DHMV selection exceeds the configured raster cell limit",
|
|
details={"pixel_count": width * height, "max_pixels": settings.dhmv_max_pixels},
|
|
status_code=422,
|
|
)
|
|
|
|
bbox_4326 = [float(value) for value in values]
|
|
bbox_31370 = [float(value) for value in lambert_bounds]
|
|
request_identity = {
|
|
"provider": DhmvAcquisitionService.PROVIDER,
|
|
"coverage_id": product.coverage_id,
|
|
"bbox_epsg4326": [round(value, 8) for value in bbox_4326],
|
|
"bbox_epsg31370": [round(value, 3) for value in bbox_31370],
|
|
"resolution_m": resolution_m,
|
|
"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()
|
|
params = {
|
|
"SERVICE": "WCS",
|
|
"VERSION": "2.0.1",
|
|
"REQUEST": "GetCoverage",
|
|
"COVERAGEID": product.coverage_id,
|
|
"FORMAT": "image/tiff",
|
|
"SUBSET": [
|
|
f"x({bbox_31370[0]:.3f},{bbox_31370[2]:.3f})",
|
|
f"y({bbox_31370[1]:.3f},{bbox_31370[3]:.3f})",
|
|
],
|
|
"SCALEFACTOR": f"{resolution_m / product.native_resolution_m:g}",
|
|
}
|
|
query = [
|
|
("SERVICE", params["SERVICE"]),
|
|
("VERSION", params["VERSION"]),
|
|
("REQUEST", params["REQUEST"]),
|
|
("COVERAGEID", params["COVERAGEID"]),
|
|
("FORMAT", params["FORMAT"]),
|
|
("SUBSET", params["SUBSET"][0]),
|
|
("SUBSET", params["SUBSET"][1]),
|
|
("SCALEFACTOR", params["SCALEFACTOR"]),
|
|
]
|
|
return {
|
|
**request_identity,
|
|
"product": product,
|
|
"request_hash": request_hash,
|
|
"request_url": f"{settings.dhmv_wcs_url}?{urlencode(query)}",
|
|
"params": params,
|
|
"bbox_epsg4326": bbox_4326,
|
|
"bbox_epsg31370": bbox_31370,
|
|
"width": width,
|
|
"height": height,
|
|
}
|
|
|
|
@staticmethod
|
|
def _wcs_request_url(
|
|
settings: Settings,
|
|
product: DhmvProduct,
|
|
bounds: tuple[float, float, float, float],
|
|
resolution_m: float,
|
|
) -> str:
|
|
query = [
|
|
("SERVICE", "WCS"),
|
|
("VERSION", "2.0.1"),
|
|
("REQUEST", "GetCoverage"),
|
|
("COVERAGEID", product.coverage_id),
|
|
("FORMAT", "image/tiff"),
|
|
("SUBSET", f"x({bounds[0]:.3f},{bounds[2]:.3f})"),
|
|
("SUBSET", f"y({bounds[1]:.3f},{bounds[3]:.3f})"),
|
|
("SCALEFACTOR", f"{resolution_m / product.native_resolution_m:g}"),
|
|
]
|
|
return f"{settings.dhmv_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"]
|
|
tiles: list[tuple[float, float, float, float]] = []
|
|
y = min_y
|
|
while y < max_y:
|
|
tile_max_y = min(y + DhmvAcquisitionService.WCS_TILE_SIDE_M, max_y)
|
|
x = min_x
|
|
while x < max_x:
|
|
tile_max_x = min(x + DhmvAcquisitionService.WCS_TILE_SIDE_M, 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="DHMV_SELECTION_OUTSIDE_AREA",
|
|
message="The DHMV selection does not overlap the selected work area",
|
|
status_code=422,
|
|
)
|
|
return intersection
|
|
|
|
@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 == DhmvAcquisitionService.PROVIDER,
|
|
Dataset.status == "ready",
|
|
)
|
|
.order_by(Dataset.imported_at.desc())
|
|
.first()
|
|
)
|
|
if candidate and candidate.storage_path and Path(candidate.storage_path).is_file():
|
|
return candidate
|
|
return None
|
|
|
|
@staticmethod
|
|
def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]:
|
|
request = Request(
|
|
request_url,
|
|
headers={
|
|
"Accept": "*/*",
|
|
"User-Agent": "GeoIntel/0.1 bounded-dhmv-acquisition",
|
|
},
|
|
)
|
|
max_bytes = settings.dhmv_max_response_mb * 1024 * 1024
|
|
try:
|
|
with (opener or guarded_opener(request_url))(request, timeout=settings.dhmv_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="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502)
|
|
content = response.read(max_bytes + 1)
|
|
except AppError:
|
|
raise
|
|
except HTTPError as exc:
|
|
preview = exc.read(300).decode("utf-8", errors="replace")
|
|
raise AppError(
|
|
code="DHMV_PROVIDER_UNAVAILABLE",
|
|
message="The official DHMV WCS could not complete the bounded request",
|
|
details={
|
|
"reason": str(exc),
|
|
"provider_status_code": int(exc.code),
|
|
"response_preview": preview,
|
|
},
|
|
status_code=502,
|
|
) from exc
|
|
except (URLError, TimeoutError, OSError) as exc:
|
|
raise AppError(
|
|
code="DHMV_PROVIDER_UNAVAILABLE",
|
|
message="The official DHMV WCS could not complete the bounded request",
|
|
details={"reason": str(exc)},
|
|
status_code=502,
|
|
) from exc
|
|
if len(content) > max_bytes:
|
|
raise AppError(code="DHMV_RESPONSE_TOO_LARGE", message="Official DHMV response exceeds the configured size limit", status_code=502)
|
|
return content, content_type
|
|
|
|
@staticmethod
|
|
def _extract_geotiff(content: bytes, content_type: str) -> bytes:
|
|
if content.startswith((b"II*\x00", b"MM\x00*")):
|
|
return content
|
|
if "multipart" not in content_type.lower():
|
|
preview = content[:300].decode("utf-8", errors="replace")
|
|
raise AppError(
|
|
code="DHMV_PROVIDER_INVALID_RESPONSE",
|
|
message="The official DHMV service did not return a GeoTIFF coverage",
|
|
details={"content_type": content_type, "response_preview": preview},
|
|
status_code=502,
|
|
)
|
|
message = BytesParser(policy=default).parsebytes(
|
|
f"Content-Type: {content_type}\r\nMIME-Version: 1.0\r\n\r\n".encode() + content
|
|
)
|
|
for part in message.iter_parts():
|
|
payload = part.get_payload(decode=True) or b""
|
|
if part.get_content_type() == "image/tiff" and payload.startswith((b"II*\x00", b"MM\x00*")):
|
|
return payload
|
|
raise AppError(
|
|
code="DHMV_PROVIDER_INVALID_RESPONSE",
|
|
message="The official DHMV multipart response contains no valid GeoTIFF coverage",
|
|
status_code=502,
|
|
)
|
|
|
|
@staticmethod
|
|
def _mosaic_geotiffs(
|
|
coverages: list[bytes],
|
|
expected_resolution_m: float | None = None,
|
|
diagnostics: dict[str, Any] | None = None,
|
|
) -> bytes:
|
|
if len(coverages) == 1:
|
|
return coverages[0]
|
|
try:
|
|
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 tiled DHMV coverages",
|
|
status_code=503,
|
|
) from exc
|
|
|
|
memories = [MemoryFile(content) for content in coverages]
|
|
sources = []
|
|
try:
|
|
sources = [memory.open() for memory in memories]
|
|
target_resolution = float(expected_resolution_m or abs(float(sources[0].res[0])))
|
|
invalid_crs = [index for index, source in enumerate(sources) if source.crs is None or source.crs.to_epsg() != 31370]
|
|
invalid_bands = [index for index, source in enumerate(sources) if source.count != 1]
|
|
tile_resolutions = [
|
|
[abs(float(source.res[0])), abs(float(source.res[1]))]
|
|
for source in sources
|
|
]
|
|
invalid_resolution = [
|
|
{
|
|
"tile_index": index,
|
|
"resolution": tile_resolutions[index],
|
|
}
|
|
for index, source in enumerate(sources)
|
|
if not all(
|
|
math.isclose(
|
|
abs(float(value)),
|
|
target_resolution,
|
|
rel_tol=DhmvAcquisitionService.WCS_EDGE_RESOLUTION_REL_TOLERANCE,
|
|
abs_tol=DhmvAcquisitionService.WCS_EDGE_RESOLUTION_ABS_TOLERANCE_M,
|
|
)
|
|
for value in source.res
|
|
)
|
|
]
|
|
if invalid_crs or invalid_bands or invalid_resolution:
|
|
raise AppError(
|
|
code="DHMV_TILE_MISMATCH",
|
|
message="DHMV coverage tiles do not match the governed CRS, band layout and resolution",
|
|
details={
|
|
"invalid_crs_tile_indexes": invalid_crs,
|
|
"invalid_band_tile_indexes": invalid_bands,
|
|
"invalid_resolution_tiles": invalid_resolution,
|
|
"expected_resolution_m": target_resolution,
|
|
},
|
|
status_code=502,
|
|
)
|
|
harmonized_tile_indexes = [
|
|
index
|
|
for index, resolution in enumerate(tile_resolutions)
|
|
if not all(
|
|
math.isclose(value, target_resolution, rel_tol=0.02, abs_tol=0.05)
|
|
for value in resolution
|
|
)
|
|
]
|
|
if diagnostics is not None:
|
|
diagnostics.update(
|
|
{
|
|
"source_tile_resolutions_m": tile_resolutions,
|
|
"target_resolution_m": target_resolution,
|
|
"harmonized_tile_indexes": harmonized_tile_indexes,
|
|
"harmonization_method": "rasterio_merge_target_resolution" if harmonized_tile_indexes else None,
|
|
}
|
|
)
|
|
mosaic, transform = merge(
|
|
sources,
|
|
res=(target_resolution, target_resolution),
|
|
nodata=DhmvAcquisitionService.NODATA,
|
|
dtype="float32",
|
|
)
|
|
profile = sources[0].profile.copy()
|
|
profile.pop("blockxsize", None)
|
|
profile.pop("blockysize", None)
|
|
profile.update(
|
|
driver="GTiff",
|
|
width=int(mosaic.shape[2]),
|
|
height=int(mosaic.shape[1]),
|
|
count=1,
|
|
dtype="float32",
|
|
crs=DhmvAcquisitionService.SOURCE_CRS,
|
|
transform=transform,
|
|
nodata=DhmvAcquisitionService.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="DHMV_TILE_MOSAIC_FAILED",
|
|
message="DHMV coverage tiles could not be assembled into one georeferenced raster",
|
|
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: DhmvProduct = prepared["product"]
|
|
request_urls = [
|
|
DhmvAcquisitionService._wcs_request_url(
|
|
settings,
|
|
product,
|
|
bounds,
|
|
prepared["resolution_m"],
|
|
)
|
|
for bounds in DhmvAcquisitionService._tile_bounds(prepared)
|
|
]
|
|
raw_hash = hashlib.sha256()
|
|
coverage_hash = hashlib.sha256()
|
|
content_types: list[str] = []
|
|
coverages: list[bytes] = []
|
|
for index, request_url in enumerate(request_urls):
|
|
if index > 0 and opener is None:
|
|
time.sleep(DhmvAcquisitionService.WCS_REQUEST_INTERVAL_SECONDS)
|
|
try:
|
|
raw_content, content_type = DhmvAcquisitionService._fetch(request_url, settings, opener)
|
|
except AppError as exc:
|
|
provider_status = (exc.details or {}).get("provider_status_code")
|
|
if opener is not None or provider_status not in DhmvAcquisitionService.WCS_TRANSIENT_STATUS_CODES:
|
|
raise
|
|
time.sleep(DhmvAcquisitionService.WCS_RETRY_DELAY_SECONDS)
|
|
raw_content, content_type = DhmvAcquisitionService._fetch(request_url, settings, opener)
|
|
coverage_content = DhmvAcquisitionService._extract_geotiff(raw_content, content_type)
|
|
raw_hash.update(len(raw_content).to_bytes(8, "big"))
|
|
raw_hash.update(raw_content)
|
|
coverage_hash.update(len(coverage_content).to_bytes(8, "big"))
|
|
coverage_hash.update(coverage_content)
|
|
content_types.append(content_type)
|
|
coverages.append(coverage_content)
|
|
mosaic_diagnostics: dict[str, Any] = {}
|
|
mosaic = DhmvAcquisitionService._mosaic_geotiffs(
|
|
coverages,
|
|
prepared["resolution_m"],
|
|
diagnostics=mosaic_diagnostics,
|
|
)
|
|
return mosaic, {
|
|
"tile_count": len(request_urls),
|
|
"request_urls": request_urls,
|
|
"response_content_types": content_types,
|
|
"response_sha256": raw_hash.hexdigest(),
|
|
"coverage_sha256": coverage_hash.hexdigest(),
|
|
"grid_harmonization": mosaic_diagnostics,
|
|
}
|
|
|
|
@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 DHMV validation", status_code=503) from exc
|
|
|
|
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="DHMV_INVALID_CRS", message="DHMV coverage must use EPSG:31370", status_code=502)
|
|
if source.count != 1:
|
|
raise AppError(code="DHMV_INVALID_BANDS", message="DHMV coverage must contain exactly one elevation band", status_code=502)
|
|
resolution = max(abs(float(source.res[0])), abs(float(source.res[1])))
|
|
if not math.isclose(resolution, prepared["resolution_m"], rel_tol=0.02, abs_tol=0.05):
|
|
raise AppError(
|
|
code="DHMV_INVALID_RESOLUTION",
|
|
message="DHMV coverage resolution differs from the governed request",
|
|
details={"expected_m": prepared["resolution_m"], "actual_m": resolution},
|
|
status_code=502,
|
|
)
|
|
transformer = Transformer.from_crs("EPSG:4326", DhmvAcquisitionService.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")
|
|
nodata = float(source.nodata if source.nodata is not None else DhmvAcquisitionService.NODATA)
|
|
invalid = ~np.isfinite(np.asarray(band.filled(np.nan), dtype="float64"))
|
|
combined_mask = np.ma.getmaskarray(band) | invalid | (np.asarray(band) == nodata)
|
|
normalized = np.ma.array(np.asarray(band, dtype="float32"), mask=combined_mask)
|
|
valid_pixel_count = int(normalized.count())
|
|
if valid_pixel_count == 0:
|
|
raise AppError(code="DHMV_NO_VALID_DATA", message="DHMV coverage contains no valid elevation cells in this selection", status_code=422)
|
|
profile = source.profile.copy()
|
|
profile.pop("blockxsize", None)
|
|
profile.pop("blockysize", None)
|
|
profile.update(
|
|
driver="GTiff",
|
|
width=int(normalized.shape[1]),
|
|
height=int(normalized.shape[0]),
|
|
count=1,
|
|
dtype="float32",
|
|
crs=DhmvAcquisitionService.SOURCE_CRS,
|
|
transform=transform,
|
|
nodata=DhmvAcquisitionService.NODATA,
|
|
compress="deflate",
|
|
predictor=3,
|
|
)
|
|
with MemoryFile() as output_memory:
|
|
with output_memory.open(**profile) as output:
|
|
output.write(normalized.filled(DhmvAcquisitionService.NODATA), 1)
|
|
normalized_content = output_memory.read()
|
|
valid_values = normalized.compressed().astype("float64")
|
|
return normalized_content, {
|
|
"width": int(normalized.shape[1]),
|
|
"height": int(normalized.shape[0]),
|
|
"valid_pixel_count": valid_pixel_count,
|
|
"nodata_value": DhmvAcquisitionService.NODATA,
|
|
"resolution_m": resolution,
|
|
"minimum_m_taw": float(valid_values.min()),
|
|
"maximum_m_taw": float(valid_values.max()),
|
|
}
|
|
except AppError:
|
|
raise
|
|
except Exception as exc:
|
|
raise AppError(
|
|
code="DHMV_RASTER_INVALID",
|
|
message="The official DHMV response could not be validated as a georeferenced elevation raster",
|
|
details={"reason": str(exc)},
|
|
status_code=502,
|
|
) from exc
|
|
|
|
@staticmethod
|
|
def acquire(
|
|
db,
|
|
project_id: UUID,
|
|
payload: DhmvAcquireRequest,
|
|
*,
|
|
settings: Settings | None = None,
|
|
opener: Callable[..., Any] | None = None,
|
|
) -> dict[str, Any]:
|
|
resolved_settings = settings or get_settings()
|
|
prepared = DhmvAcquisitionService._prepared_request(payload, resolved_settings)
|
|
product: DhmvProduct = prepared["product"]
|
|
scope_geometry = DhmvAcquisitionService._scope_geometry(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
|
|
resolution_token = f"{prepared['resolution_m']:g}".replace(".", "p")
|
|
filename = f"dhmvii_{product.surface_model}_{resolution_token}m_{prepared['request_hash'][:12]}.tif"
|
|
if not payload.force_refresh:
|
|
cached = DhmvAcquisitionService._cached_dataset(db, project_id, filename)
|
|
if cached is not None:
|
|
source_metadata = cached.source_metadata or {}
|
|
raster_metadata = cached.metadata_json or {}
|
|
return DhmvAcquisitionResult(
|
|
output_dataset_id=cached.id,
|
|
reused=True,
|
|
provider=DhmvAcquisitionService.PROVIDER,
|
|
product_key=product.key,
|
|
display_name=product.display_name,
|
|
surface_model=product.surface_model,
|
|
coverage_id=product.coverage_id,
|
|
native_resolution_m=product.native_resolution_m,
|
|
resolution_m=float(source_metadata.get("analysis_resolution_m", prepared["resolution_m"])),
|
|
width=int(raster_metadata.get("width", prepared["width"])),
|
|
height=int(raster_metadata.get("height", prepared["height"])),
|
|
valid_pixel_count=int(source_metadata.get("valid_pixel_count", 0)),
|
|
nodata_value=float(raster_metadata.get("nodata", DhmvAcquisitionService.NODATA)),
|
|
bbox_epsg4326=prepared["bbox_epsg4326"],
|
|
bbox_epsg31370=prepared["bbox_epsg31370"],
|
|
vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
|
|
acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
|
|
attribution=DhmvAcquisitionService.ATTRIBUTION,
|
|
limitation_message=product.limitation_message,
|
|
).model_dump(mode="json")
|
|
|
|
coverage_content, transfer = DhmvAcquisitionService._fetch_coverage(prepared, resolved_settings, opener)
|
|
normalized_content, validation = DhmvAcquisitionService._normalize_raster(coverage_content, scope_geometry, prepared)
|
|
acquired_at = datetime.now(UTC)
|
|
dataset = DatasetService.import_raster_bytes(
|
|
db,
|
|
project_id=project_id,
|
|
area_id=payload.area_id,
|
|
filename=filename,
|
|
content=normalized_content,
|
|
source=f"Digitaal Vlaanderen WCS {product.coverage_id}",
|
|
source_name=DhmvAcquisitionService.PROVIDER,
|
|
temporal_series_key=f"digitaal-vlaanderen:dhmvii:{product.key}:{prepared['request_hash'][:24]}",
|
|
observed_at=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC),
|
|
valid_from=datetime(2013, 1, 1, tzinfo=UTC),
|
|
valid_to=datetime(2015, 12, 31, 23, 59, 59, tzinfo=UTC),
|
|
temporal_granularity="period",
|
|
source_version=DhmvAcquisitionService.SOURCE_VERSION,
|
|
content_type="image/tiff",
|
|
source_metadata={
|
|
"provider": DhmvAcquisitionService.PROVIDER,
|
|
"service": "WCS",
|
|
"service_version": "2.0.1",
|
|
"product_key": product.key,
|
|
"product_display_name": product.display_name,
|
|
"surface_model": product.surface_model,
|
|
"coverage_id": product.coverage_id,
|
|
"native_resolution_m": product.native_resolution_m,
|
|
"analysis_resolution_m": validation["resolution_m"],
|
|
"source_crs": DhmvAcquisitionService.SOURCE_CRS,
|
|
"vertical_reference": DhmvAcquisitionService.VERTICAL_REFERENCE,
|
|
"vertical_unit": "m",
|
|
"acquisition_period": DhmvAcquisitionService.ACQUISITION_PERIOD,
|
|
"observation_date_precision": "period",
|
|
"nodata_value": validation["nodata_value"],
|
|
"valid_pixel_count": validation["valid_pixel_count"],
|
|
"minimum_m_taw": validation["minimum_m_taw"],
|
|
"maximum_m_taw": validation["maximum_m_taw"],
|
|
"bbox_epsg4326": prepared["bbox_epsg4326"],
|
|
"bbox_epsg31370": prepared["bbox_epsg31370"],
|
|
"catalog_url": product.catalog_url,
|
|
"attribution": DhmvAcquisitionService.ATTRIBUTION,
|
|
"license_note": DhmvAcquisitionService.LICENSE_NOTE,
|
|
"theme": "elevation",
|
|
"coverage_scope": "municipality" if payload.area_id else "bounded_selection",
|
|
},
|
|
provenance_metadata={
|
|
"acquisition": "explicit_bounded_tiled_wcs_coverage",
|
|
"acquired_at": acquired_at.isoformat(),
|
|
"request_hash": prepared["request_hash"],
|
|
"request_url": prepared["request_url"],
|
|
"tile_count": transfer["tile_count"],
|
|
"tile_request_urls": transfer["request_urls"],
|
|
"response_content_types": transfer["response_content_types"],
|
|
"response_sha256": transfer["response_sha256"],
|
|
"coverage_sha256": transfer["coverage_sha256"],
|
|
"grid_harmonization": transfer["grid_harmonization"],
|
|
"normalized_sha256": hashlib.sha256(normalized_content).hexdigest(),
|
|
"bbox_epsg4326": prepared["bbox_epsg4326"],
|
|
"bbox_epsg31370": prepared["bbox_epsg31370"],
|
|
"requested_resolution_m": prepared["resolution_m"],
|
|
"clipped_to_area_id": str(payload.area_id) if payload.area_id else None,
|
|
"validation": validation,
|
|
"limitation_message": product.limitation_message,
|
|
"water_depth_available": False,
|
|
"water_volume_available": False,
|
|
},
|
|
)
|
|
return DhmvAcquisitionResult(
|
|
output_dataset_id=dataset.id,
|
|
reused=False,
|
|
provider=DhmvAcquisitionService.PROVIDER,
|
|
product_key=product.key,
|
|
display_name=product.display_name,
|
|
surface_model=product.surface_model,
|
|
coverage_id=product.coverage_id,
|
|
native_resolution_m=product.native_resolution_m,
|
|
resolution_m=validation["resolution_m"],
|
|
width=validation["width"],
|
|
height=validation["height"],
|
|
valid_pixel_count=validation["valid_pixel_count"],
|
|
nodata_value=validation["nodata_value"],
|
|
bbox_epsg4326=prepared["bbox_epsg4326"],
|
|
bbox_epsg31370=prepared["bbox_epsg31370"],
|
|
vertical_reference=DhmvAcquisitionService.VERTICAL_REFERENCE,
|
|
acquisition_period=DhmvAcquisitionService.ACQUISITION_PERIOD,
|
|
attribution=DhmvAcquisitionService.ATTRIBUTION,
|
|
limitation_message=product.limitation_message,
|
|
).model_dump(mode="json")
|