Files
geointel/backend/app/services/orthophoto_acquisition_service.py
T
Codex daccd3869a
GeoIntel CI / docs-smoke (push) Canceled after 0s
GeoIntel CI / contract-smoke (push) Canceled after 0s
feat: add map-driven orthophoto analysis
2026-07-15 02:01:03 +02:00

316 lines
16 KiB
Python

from __future__ import annotations
import hashlib
import json
import math
import warnings
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
from app.services.dataset_service import DatasetService
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."
@staticmethod
def _prepared_request(
payload: OrthophotoAcquireRequest,
settings: Settings,
) -> dict[str, Any]:
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,
"wms_url": settings.orthophoto_wms_url,
"layer": settings.orthophoto_wms_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()
params = {
"SERVICE": "WMS",
"VERSION": "1.3.0",
"REQUEST": "GetMap",
"LAYERS": settings.orthophoto_wms_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,
"request_hash": request_hash,
"request_url": f"{settings.orthophoto_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) -> Dataset | None:
if 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 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:
if source.width != prepared["width"] or source.height != prepared["height"] or source.count < 3:
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="Digitaal Vlaanderen OMWRGBMRVL WMS Ortho 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)
OrthophotoAcquisitionService._validate_area_scope(db, project_id, payload.area_id, prepared["bbox_epsg4326"])
filename = f"orthofoto_selectie_{prepared['request_hash'][:12]}.tif"
cached = None if payload.force_refresh else OrthophotoAcquisitionService._cached_dataset(db, project_id, filename, resolved_settings)
if cached is not None:
return OrthophotoAcquisitionResult(
output_dataset_id=cached.id,
reused=True,
provider=OrthophotoAcquisitionService.PROVIDER,
layer=resolved_settings.orthophoto_wms_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=OrthophotoAcquisitionService.LIMITATION,
).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)
dataset = DatasetService.import_raster_bytes(
db,
project_id=project_id,
area_id=payload.area_id,
filename=filename,
content=geotiff_content,
source="Digitaal Vlaanderen OMWRGBMRVL WMS",
source_name=OrthophotoAcquisitionService.PROVIDER,
source_version=f"most_recent_at_{acquired_at.date().isoformat()}",
content_type="image/tiff",
source_metadata={
"provider": OrthophotoAcquisitionService.PROVIDER,
"service": "WMS",
"service_version": "1.3.0",
"layer": resolved_settings.orthophoto_wms_layer,
"catalog_url": OrthophotoAcquisitionService.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"],
"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": OrthophotoAcquisitionService.LIMITATION,
},
)
return OrthophotoAcquisitionResult(
output_dataset_id=dataset.id,
reused=False,
provider=OrthophotoAcquisitionService.PROVIDER,
layer=resolved_settings.orthophoto_wms_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=OrthophotoAcquisitionService.LIMITATION,
).model_dump(mode="json")