Files
geointel/backend/app/services/orthophoto_acquisition_service.py
T
Codex fb38eb3e91
GeoIntel CI / docs-smoke (push) Canceled after 0s
GeoIntel CI / contract-smoke (push) Canceled after 0s
feat: add governed hydrology and historical imagery
2026-07-15 12:17:45 +02:00

566 lines
27 KiB
Python

from __future__ import annotations
import hashlib
import io
import json
import math
import warnings
from dataclasses import dataclass
from datetime import UTC, datetime, timedelta
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
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.orthophoto import OrthophotoAcquireRequest, OrthophotoAcquisitionResult, OrthophotoProductRead
from app.services.dataset_service import DatasetService
@dataclass(frozen=True)
class OrthophotoProduct:
key: str
display_name: str
observation_label: str
temporal_granularity: str
native_resolution_m: float
wms_url: str
layer: str
catalog_url: str
limitation_message: str
supports_detection: bool = False
color_mode: str = "rgb"
observed_at: datetime | None = None
valid_from: datetime | None = None
valid_to: datetime | None = None
class OrthophotoAcquisitionService:
PROVIDER = "digitaal_vlaanderen_orthophoto"
ATTRIBUTION = "Bron: Orthofotomozaiek Vlaanderen, Digitaal Vlaanderen"
CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/orthofotomozaiek-middenschalig-winteropnamen-kleur-meest-recent-vlaanderen"
LIMITATION = "Meest recente samengestelde winterorthofoto op het moment van de aanvraag; geen historische opnamedatum per pixel."
HISTORICAL_WINTER_WMS_URL = "https://geo.api.vlaanderen.be/OMW/wms"
HISTORICAL_WINTER_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/wmts-orthofotomozaiek-middenschalig-winteropnamen"
HISTORICAL_SUMMER_WMS_URL = "https://geo.api.vlaanderen.be/OKZ/wms"
HISTORICAL_SUMMER_CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/orthofotomozaiek-kleinschalig-zomeropnamen"
@staticmethod
def _products(settings: Settings) -> dict[str, OrthophotoProduct]:
products: list[OrthophotoProduct] = [
OrthophotoProduct(
key="most_recent",
display_name="Meest recente winterluchtbeeld",
observation_label="Meest recent beschikbaar",
temporal_granularity="snapshot",
native_resolution_m=0.15,
wms_url=settings.orthophoto_wms_url,
layer=settings.orthophoto_wms_layer,
catalog_url=OrthophotoAcquisitionService.CATALOG_URL,
limitation_message=OrthophotoAcquisitionService.LIMITATION,
supports_detection=True,
)
]
for year in range(2025, 2011, -1):
products.append(
OrthophotoProduct(
key=str(year),
display_name=f"Winterluchtbeeld {year}",
observation_label=str(year),
temporal_granularity="year",
native_resolution_m=0.15 if year >= 2022 else 0.25,
wms_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_WMS_URL,
layer=f"OMWRGB{year % 100:02d}VL",
catalog_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_CATALOG_URL,
limitation_message=(
"Officiële samengestelde winterorthofoto voor deze jaargang; de exacte opnamedatum kan per tegel verschillen. "
"Historische beelden worden niet met de actuele GRB-toestand gevalideerd."
),
observed_at=datetime(year, 1, 1, tzinfo=UTC),
valid_from=datetime(year, 1, 1, tzinfo=UTC),
valid_to=datetime(year, 12, 31, 23, 59, 59, tzinfo=UTC),
)
)
for key, start_year, end_year, layer in (
("2008_2011", 2008, 2011, "OMWRGB08_11VL"),
("2005_2007", 2005, 2007, "OMWRGB05_07VL"),
("2000_2003", 2000, 2003, "OMWRGB00_03VL"),
):
products.append(
OrthophotoProduct(
key=key,
display_name=f"Winterluchtbeeld {start_year}-{end_year}",
observation_label=f"{start_year}-{end_year}",
temporal_granularity="period",
native_resolution_m=0.25,
wms_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_WMS_URL,
layer=layer,
catalog_url=OrthophotoAcquisitionService.HISTORICAL_WINTER_CATALOG_URL,
limitation_message=(
"Officiële samengestelde winterorthofoto uit een meerjarige opnameperiode; dit is geen exacte jaaropname. "
"Historische beelden worden niet met de actuele GRB-toestand gevalideerd."
),
observed_at=datetime(start_year, 1, 1, tzinfo=UTC),
valid_from=datetime(start_year, 1, 1, tzinfo=UTC),
valid_to=datetime(end_year, 12, 31, 23, 59, 59, tzinfo=UTC),
)
)
products.extend(
[
OrthophotoProduct(
key="1979_1990",
display_name="Zomerluchtbeeld 1979-1990",
observation_label="1979-1990",
temporal_granularity="period",
native_resolution_m=1.0,
wms_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_WMS_URL,
layer="OKZRGB79_90VL",
catalog_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_CATALOG_URL,
limitation_message="Kleinschalig RGB-mozaïek uit meerdere zomervluchten tussen 1979 en 1990; geen exacte jaartoestand.",
observed_at=datetime(1979, 1, 1, tzinfo=UTC),
valid_from=datetime(1979, 1, 1, tzinfo=UTC),
valid_to=datetime(1990, 12, 31, 23, 59, 59, tzinfo=UTC),
),
OrthophotoProduct(
key="1971",
display_name="Zomerluchtbeeld 1971",
observation_label="1971",
temporal_granularity="year",
native_resolution_m=1.0,
wms_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_WMS_URL,
layer="OKZPAN71VL",
catalog_url=OrthophotoAcquisitionService.HISTORICAL_SUMMER_CATALOG_URL,
limitation_message="Kleinschalig panchromatisch mozaïek uit 1971; zwart-wit en niet geschikt voor het huidige RGB-detectiemodel.",
color_mode="panchromatic",
observed_at=datetime(1971, 1, 1, tzinfo=UTC),
valid_from=datetime(1971, 1, 1, tzinfo=UTC),
valid_to=datetime(1971, 12, 31, 23, 59, 59, tzinfo=UTC),
),
]
)
return {product.key: product for product in products}
@staticmethod
def list_products(settings: Settings | None = None) -> list[dict[str, Any]]:
resolved_settings = settings or get_settings()
return [
OrthophotoProductRead(
key=product.key,
display_name=product.display_name,
observation_label=product.observation_label,
temporal_granularity=product.temporal_granularity,
native_resolution_m=product.native_resolution_m,
supports_detection=product.supports_detection,
color_mode=product.color_mode,
catalog_url=product.catalog_url,
limitation_message=product.limitation_message,
).model_dump()
for product in OrthophotoAcquisitionService._products(resolved_settings).values()
]
@staticmethod
def _product(product_key: str, settings: Settings) -> OrthophotoProduct:
product = OrthophotoAcquisitionService._products(settings).get(product_key.strip().lower())
if product is None:
raise AppError(
code="ORTHOPHOTO_PRODUCT_NOT_SUPPORTED",
message="Select an orthophoto product from the official product registry",
details={"product_key": product_key},
status_code=422,
)
return product
@staticmethod
def _prepared_request(
payload: OrthophotoAcquireRequest,
settings: Settings,
) -> dict[str, Any]:
product = OrthophotoAcquisitionService._product(payload.product_key, settings)
if payload.bbox.crs.upper() != "EPSG:4326":
raise AppError(code="INVALID_CRS", message="Orthophoto selection bbox must use EPSG:4326", status_code=400)
min_x = float(payload.bbox.min_x)
min_y = float(payload.bbox.min_y)
max_x = float(payload.bbox.max_x)
max_y = float(payload.bbox.max_y)
if not all(math.isfinite(value) for value in (min_x, min_y, max_x, max_y)) or min_x >= max_x or min_y >= max_y:
raise AppError(code="INVALID_BBOX", message="Orthophoto selection must be a finite non-empty rectangle", status_code=400)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True)
lambert_bounds = transformer.transform_bounds(min_x, min_y, max_x, max_y, densify_pts=21)
width_m = lambert_bounds[2] - lambert_bounds[0]
height_m = lambert_bounds[3] - lambert_bounds[1]
if width_m < settings.orthophoto_min_side_m or height_m < settings.orthophoto_min_side_m:
raise AppError(
code="ORTHOPHOTO_SELECTION_TOO_SMALL",
message=f"Select an area of at least {settings.orthophoto_min_side_m:.0f} by {settings.orthophoto_min_side_m:.0f} metres",
status_code=422,
)
if width_m > settings.orthophoto_max_side_m or height_m > settings.orthophoto_max_side_m:
raise AppError(
code="ORTHOPHOTO_SELECTION_TOO_LARGE",
message=f"Select an area no larger than {settings.orthophoto_max_side_m:.0f} by {settings.orthophoto_max_side_m:.0f} metres",
details={"width_m": width_m, "height_m": height_m},
status_code=422,
)
width = max(1, math.ceil(width_m / settings.orthophoto_resolution_m))
height = max(1, math.ceil(height_m / settings.orthophoto_resolution_m))
bbox_4326 = [min_x, min_y, max_x, max_y]
bbox_31370 = [float(value) for value in lambert_bounds]
request_identity = {
"provider": OrthophotoAcquisitionService.PROVIDER,
"product_key": product.key,
"wms_url": product.wms_url,
"layer": product.layer,
"bbox_epsg4326": [round(value, 8) for value in bbox_4326],
"bbox_epsg31370": [round(value, 3) for value in bbox_31370],
"width": width,
"height": height,
"resolution_m": settings.orthophoto_resolution_m,
}
request_hash = hashlib.sha256(json.dumps(request_identity, sort_keys=True).encode("utf-8")).hexdigest()
spatial_identity = {
"bbox_epsg4326": request_identity["bbox_epsg4326"],
"width": width,
"height": height,
"resolution_m": settings.orthophoto_resolution_m,
}
spatial_hash = hashlib.sha256(json.dumps(spatial_identity, sort_keys=True).encode("utf-8")).hexdigest()
params = {
"SERVICE": "WMS",
"VERSION": "1.3.0",
"REQUEST": "GetMap",
"LAYERS": product.layer,
"STYLES": "",
"FORMAT": "image/tiff",
"CRS": "EPSG:31370",
"BBOX": ",".join(f"{value:.3f}" for value in bbox_31370),
"WIDTH": str(width),
"HEIGHT": str(height),
}
return {
**request_identity,
"product": product,
"spatial_hash": spatial_hash,
"request_hash": request_hash,
"request_url": f"{product.wms_url}?{urlencode(params)}",
"params": params,
"bbox_epsg4326": bbox_4326,
"bbox_epsg31370": bbox_31370,
}
@staticmethod
def _validate_area_scope(db, project_id: UUID, area_id: UUID | None, bbox_epsg4326: list[float]) -> None:
if not db.get(Project, project_id):
raise AppError(code="PROJECT_NOT_FOUND", message="Project not found", status_code=404)
if area_id is None:
return
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)
selection = box(*bbox_epsg4326)
area_geometry = to_shape(area.geometry)
transformer = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True)
selection_metric = shapely_transform(transformer.transform, selection)
area_metric = shapely_transform(transformer.transform, area_geometry)
overlap_ratio = area_metric.intersection(selection_metric).area / selection_metric.area
if overlap_ratio < 0.99:
raise AppError(
code="ORTHOPHOTO_SELECTION_OUTSIDE_AREA",
message="Keep the orthophoto rectangle inside the selected work area",
details={"coverage_ratio": overlap_ratio},
status_code=422,
)
@staticmethod
def _cached_dataset(
db,
project_id: UUID,
filename: str,
settings: Settings,
product: OrthophotoProduct,
) -> Dataset | None:
if product.key == "most_recent" and settings.orthophoto_cache_ttl_hours <= 0:
return None
candidate = (
db.query(Dataset)
.filter(
Dataset.project_id == project_id,
Dataset.name == filename,
Dataset.source_name == OrthophotoAcquisitionService.PROVIDER,
Dataset.status == "ready",
)
.order_by(Dataset.imported_at.desc())
.first()
)
if not candidate or not candidate.storage_path or not Path(candidate.storage_path).is_file():
return None
imported_at = candidate.imported_at
if imported_at is None:
return None
if imported_at.tzinfo is None:
imported_at = imported_at.replace(tzinfo=UTC)
if product.key == "most_recent" and datetime.now(UTC) - imported_at > timedelta(hours=settings.orthophoto_cache_ttl_hours):
return None
return candidate
@staticmethod
def _fetch(request_url: str, settings: Settings, opener: Callable[..., Any] | None = None) -> tuple[bytes, str]:
request = Request(request_url, headers={"User-Agent": "GeoIntel/0.1 bounded-orthophoto-acquisition"})
open_request = opener or urlopen
try:
with open_request(request, timeout=settings.orthophoto_timeout_seconds) as response:
content_type = str(response.headers.get("Content-Type", ""))
content_length = response.headers.get("Content-Length")
max_bytes = settings.orthophoto_max_response_mb * 1024 * 1024
if content_length and int(content_length) > max_bytes:
raise AppError(code="ORTHOPHOTO_RESPONSE_TOO_LARGE", message="Official orthophoto response exceeds the configured size limit", status_code=502)
content = response.read(max_bytes + 1)
except AppError:
raise
except (HTTPError, URLError, TimeoutError, OSError) as exc:
raise AppError(
code="ORTHOPHOTO_PROVIDER_UNAVAILABLE",
message="The official orthophoto service could not complete the bounded request",
details={"reason": str(exc)},
status_code=502,
) from exc
if len(content) > settings.orthophoto_max_response_mb * 1024 * 1024:
raise AppError(code="ORTHOPHOTO_RESPONSE_TOO_LARGE", message="Official orthophoto response exceeds the configured size limit", status_code=502)
if "image" not in content_type.lower() and "tiff" not in content_type.lower():
preview = content[:300].decode("utf-8", errors="replace")
raise AppError(
code="ORTHOPHOTO_PROVIDER_INVALID_RESPONSE",
message="The official orthophoto service did not return an image",
details={"content_type": content_type, "response_preview": preview},
status_code=502,
)
return content, content_type
@staticmethod
def _georeference_tiff(content: bytes, prepared: dict[str, Any]) -> bytes:
try:
from rasterio.io import MemoryFile
from rasterio.errors import NotGeoreferencedWarning
from rasterio.transform import from_bounds
except ImportError as exc:
raise AppError(code="RASTER_PROCESSING_UNAVAILABLE", message="Rasterio is required for orthophoto acquisition", status_code=503) from exc
try:
with MemoryFile(content) as source_memory:
with warnings.catch_warnings():
warnings.simplefilter("ignore", NotGeoreferencedWarning)
with source_memory.open() as source:
product: OrthophotoProduct = prepared["product"]
minimum_band_count = 1 if product.color_mode == "panchromatic" else 3
if source.width != prepared["width"] or source.height != prepared["height"] or source.count < minimum_band_count:
raise AppError(
code="ORTHOPHOTO_PROVIDER_INVALID_RESPONSE",
message="Official orthophoto dimensions or RGB bands do not match the bounded request",
details={"width": source.width, "height": source.height, "bands": source.count},
status_code=502,
)
image = source.read()
profile = source.profile.copy()
profile.update(
driver="GTiff",
crs="EPSG:31370",
transform=from_bounds(*prepared["bbox_epsg31370"], source.width, source.height),
compress="deflate",
tiled=False,
)
with MemoryFile() as output_memory:
with output_memory.open(**profile) as output:
output.write(image)
output.update_tags(
source=f"Digitaal Vlaanderen WMS {product.layer}",
source_url=prepared["request_url"],
attribution=OrthophotoAcquisitionService.ATTRIBUTION,
acquisition="explicit_bounded_map_selection",
)
return output_memory.read()
except AppError:
raise
except Exception as exc:
raise AppError(
code="ORTHOPHOTO_PROVIDER_INVALID_RESPONSE",
message="The official orthophoto response is not a readable GeoTIFF",
details={"reason": str(exc)},
status_code=502,
) from exc
@staticmethod
def acquire(
db,
project_id: UUID,
payload: OrthophotoAcquireRequest,
*,
settings: Settings | None = None,
opener: Callable[..., Any] | None = None,
) -> dict[str, Any]:
resolved_settings = settings or get_settings()
if not resolved_settings.orthophoto_enabled:
raise AppError(code="ORTHOPHOTO_NOT_CONFIGURED", message="Official orthophoto acquisition is disabled", status_code=503)
prepared = OrthophotoAcquisitionService._prepared_request(payload, resolved_settings)
product: OrthophotoProduct = prepared["product"]
OrthophotoAcquisitionService._validate_area_scope(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
filename = f"orthofoto_{product.key}_{prepared['request_hash'][:12]}.tif"
cached = None if payload.force_refresh else OrthophotoAcquisitionService._cached_dataset(
db,
project_id,
filename,
resolved_settings,
product,
)
if cached is not None:
return OrthophotoAcquisitionResult(
output_dataset_id=cached.id,
reused=True,
provider=OrthophotoAcquisitionService.PROVIDER,
product_key=product.key,
display_name=product.display_name,
observation_label=product.observation_label,
temporal_granularity=product.temporal_granularity,
supports_detection=product.supports_detection,
layer=product.layer,
width=prepared["width"],
height=prepared["height"],
resolution_m=resolved_settings.orthophoto_resolution_m,
bbox_epsg4326=prepared["bbox_epsg4326"],
bbox_epsg31370=prepared["bbox_epsg31370"],
attribution=OrthophotoAcquisitionService.ATTRIBUTION,
limitation_message=product.limitation_message,
).model_dump(mode="json")
raw_content, response_content_type = OrthophotoAcquisitionService._fetch(prepared["request_url"], resolved_settings, opener)
geotiff_content = OrthophotoAcquisitionService._georeference_tiff(raw_content, prepared)
acquired_at = datetime.now(UTC)
observed_at = product.observed_at or acquired_at
dataset = DatasetService.import_raster_bytes(
db,
project_id=project_id,
area_id=payload.area_id,
filename=filename,
content=geotiff_content,
source=f"Digitaal Vlaanderen WMS {product.layer}",
source_name=OrthophotoAcquisitionService.PROVIDER,
temporal_series_key=f"digitaal-vlaanderen:orthophoto:{prepared['spatial_hash'][:24]}",
observed_at=observed_at,
valid_from=product.valid_from or observed_at,
valid_to=product.valid_to,
temporal_granularity=product.temporal_granularity,
source_version=(
f"most_recent_at_{acquired_at.date().isoformat()}"
if product.key == "most_recent"
else product.key
),
content_type="image/tiff",
source_metadata={
"provider": OrthophotoAcquisitionService.PROVIDER,
"service": "WMS",
"service_version": "1.3.0",
"product_key": product.key,
"product_display_name": product.display_name,
"observation_label": product.observation_label,
"observation_date_precision": product.temporal_granularity,
"native_resolution_m": product.native_resolution_m,
"requested_resolution_m": resolved_settings.orthophoto_resolution_m,
"color_mode": product.color_mode,
"supports_detection": product.supports_detection,
"layer": product.layer,
"catalog_url": product.catalog_url,
"attribution": OrthophotoAcquisitionService.ATTRIBUTION,
"license_note": "Gebruik volgens het gebruiksrecht geografische webdiensten van Digitaal Vlaanderen.",
},
provenance_metadata={
"acquisition": "explicit_bounded_map_selection",
"acquired_at": acquired_at.isoformat(),
"request_hash": prepared["request_hash"],
"spatial_hash": prepared["spatial_hash"],
"request_url": prepared["request_url"],
"response_content_type": response_content_type,
"bbox_epsg4326": prepared["bbox_epsg4326"],
"bbox_epsg31370": prepared["bbox_epsg31370"],
"width": prepared["width"],
"height": prepared["height"],
"resolution_m": resolved_settings.orthophoto_resolution_m,
"limitation_message": product.limitation_message,
},
)
return OrthophotoAcquisitionResult(
output_dataset_id=dataset.id,
reused=False,
provider=OrthophotoAcquisitionService.PROVIDER,
product_key=product.key,
display_name=product.display_name,
observation_label=product.observation_label,
temporal_granularity=product.temporal_granularity,
supports_detection=product.supports_detection,
layer=product.layer,
width=prepared["width"],
height=prepared["height"],
resolution_m=resolved_settings.orthophoto_resolution_m,
bbox_epsg4326=prepared["bbox_epsg4326"],
bbox_epsg31370=prepared["bbox_epsg31370"],
attribution=OrthophotoAcquisitionService.ATTRIBUTION,
limitation_message=product.limitation_message,
).model_dump(mode="json")
@staticmethod
def render_png(db, project_id: UUID, dataset_id: UUID, *, max_dimension: int = 1600) -> bytes:
dataset = db.get(Dataset, dataset_id)
if (
dataset is None
or dataset.project_id != project_id
or dataset.source_name != OrthophotoAcquisitionService.PROVIDER
or dataset.status != "ready"
or not dataset.storage_path
or not Path(dataset.storage_path).is_file()
):
raise AppError(code="ORTHOPHOTO_NOT_FOUND", message="Orthophoto dataset not found", status_code=404)
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="Raster preview dependencies are unavailable", status_code=503) from exc
try:
with rasterio.open(dataset.storage_path) as source:
scale = min(1.0, max_dimension / max(source.width, source.height))
width = max(1, round(source.width * scale))
height = max(1, round(source.height * scale))
indexes = [1] if source.count == 1 else list(range(1, min(source.count, 3) + 1))
pixels = source.read(indexes, out_shape=(len(indexes), height, width), resampling=Resampling.bilinear)
if pixels.dtype != np.uint8:
pixels = np.clip(pixels, 0, 255).astype(np.uint8)
if len(indexes) == 1:
image = Image.fromarray(pixels[0])
else:
image = Image.fromarray(np.moveaxis(pixels[:3], 0, 2))
output = io.BytesIO()
image.save(output, format="PNG", optimize=True)
return output.getvalue()
except AppError:
raise
except Exception as exc:
raise AppError(
code="ORTHOPHOTO_PREVIEW_FAILED",
message="The persisted orthophoto could not be rendered",
details={"reason": str(exc)},
status_code=500,
) from exc