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")