from __future__ import annotations from dataclasses import dataclass from datetime import UTC, datetime, timedelta import hashlib import json import math from pathlib import Path from typing import Any, Callable from urllib.error import HTTPError, URLError from urllib.parse import parse_qsl, urlencode, urljoin, urlparse, urlunparse from urllib.request import Request from uuid import UUID from geoalchemy2.shape import to_shape from pyproj import Transformer from shapely.geometry import LineString, MultiLineString, MultiPolygon, Polygon, box, mapping, shape from shapely.ops import unary_union from shapely.validation import make_valid 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.grb import GrbAcquireRequest, GrbAcquisitionResult, GrbProductRead from app.services.dataset_service import DatasetService @dataclass(frozen=True) class GrbCollection: name: str geometry_dimension: int @dataclass(frozen=True) class GrbProduct: key: str display_name: str reference_layer_name: str layer_type: str collections: tuple[GrbCollection, ...] geometry_types: tuple[str, ...] metric_key: str metric_method: str metric_label: str metric_unit: str metric_dimension: int metric_warning: str limitation_message: str class GrbAcquisitionService: PROVIDER = "grb" SOURCE_CRS = "EPSG:4326" OGC_CRS84_URI = "http://www.opengis.net/def/crs/OGC/1.3/CRS84" AUTHORITY_LEVEL = "authoritative" ATTRIBUTION = "Bron: Grootschalig Referentie Bestand Vlaanderen, Digitaal Vlaanderen" LICENSE_NOTE = "Hergebruik volgens de open-datavoorwaarden en bronvermelding van Digitaal Vlaanderen." CATALOG_URL = "https://www.vlaanderen.be/datavindplaats/catalogus/basiskaart-vlaanderen-grb" @staticmethod def _products() -> dict[str, GrbProduct]: products = ( GrbProduct( key="buildings", display_name="GRB gebouwcontouren", reference_layer_name="buildings", layer_type="building", collections=(GrbCollection("GBG", 2),), geometry_types=("Polygon", "MultiPolygon"), metric_key="footprint_area", metric_method="intersection_area", metric_label="Bebouwde grondoppervlakte", metric_unit="ha", metric_dimension=2, metric_warning=( "Dit is de grondoppervlakte van gebouwcontouren, niet de totale vloeroppervlakte " "of het gebouwvolume." ), limitation_message=( "GRB GBG bevat gebouwcontouren uit de basiskaart. Registratie en fysieke verandering " "kunnen in tijd verschillen." ), ), GrbProduct( key="roads", display_name="GRB wegsegmenten", reference_layer_name="roads", layer_type="road", collections=(GrbCollection("Wegsegment", 1),), geometry_types=("LineString", "MultiLineString"), metric_key="road_length", metric_method="intersection_length", metric_label="Totale weglengte", metric_unit="km", metric_dimension=1, metric_warning=( "De lengte volgt GRB-wegsegmenten en zegt niets over rijstroken, verkeersvolume " "of verhardingsoppervlakte." ), limitation_message=( "GRB Wegsegment beschrijft netwerkgeometrie en is geen verkeersmodel of routeadvies." ), ), GrbProduct( key="water", display_name="GRB wateroppervlakken en waterlijnen", reference_layer_name="water", layer_type="water", collections=( GrbCollection("WTZ", 2), GrbCollection("WLAS", 1), GrbCollection("WGR", 1), ), geometry_types=("LineString", "MultiLineString", "Polygon", "MultiPolygon"), metric_key="water_area", metric_method="intersection_area", metric_label="Wateroppervlakte", metric_unit="ha", metric_dimension=2, metric_warning=( "Watervolume is niet berekenbaar zonder betrouwbare diepte- of bathymetrische gegevens. " "De GRB-bron levert alleen oppervlakte- en lijngeometrie." ), limitation_message=( "GRB-water combineert wateroppervlakken en watergerelateerde lijnen. Objectaantallen en " "oppervlakte zijn geen actueel waterpeil of watervolume." ), ), GrbProduct( key="parcels", display_name="GRB administratieve percelen", reference_layer_name="parcels", layer_type="parcel", collections=(GrbCollection("ADP", 2),), geometry_types=("Polygon", "MultiPolygon"), metric_key="parcel_area", metric_method="intersection_area", metric_label="Perceeloppervlakte", metric_unit="ha", metric_dimension=2, metric_warning=( "GRB-percelen zijn een grafische referentie en vormen geen juridische grensopmeting." ), limitation_message=( "GRB ADP toont de vermoedelijke ligging van kadastrale percelen en is geen juridische grens." ), ), ) return {product.key: product for product in products} @staticmethod def list_products() -> list[dict[str, Any]]: return [ GrbProductRead( key=product.key, display_name=product.display_name, reference_layer_name=product.reference_layer_name, collections=[collection.name for collection in product.collections], geometry_types=list(product.geometry_types), source_crs=GrbAcquisitionService.SOURCE_CRS, authority_level=GrbAcquisitionService.AUTHORITY_LEVEL, catalog_url=GrbAcquisitionService.CATALOG_URL, attribution=GrbAcquisitionService.ATTRIBUTION, license_note=GrbAcquisitionService.LICENSE_NOTE, limitation_message=product.limitation_message, ).model_dump() for product in GrbAcquisitionService._products().values() ] @staticmethod def _product(product_key: str) -> GrbProduct: product = GrbAcquisitionService._products().get(product_key.strip().lower()) if product is None: raise AppError( code="GRB_PRODUCT_NOT_SUPPORTED", message="Select buildings, roads, water or parcels from the governed GRB product registry", details={"product_key": product_key}, status_code=422, ) return product @staticmethod def _validate_scope( db, project_id: UUID, payload: GrbAcquireRequest, settings: Settings, ) -> tuple[Any, list[float], list[float]]: if not settings.grb_enabled: raise AppError(code="GRB_NOT_CONFIGURED", message="Bounded GRB acquisition is disabled", status_code=503) 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="GRB_INVALID_CRS", message="GRB 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): raise AppError(code="GRB_INVALID_BBOX", message="Bounding box values must be finite", status_code=400) if values[0] >= values[2] or values[1] >= values[3]: raise AppError(code="GRB_INVALID_BBOX", message="Bounding box has no area", status_code=400) if values[0] < -180 or values[2] > 180 or values[1] < -90 or values[3] > 90: raise AppError(code="GRB_INVALID_BBOX", message="Bounding box is outside EPSG:4326", status_code=400) transformer = Transformer.from_crs("EPSG:4326", "EPSG:31370", always_xy=True) metric_bounds = transformer.transform_bounds(*values, densify_pts=21) width_m = float(metric_bounds[2] - metric_bounds[0]) height_m = float(metric_bounds[3] - metric_bounds[1]) if width_m < settings.grb_min_side_m or height_m < settings.grb_min_side_m: raise AppError( code="GRB_SELECTION_TOO_SMALL", message=f"Select an area of at least {settings.grb_min_side_m:g} by {settings.grb_min_side_m:g} metres", status_code=422, ) if width_m > settings.grb_max_side_m or height_m > settings.grb_max_side_m: raise AppError( code="GRB_SELECTION_TOO_LARGE", message=f"Select an area no larger than {settings.grb_max_side_m:g} by {settings.grb_max_side_m:g} metres", details={"width_m": width_m, "height_m": height_m}, status_code=422, ) selection = box(*values) if payload.area_id is None: scope_geometry = selection else: area = db.get(Area, payload.area_id) if area is None: 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) scope_geometry = to_shape(area.geometry).intersection(selection) if scope_geometry.is_empty: raise AppError( code="GRB_SCOPE_EMPTY", message="The requested bounding box does not intersect the selected area", status_code=400, ) return scope_geometry, [float(value) for value in values], [float(value) for value in metric_bounds] @staticmethod def _geometry_dimension(geometry: Any) -> int: if geometry is None or geometry.is_empty: return -1 if "Polygon" in geometry.geom_type: return 2 if "LineString" in geometry.geom_type or geometry.geom_type == "LinearRing": return 1 if "Point" in geometry.geom_type: return 0 if hasattr(geometry, "geoms"): return max((GrbAcquisitionService._geometry_dimension(item) for item in geometry.geoms), default=-1) return -1 @staticmethod def _extract_dimension(geometry: Any, expected_dimension: int) -> Any | None: if geometry is None or geometry.is_empty: return None if not geometry.is_valid: geometry = make_valid(geometry) parts: list[Any] = [] def collect(candidate: Any) -> None: if candidate is None or candidate.is_empty: return if expected_dimension == 2: if isinstance(candidate, Polygon): parts.append(candidate) return if isinstance(candidate, MultiPolygon): parts.extend(item for item in candidate.geoms if not item.is_empty) return if expected_dimension == 1: if isinstance(candidate, LineString): parts.append(candidate) return if isinstance(candidate, MultiLineString): parts.extend(item for item in candidate.geoms if not item.is_empty) return if hasattr(candidate, "geoms"): for item in candidate.geoms: collect(item) collect(geometry) if not parts: return None normalized = unary_union(parts) if normalized.is_empty: return None if not normalized.is_valid: normalized = make_valid(normalized) if ( normalized.is_empty or not normalized.is_valid or GrbAcquisitionService._geometry_dimension(normalized) != expected_dimension ): return None return normalized @staticmethod def _collection_url(settings: Settings, collection: GrbCollection, bbox_values: tuple[float, ...]) -> str: base = settings.grb_ogc_api_url.rstrip("/") query = urlencode( { "f": "application/geo+json", "limit": str(settings.grb_page_size), "bbox": ",".join(f"{value:.8f}" for value in bbox_values), "bbox-crs": GrbAcquisitionService.OGC_CRS84_URI, "crs": GrbAcquisitionService.OGC_CRS84_URI, } ) return f"{base}/collections/{collection.name}/items?{query}" @staticmethod def _validated_page_url( url: str, settings: Settings, product: GrbProduct, bbox_values: tuple[float, ...], ) -> str: parsed = urlparse(url) base = urlparse(settings.grb_ogc_api_url) allowed_paths = { f"{base.path.rstrip('/')}/collections/{collection.name}/items" for collection in product.collections } if ( parsed.scheme != "https" or base.scheme != "https" or parsed.netloc.casefold() != base.netloc.casefold() or parsed.path not in allowed_paths ): raise AppError( code="GRB_PROVIDER_INVALID_PAGINATION", message="GRB returned a pagination URL outside the governed OGC API allowlist", status_code=502, ) query = dict(parse_qsl(parsed.query, keep_blank_values=True)) query.update( { "f": "application/geo+json", "limit": str(settings.grb_page_size), "bbox": ",".join(f"{value:.8f}" for value in bbox_values), "bbox-crs": GrbAcquisitionService.OGC_CRS84_URI, "crs": GrbAcquisitionService.OGC_CRS84_URI, } ) return urlunparse((parsed.scheme, parsed.netloc, parsed.path, "", urlencode(query), "")) @staticmethod def _read_page( url: str, settings: Settings, opener: Callable[..., Any] | None, ) -> tuple[dict[str, Any], str, int]: request = Request( url, headers={ "Accept": "application/geo+json, application/json", "User-Agent": "GeoIntel/1.0 bounded-grb-acquisition", }, ) try: with (opener or guarded_opener(url, allow_redirect=False))( request, timeout=settings.grb_timeout_seconds, ) as response: limit = settings.grb_max_response_mb * 1024 * 1024 content = response.read(limit + 1) except HTTPError as exc: raise AppError( code="GRB_PROVIDER_HTTP_ERROR", message="The GRB OGC API returned an HTTP error", details={"status_code": exc.code}, status_code=502, ) from exc except (TimeoutError, URLError, OSError) as exc: raise AppError( code="GRB_PROVIDER_UNAVAILABLE", message="The GRB OGC API is unavailable", status_code=502, ) from exc if len(content) > limit: raise AppError( code="GRB_PROVIDER_RESPONSE_TOO_LARGE", message="A GRB response page exceeded the configured size limit", status_code=502, ) try: payload = json.loads(content.decode("utf-8")) except (UnicodeDecodeError, json.JSONDecodeError) as exc: raise AppError( code="GRB_PROVIDER_INVALID_RESPONSE", message="The GRB OGC API returned invalid GeoJSON", status_code=502, ) from exc if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection": raise AppError( code="GRB_PROVIDER_INVALID_RESPONSE", message="The GRB OGC API returned a non-FeatureCollection response", status_code=502, ) return payload, hashlib.sha256(content).hexdigest(), len(content) @staticmethod def _next_url(payload: dict[str, Any], current_url: str) -> str | None: links = payload.get("links") if not isinstance(links, list): return None for link in links: if isinstance(link, dict) and link.get("rel") == "next" and link.get("href"): return urljoin(current_url, str(link["href"])) return None @staticmethod def _fetch_features( product: GrbProduct, scope_geometry: Any, bbox_values: tuple[float, ...], coverage_scope: str, settings: Settings, opener: Callable[..., Any] | None, ) -> tuple[list[dict[str, Any]], dict[str, Any]]: retained: list[dict[str, Any]] = [] seen_ids: set[str] = set() request_urls: list[str] = [] response_sha256: list[str] = [] candidate_feature_count = 0 total_response_bytes = 0 geometry_types: dict[str, int] = {} collection_counts: dict[str, int] = {item.name: 0 for item in product.collections} for collection in product.collections: url: str | None = GrbAcquisitionService._collection_url(settings, collection, bbox_values) seen_pages: set[str] = set() while url: url = GrbAcquisitionService._validated_page_url(url, settings, product, bbox_values) if url in seen_pages: raise AppError( code="GRB_PROVIDER_PAGINATION_LOOP", message="The GRB OGC API repeated a pagination URL", status_code=502, ) if len(request_urls) >= settings.grb_max_pages: raise AppError( code="GRB_SELECTION_TOO_LARGE", message="GRB acquisition exceeded the configured page limit", details={"max_pages": settings.grb_max_pages}, status_code=422, ) seen_pages.add(url) payload, response_hash, response_size = GrbAcquisitionService._read_page(url, settings, opener) request_urls.append(url) response_sha256.append(response_hash) total_response_bytes += response_size if total_response_bytes > settings.grb_max_total_response_mb * 1024 * 1024: raise AppError( code="GRB_PROVIDER_RESPONSE_TOO_LARGE", message="The complete GRB response exceeded the configured transfer limit", status_code=502, ) source_features = payload.get("features") if not isinstance(source_features, list): raise AppError( code="GRB_PROVIDER_INVALID_RESPONSE", message="The GRB FeatureCollection has no valid feature list", status_code=502, ) for source_feature in source_features: candidate_feature_count += 1 if not isinstance(source_feature, dict): continue raw_id = str(source_feature.get("id") or "").strip() if not raw_id: raise AppError( code="GRB_PROVIDER_INVALID_RESPONSE", message=f"GRB {collection.name} returned a feature without an official identity", status_code=502, ) feature_id = f"{collection.name}:{raw_id}" if feature_id in seen_ids: continue seen_ids.add(feature_id) try: source_geometry = GrbAcquisitionService._extract_dimension( shape(source_feature.get("geometry")), collection.geometry_dimension, ) except Exception as exc: raise AppError( code="GRB_PROVIDER_INVALID_GEOMETRY", message=f"GRB {collection.name} returned invalid geometry", status_code=502, ) from exc if source_geometry is None or not source_geometry.intersects(scope_geometry): continue retained_geometry = GrbAcquisitionService._extract_dimension( source_geometry.intersection(scope_geometry), collection.geometry_dimension, ) if retained_geometry is None: continue if len(retained) >= settings.grb_max_features: raise AppError( code="GRB_SELECTION_TOO_LARGE", message="GRB selection exceeds the configured feature limit; draw a smaller rectangle", details={"max_features": settings.grb_max_features}, status_code=422, ) properties = dict(source_feature.get("properties") or {}) properties.update( { "source_name": GrbAcquisitionService.PROVIDER, "source_collection": collection.name, "source_feature_id": feature_id, "reference_layer_name": product.reference_layer_name, "layer_type": product.layer_type, "theme": product.key, "authority_level": GrbAcquisitionService.AUTHORITY_LEVEL, "coverage_scope": coverage_scope, "geometry_clipped_to_selection": not scope_geometry.covers(source_geometry), "attribution": GrbAcquisitionService.ATTRIBUTION, } ) retained.append( { "type": "Feature", "id": feature_id, "geometry": mapping(retained_geometry), "properties": properties, } ) collection_counts[collection.name] += 1 geometry_types[retained_geometry.geom_type] = geometry_types.get(retained_geometry.geom_type, 0) + 1 url = GrbAcquisitionService._next_url(payload, url) return retained, { "candidate_feature_count": candidate_feature_count, "feature_count": len(retained), "page_count": len(request_urls), "request_urls": request_urls, "response_sha256": response_sha256, "response_size_bytes": total_response_bytes, "collection_feature_counts": collection_counts, "geometry_types": geometry_types, "output_crs": GrbAcquisitionService.OGC_CRS84_URI, "reference_truncated": False, } @staticmethod def _cached_dataset( db, project_id: UUID, product: GrbProduct, request_hash: str, settings: Settings, ) -> Dataset | None: if settings.grb_cache_ttl_hours <= 0: return None candidates = ( db.query(Dataset) .filter( Dataset.project_id == project_id, Dataset.source_name == GrbAcquisitionService.PROVIDER, Dataset.reference_layer_name == product.reference_layer_name, Dataset.status == "ready", ) .order_by(Dataset.imported_at.desc()) .all() ) cutoff = datetime.now(UTC) - timedelta(hours=settings.grb_cache_ttl_hours) for candidate in candidates: provenance = candidate.provenance_metadata if isinstance(candidate.provenance_metadata, dict) else {} imported_at = candidate.imported_at if ( provenance.get("request_hash") == request_hash and candidate.storage_path and Path(candidate.storage_path).is_file() and imported_at is not None and imported_at >= cutoff ): return candidate return None @staticmethod def _result( dataset: Dataset, product: GrbProduct, *, reused: bool, bbox_values: list[float], ) -> dict[str, Any]: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} provenance = dataset.provenance_metadata if isinstance(dataset.provenance_metadata, dict) else {} metadata = dataset.metadata_json if isinstance(dataset.metadata_json, dict) else {} return GrbAcquisitionResult( output_dataset_id=dataset.id, reused=reused, provider=GrbAcquisitionService.PROVIDER, product_key=product.key, display_name=product.display_name, reference_layer_name=product.reference_layer_name, collections=[collection.name for collection in product.collections], feature_count=int(metadata.get("feature_count", source_metadata.get("feature_count", 0))), candidate_feature_count=int(provenance.get("candidate_feature_count", 0)), page_count=int(provenance.get("page_count", 0)), bbox_epsg4326=bbox_values, source_version=str(dataset.source_version or ""), attribution=GrbAcquisitionService.ATTRIBUTION, limitation_message=product.limitation_message, ).model_dump(mode="json") @staticmethod def acquire( db, project_id: UUID, payload: GrbAcquireRequest, *, settings: Settings | None = None, opener: Callable[..., Any] | None = None, ) -> dict[str, Any]: resolved_settings = settings or get_settings() product = GrbAcquisitionService._product(payload.product_key) scope_geometry, bbox_values, metric_bounds = GrbAcquisitionService._validate_scope( db, project_id, payload, resolved_settings, ) request_identity = { "provider": GrbAcquisitionService.PROVIDER, "product_key": product.key, "bbox_epsg4326": [round(value, 8) for value in bbox_values], "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() if not payload.force_refresh: cached = GrbAcquisitionService._cached_dataset( db, project_id, product, request_hash, resolved_settings, ) if cached is not None: return GrbAcquisitionService._result(cached, product, reused=True, bbox_values=bbox_values) area = db.get(Area, payload.area_id) if payload.area_id else None coverage_scope = ( "municipality" if area is not None and area.name.strip().lower().startswith("gemeente ") else "bounded_selection" ) features, transfer = GrbAcquisitionService._fetch_features( product, scope_geometry, tuple(scope_geometry.bounds), coverage_scope, resolved_settings, opener, ) acquired_at = datetime.now(UTC) source_version = acquired_at.date().isoformat() feature_collection = { "type": "FeatureCollection", "name": product.display_name, "crs": {"type": "name", "properties": {"name": "EPSG:4326"}}, "features": features, } artifact = json.dumps(feature_collection, ensure_ascii=False, separators=(",", ":")).encode("utf-8") filename = f"grb_{product.key}_{source_version}_{request_hash[:12]}.geojson" source_metadata: dict[str, Any] = { "provider": "Digitaal Vlaanderen", "service": "OGC API Features", "product_key": product.key, "product_display_name": product.display_name, "collections": [collection.name for collection in product.collections], "authority_level": GrbAcquisitionService.AUTHORITY_LEVEL, "theme": product.key, "layer_type": product.layer_type, "coverage_scope": coverage_scope, "geometry_clipped_to_area": payload.area_id is not None, "geometry_clipped_to_selection": True, "bbox_epsg4326": bbox_values, "bbox_epsg31370": metric_bounds, "feature_count": len(features), "collection_feature_counts": transfer["collection_feature_counts"], "identity_stable": True, "identity_scheme": "grb_ogc_feature_id", "source_storage_crs": "EPSG:31370", "requested_output_crs": GrbAcquisitionService.OGC_CRS84_URI, "selection_aggregation": { "metric_key": product.metric_key, "method": product.metric_method, "label": product.metric_label, "unit": product.metric_unit, "geometry_dimension": product.metric_dimension, "is_estimate": False, "warning": product.metric_warning, }, "attribution": GrbAcquisitionService.ATTRIBUTION, "license_note": GrbAcquisitionService.LICENSE_NOTE, "catalog_url": GrbAcquisitionService.CATALOG_URL, "limitation_message": product.limitation_message, } if product.key == "water": source_metadata["selection_metrics"] = [ { "metric_key": "water_length", "method": "intersection_length", "label": "Lengte watergerelateerde lijnen", "unit": "km", "geometry_dimension": 1, "is_estimate": False, } ] try: dataset_response = DatasetService.import_vector_bytes( db, project_id=project_id, area_id=payload.area_id, filename=filename, content=artifact, source="Digitaal Vlaanderen GRB OGC API Features", source_name=GrbAcquisitionService.PROVIDER, dataset_role="reference", reference_layer_name=product.reference_layer_name, temporal_series_key=f"grb:{product.key}:{request_hash[:24]}", observed_at=datetime( acquired_at.year, acquired_at.month, acquired_at.day, tzinfo=UTC, ), temporal_granularity="snapshot", source_version=source_version, source_metadata=source_metadata, provenance_metadata={ "acquisition": "explicit_bounded_ogc_api_features", "acquired_at": acquired_at.isoformat(), "request_hash": request_hash, "request_urls": transfer["request_urls"], "response_sha256": transfer["response_sha256"], "response_size_bytes": transfer["response_size_bytes"], "page_count": transfer["page_count"], "candidate_feature_count": transfer["candidate_feature_count"], "exact_feature_count": transfer["feature_count"], "reference_truncated": False, "artifact_sha256": hashlib.sha256(artifact).hexdigest(), "clipped_to_area_id": str(payload.area_id) if payload.area_id else None, "scope_geometry_type": scope_geometry.geom_type, "limitation_message": product.limitation_message, }, ) except AppError: raise except Exception as exc: raise AppError( code="GRB_PERSISTENCE_FAILED", message="The validated GRB selection could not be persisted", details={"reason": str(exc)}, status_code=500, ) from exc persisted = db.get(Dataset, dataset_response.id) if persisted is None: raise AppError( code="GRB_PERSISTENCE_FAILED", message="The persisted GRB dataset could not be reloaded", status_code=500, ) return GrbAcquisitionService._result(persisted, product, reused=False, bbox_values=bbox_values)