from __future__ import annotations import json import math from pathlib import Path from typing import Any, Iterable from uuid import UUID from geoalchemy2.functions import ST_Intersects, ST_MakeEnvelope from geoalchemy2.shape import from_shape from geoalchemy2.shape import to_shape from pyproj import CRS, Transformer from shapely.geometry import box, mapping, shape from shapely.ops import transform as transform_geometry from shapely.validation import make_valid from sqlalchemy import Float, String, case, cast, func from app.core.errors import AppError from app.models import Dataset, VectorFeature FULL_AREA_CLIPPED_OPERATOR_TOOLS = { "provision_mol_population_history.py", "provision_official_landuse_timeseries.py", "provision_regional_grb_buildings.py", "provision_regional_grb_context.py", "provision_regional_historical_landuse.py", "provision_waterinfo_station_history.py", "provision_mol_bwk_natura2000.py", "provision_regional_bwk_natura2000.py", "provision_agricultural_parcel_history.py", "provision_buildings_addresses_register.py", "provision_mol_soil_map.py", } SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS = { # Historical land-use themes are polygon map classes. Generic live-theme # line metrics (road/watercourse length) would therefore be meaningless. "provision_regional_historical_landuse.py", } PROPERTY_AGGREGATION_METHODS = {"sum", "mean", "area_weighted_sum"} PROPERTY_EXTREMA_METHODS = {"min", "max"} PRECLIPPED_MUNICIPALITY_PARTITION_OPERATOR_TOOLS = { "provision_regional_bwk_natura2000.py", } SEMANTIC_SELECTION_METRICS: dict[str, tuple[dict[str, Any], ...]] = { "administrative": ( { "metric_key": "covered_area", "method": "intersection_area", "label": "Bestuurlijk ingedeelde oppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": ( "Dit is de doorsnede met één bestuurlijk schaalniveau uit de gekozen NGI-laag; " "het is geen kadastrale of juridische grensopmeting." ), }, ), "buildings": ( { "metric_key": "footprint_area", "method": "intersection_area", "label": "Bebouwde grondoppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": "Dit is de grondoppervlakte van gebouwcontouren, niet de totale vloeroppervlakte of het gebouwvolume.", }, ), "forest": ( { "metric_key": "forest_area", "method": "intersection_area", "label": "Bosoppervlakte", "unit": "ha", "geometry_dimension": 2, }, ), "water": ( { "metric_key": "water_area", "method": "intersection_area", "label": "Wateroppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": "Watervolume is niet berekenbaar zonder betrouwbare diepte- of bathymetrische gegevens. De kaartbron levert alleen oppervlakte- en lijngeometrie.", }, { "metric_key": "watercourse_length", "method": "intersection_length", "label": "Lengte waterlopen", "unit": "km", "geometry_dimension": 1, }, ), "roads": ( { "metric_key": "road_length", "method": "intersection_length", "label": "Totale weglengte", "unit": "km", "geometry_dimension": 1, "warning": "De lengte volgt de GRB-wegsegmenten en zegt niets over rijstroken, verkeersvolume of verhardingsoppervlakte.", }, ), "parcels": ( { "metric_key": "parcel_area", "method": "intersection_area", "label": "Perceeloppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": "GRB-percelen zijn een grafische referentie en vormen geen juridische grensopmeting.", }, ), "nature_value": (), "agriculture": (), "soil": ( { "metric_key": "soil_mapped_area", "method": "intersection_area", "label": "Bodemkaartoppervlakte", "unit": "ha", "geometry_dimension": 2, "warning": "Historische bodemkartering op schaal 1:20.000; actuele lokale bodem- en drainagetoestand kan afwijken.", }, ), # Maritieme plan- en rapportagezones kunnen elkaar overlappen. Een # opgetelde oppervlakte zou daarom geen unieke zeeoppervlakte voorstellen. "maritime_planning": (), "marine_environment": (), } SEMANTIC_COUNT_LABELS = { "administrative": "Bestuursgebieden", "buildings": "Gebouwen", "population": "Statistische sectoren", "forest": "Bosvlakken", "water": "Waterobjecten", "roads": "Wegsegmenten", "parcels": "Percelen", "nature_value": "BWK-kaartvlakken", "agriculture": "Landbouwgebruikspercelen", "soil": "Bodemkaartvlakken", "maritime_planning": "Maritieme planobjecten", "marine_environment": "Mariene rapportagezones", } # Sprint 205 initially normalized two official comma-separated ALZ group labels # mechanically. Keep those persisted values queryable while new artifacts use # the explicit controlled keys. SELECTION_FILTER_VALUE_ALIASES: dict[tuple[str, str], tuple[str, ...]] = { ("main_crop_group_key", "grains_seeds_legumes"): ("granen,_zaden_en_peulvruchten",), ("main_crop_group_key", "horticulture"): ("groenten,_kruiden_en_sierplanten",), } class VectorFeatureService: MAX_VECTOR_PARTITIONS = 500 @staticmethod def _expanded_selection_filter_values(filter_property: str, filter_values: list[Any]) -> list[str]: expanded: list[str] = [] for value in filter_values: normalized = str(value) expanded.append(normalized) expanded.extend(SELECTION_FILTER_VALUE_ALIASES.get((filter_property, normalized), ())) return list(dict.fromkeys(expanded)) @staticmethod def _dataset_theme(dataset: Dataset) -> str | None: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} candidates = ( source_metadata.get("theme"), dataset.reference_layer_name, source_metadata.get("layer_type"), ) aliases = { "belgium_land_boundary": "administrative", "belgium_regions": "administrative", "belgium_provinces": "administrative", "belgium_municipalities": "administrative", "marine_spatial_plan_2026": "maritime_planning", "marine_legal_scopes": "marine_environment", "building": "buildings", "bebouwing": "buildings", "population": "population", "forest": "forest", "forestry": "forest", "waterways": "water", "road": "roads", "parcel": "parcels", "nature": "nature_value", "biodiversity": "nature_value", "bwk": "nature_value", "natura2000": "nature_value", "agricultural": "agriculture", "landbouw": "agriculture", "landbouwgebruik": "agriculture", "building_registry": "buildings", "soil_map": "soil", "bodem": "soil", } for candidate in candidates: if not isinstance(candidate, str) or not candidate.strip(): continue normalized = candidate.strip().lower() if normalized.startswith("regional_"): normalized = normalized.removeprefix("regional_") normalized = aliases.get(normalized, normalized) if normalized in {*SEMANTIC_SELECTION_METRICS, "population"}: return normalized return None @staticmethod def supports_selection_summary(dataset: Dataset) -> bool: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} return isinstance(source_metadata.get("selection_aggregation"), dict) or VectorFeatureService._dataset_theme(dataset) is not None @staticmethod def deduplicate_rows(rows: list[Any]) -> list[Any]: """Collapse rows that describe one source feature across partitions. Municipal partitions of one product overlap at their shared boundary, so a rectangle drawn across it returns the same building from both. An empty or missing ``source_feature_id`` is not a shared identity — two rows without one are two features, not a duplicate pair. """ seen: set[str] = set() kept: list[Any] = [] for row in rows: source_feature_id = getattr(row, "source_feature_id", None) identity = str(source_feature_id).strip() if source_feature_id is not None else "" if not identity: kept.append(row) continue if identity in seen: continue seen.add(identity) kept.append(row) return kept @staticmethod def count_disclosure( *, total_feature_count: int, fully_covered_feature_count: int | None, ) -> dict[str, Any]: """Describe how much of the counted population the selection cuts. A feature that merely touches the drawn rectangle is counted whole, while ``intersection_area`` clips it. Reporting both numbers without saying so puts two figures for different populations side by side. The count stays whole-feature — that is what an operator expects from "objecten" — but says how many of them the edge cuts, and is flagged as an estimate when it does. ``fully_covered_feature_count`` is ``None`` when the selection covers a pre-clipped whole work area, where no edge effect exists. """ if fully_covered_feature_count is None: return { "partially_covered_feature_count": None, "is_estimate": False, "warning": None, } partial = max(0, int(total_feature_count) - int(fully_covered_feature_count)) if partial <= 0: return { "partially_covered_feature_count": 0, "is_estimate": False, "warning": None, } return { "partially_covered_feature_count": partial, "is_estimate": True, "warning": ( f"{partial} van de {int(total_feature_count)} objecten liggen deels buiten de selectie en zijn " "aan de rand doorgesneden. Ze tellen volledig mee in het aantal; oppervlakte- en lengtematen " "gebruiken alleen het deel binnen de selectie." ), } @staticmethod def constrain_bbox_to_area( bbox: dict[str, Any], area_geometry: Any, ) -> tuple[Any, bool]: bbox_geometry = box( float(bbox["min_x"]), float(bbox["min_y"]), float(bbox["max_x"]), float(bbox["max_y"]), ) area_shape = to_shape(area_geometry) constrained_geometry = bbox_geometry.intersection(area_shape) if constrained_geometry.is_empty or constrained_geometry.area <= 0: raise AppError( code="VECTOR_SELECTION_OUTSIDE_AREA", message="Selection does not overlap the selected work area", status_code=422, ) return from_shape(constrained_geometry, srid=4326), constrained_geometry.equals(area_shape) @staticmethod def can_use_full_area_fast_path(dataset: Dataset, selection_area_id: UUID | None) -> bool: if selection_area_id is None or dataset.area_id != selection_area_id: return False source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} if source_metadata.get("geometry_clipped_to_area") is True: return True provenance = dataset.provenance_metadata if isinstance(dataset.provenance_metadata, dict) else {} return provenance.get("operator_tool") in FULL_AREA_CLIPPED_OPERATOR_TOOLS @staticmethod def preclipped_partition_filter(dataset: Dataset, selection_area_name: str | None) -> tuple[str, str] | None: provenance = dataset.provenance_metadata if isinstance(dataset.provenance_metadata, dict) else {} if provenance.get("operator_tool") not in PRECLIPPED_MUNICIPALITY_PARTITION_OPERATOR_TOOLS: return None source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} if ( source_metadata.get("partitioned_source_audit") is not True or source_metadata.get("geometry_clipped_to_area") is not True ): return None normalized_name = str(selection_area_name or "").strip() prefix = "Gemeente " if not normalized_name.startswith(prefix): return None municipality = normalized_name[len(prefix):].split(" - ", 1)[0].strip() return ("municipality", municipality) if municipality else None @staticmethod def _feature_row( dataset_id: UUID, feature: dict[str, Any], index: int, feature_class: str | None, *, source_crs: str = "EPSG:4326", ) -> VectorFeature | None: geometry_payload = feature.get("geometry") if geometry_payload is None: return None try: geometry = shape(geometry_payload) except Exception as exc: raise AppError(code="INVALID_GEOJSON", message=f"Invalid feature geometry at index {index}", status_code=400) from exc if geometry.is_empty: return None if not geometry.is_valid: geometry = make_valid(geometry) if geometry.is_empty or not geometry.is_valid: raise AppError(code="INVALID_GEOMETRY", message=f"Invalid feature geometry at index {index}", status_code=400) geometry = VectorFeatureService._canonical_geometry(geometry, source_crs=source_crs, index=index) properties = feature.get("properties") if isinstance(feature.get("properties"), dict) else {} source_feature_id = feature.get("id") if source_feature_id is None: source_feature_id = properties.get("id") or properties.get("source_feature_id") return VectorFeature( dataset_id=dataset_id, feature_class=feature_class, source_feature_id=str(source_feature_id) if source_feature_id is not None else None, properties_json=properties, geometry=from_shape(geometry, srid=4326), ) @staticmethod def _canonical_geometry(geometry: Any, *, source_crs: str, index: int): """Transform one source geometry to canonical EPSG:4326 safely.""" if geometry.has_z: geometry = transform_geometry(lambda x, y, z=None: (x, y), geometry) # VectorFeature is deliberately canonical WGS84 storage. Treating # Lambert or another source CRS as EPSG:4326 produces geometries that # look syntactically valid but are spatially wrong. All governed # import callers therefore pass the declared source CRS; the default # only preserves compatibility for legacy, already-WGS84 call sites. try: parsed_source_crs = CRS.from_user_input(source_crs) target_crs = CRS.from_epsg(4326) except Exception as exc: raise AppError( code="INVALID_DATASET_CRS", message=f"Invalid source CRS for vector feature at index {index}", details={"source_crs": source_crs}, status_code=400, ) from exc if not parsed_source_crs.equals(target_crs): try: transformer = Transformer.from_crs(parsed_source_crs, target_crs, always_xy=True) geometry = transform_geometry(transformer.transform, geometry) except Exception as exc: raise AppError( code="VECTOR_CRS_TRANSFORMATION_FAILED", message=f"Could not transform vector feature at index {index} to EPSG:4326", details={"source_crs": source_crs}, status_code=400, ) from exc if geometry.is_empty or not geometry.is_valid: geometry = make_valid(geometry) if geometry.is_empty or not geometry.is_valid: raise AppError( code="INVALID_GEOMETRY", message=f"Invalid transformed feature geometry at index {index}", status_code=400, ) min_x, min_y, max_x, max_y = geometry.bounds if ( not all(math.isfinite(value) for value in (min_x, min_y, max_x, max_y)) or min_x < -180 or max_x > 180 or min_y < -90 or max_y > 90 ): raise AppError( code="VECTOR_GEOMETRY_OUTSIDE_EPSG4326", message=f"Transformed feature geometry at index {index} is outside EPSG:4326 bounds", details={"source_crs": source_crs, "bounds": [min_x, min_y, max_x, max_y]}, status_code=400, ) return geometry @staticmethod def canonicalize_geojson_payload(payload: dict[str, Any], *, source_crs: str) -> dict[str, Any]: """Return a canonical-WGS84 feature collection without losing source attributes. Callers use this payload for validation, spatial indexing and map-safe consumption storage. A non-canonical source file, when retained, belongs to explicit provenance evidence rather than the Dataset consumption path; no implicit CRS assumption is recorded. """ features = payload.get("features") if payload.get("type") != "FeatureCollection" or not isinstance(features, list): raise AppError(code="INVALID_GEOJSON", message="GeoJSON payload must be a FeatureCollection", status_code=400) canonical_features: list[dict[str, Any]] = [] for index, feature in enumerate(features): if not isinstance(feature, dict): raise AppError(code="INVALID_GEOJSON", message=f"Feature {index} must be an object", status_code=400) canonical_feature = dict(feature) geometry_payload = feature.get("geometry") if geometry_payload is not None: try: geometry = shape(geometry_payload) except Exception as exc: raise AppError( code="INVALID_GEOJSON", message=f"Invalid feature geometry at index {index}", status_code=400, ) from exc if not geometry.is_empty: if not geometry.is_valid: geometry = make_valid(geometry) if geometry.is_empty or not geometry.is_valid: raise AppError( code="INVALID_GEOMETRY", message=f"Invalid feature geometry at index {index}", status_code=400, ) canonical_feature["geometry"] = mapping( VectorFeatureService._canonical_geometry(geometry, source_crs=source_crs, index=index) ) canonical_features.append(canonical_feature) return { **{key: value for key, value in payload.items() if key not in {"crs", "features"}}, "type": "FeatureCollection", "crs": {"type": "name", "properties": {"name": "EPSG:4326"}}, "features": canonical_features, } @staticmethod def _normalize_selection_bbox(bbox: dict[str, Any]) -> dict[str, float | str]: try: min_x = float(bbox["min_x"]) min_y = float(bbox["min_y"]) max_x = float(bbox["max_x"]) max_y = float(bbox["max_y"]) except (KeyError, TypeError, ValueError) as exc: raise AppError( code="INVALID_SELECTION_BBOX", message="Selection bbox must include numeric min_x, min_y, max_x and max_y values", status_code=400, ) from exc crs = str(bbox.get("crs") or "EPSG:4326").upper() if crs != "EPSG:4326": raise AppError( code="UNSUPPORTED_SELECTION_CRS", message="Map selection currently supports EPSG:4326 bbox coordinates only", details={"crs": crs}, status_code=400, ) if min_x >= max_x or min_y >= max_y: raise AppError( code="INVALID_SELECTION_BBOX", message="Selection bbox must have min_x < max_x and min_y < max_y", status_code=400, ) if min_x < -180 or max_x > 180 or min_y < -90 or max_y > 90: raise AppError( code="INVALID_SELECTION_BBOX", message="Selection bbox is outside EPSG:4326 longitude/latitude bounds", status_code=400, ) return {"min_x": min_x, "min_y": min_y, "max_x": max_x, "max_y": max_y, "crs": "EPSG:4326"} @staticmethod def _dataset_bbox_intersects( dataset: Dataset, bbox: dict[str, float | str], ) -> bool: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} source_bbox = source_metadata.get("bbox_epsg4326") if not isinstance(source_bbox, list) or len(source_bbox) != 4: return True try: min_x, min_y, max_x, max_y = (float(value) for value in source_bbox) except (TypeError, ValueError): return True return not ( max_x <= float(bbox["min_x"]) or min_x >= float(bbox["max_x"]) or max_y <= float(bbox["min_y"]) or min_y >= float(bbox["max_y"]) ) @staticmethod def _latest_complete_partition_manifest( datasets: Iterable[Dataset], *, source_name: str, partition_scope_key: str, ) -> list[Dataset]: groups: dict[str, list[Dataset]] = {} for dataset in datasets: source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} manifest_sha256 = str(source_metadata.get("partition_manifest_sha256") or "") if ( dataset.source_name != source_name or dataset.dataset_type not in {"vector", "geojson"} or dataset.status != "ready" or dataset.area_id is None or source_metadata.get("regional_partitions_complete") is not True or source_metadata.get("partition_scope_key") != partition_scope_key or len(manifest_sha256) != 64 ): continue groups.setdefault(manifest_sha256, []).append(dataset) complete_groups: list[list[Dataset]] = [] for group in groups.values(): area_ids = {dataset.area_id for dataset in group} expected_data_count = max( int((dataset.source_metadata or {}).get("data_partition_count") or 0) for dataset in group ) if expected_data_count > 0 and len(group) == expected_data_count and len(area_ids) == len(group): complete_groups.append(group) if not complete_groups: return [] def manifest_priority(group: list[Dataset]) -> tuple[str, int, str]: observed_at = max( str((dataset.source_metadata or {}).get("partition_manifest_observed_at") or "") for dataset in group ) manifest_sha256 = str((group[0].source_metadata or {}).get("partition_manifest_sha256") or "") return observed_at, len(group), manifest_sha256 selected = max(complete_groups, key=manifest_priority) return sorted(selected, key=lambda dataset: (str(dataset.area_id), str(dataset.id))) @staticmethod def select_partitioned_features_by_bbox( db, *, project_id: UUID, source_name: str, partition_scope_key: str, bbox: dict[str, Any], limit: int = 100, selection_geometry: Any | None = None, selection_area_id: UUID | None = None, partition_area_id: UUID | None = None, ) -> dict[str, Any]: normalized_bbox = VectorFeatureService._normalize_selection_bbox(bbox) safe_limit = max(1, min(int(limit), 1000)) project_datasets = ( db.query(Dataset) .filter( Dataset.project_id == project_id, Dataset.source_name == source_name, Dataset.status == "ready", ) .all() ) manifest_datasets = VectorFeatureService._latest_complete_partition_manifest( project_datasets, source_name=source_name, partition_scope_key=partition_scope_key, ) if not manifest_datasets: raise AppError( code="VECTOR_PARTITIONS_NOT_READY", message="No complete persisted vector partition manifest is available", details={"source_name": source_name, "partition_scope_key": partition_scope_key}, status_code=409, ) if len(manifest_datasets) > VectorFeatureService.MAX_VECTOR_PARTITIONS: raise AppError( code="VECTOR_PARTITION_LIMIT_EXCEEDED", message="The complete vector partition manifest exceeds the safety limit", details={ "partition_count": len(manifest_datasets), "max_partitions": VectorFeatureService.MAX_VECTOR_PARTITIONS, }, status_code=422, ) scoped_datasets = [ dataset for dataset in manifest_datasets if (partition_area_id is None or dataset.area_id == partition_area_id) and VectorFeatureService._dataset_bbox_intersects(dataset, normalized_bbox) ] dataset_ids = [dataset.id for dataset in scoped_datasets] representative = scoped_datasets[0] if scoped_datasets else manifest_datasets[0] selection_shape = selection_geometry if selection_shape is None: selection_shape = ST_MakeEnvelope( normalized_bbox["min_x"], normalized_bbox["min_y"], normalized_bbox["max_x"], normalized_bbox["max_y"], 4326, ) query = db.query(VectorFeature).filter( VectorFeature.dataset_id.in_(dataset_ids), ST_Intersects(VectorFeature.geometry, selection_shape), ) total_feature_count = int(query.count()) rows = ( query.order_by(VectorFeature.created_at.asc(), VectorFeature.id.asc()) .limit(safe_limit + 1) .all() ) features = [ VectorFeatureService._row_to_geojson_feature(row) for row in rows[:safe_limit] ] result = { "selection_bbox": normalized_bbox, "feature_count": len(features), "total_feature_count": total_feature_count, "limit": safe_limit, "truncated": total_feature_count > safe_limit, "geojson": {"type": "FeatureCollection", "features": features}, "partition_count": len(scoped_datasets), "available_partition_count": len(manifest_datasets), "partition_scope_key": partition_scope_key, "source_name": source_name, "dataset_ids": dataset_ids, } if selection_area_id is not None: result["selection_area_id"] = str(selection_area_id) if VectorFeatureService.supports_selection_summary(representative): result["summary"] = VectorFeatureService.summarize_features_by_bbox( db, dataset=representative, dataset_ids=dataset_ids, bbox=normalized_bbox, total_feature_count=total_feature_count, selection_geometry=selection_shape, ) return result @staticmethod def _row_to_geojson_feature(row: VectorFeature) -> dict[str, Any]: geometry_value = row.geometry try: geometry = geometry_value if hasattr(geometry_value, "__geo_interface__") else to_shape(geometry_value) except Exception as exc: raise AppError( code="INVALID_VECTOR_FEATURE_GEOMETRY", message="Persisted vector feature geometry could not be converted to GeoJSON", details={"vector_feature_id": str(row.id)}, status_code=500, ) from exc properties = dict(row.properties_json or {}) properties.update( { "vector_feature_id": str(row.id), "dataset_id": str(row.dataset_id), "source_feature_id": row.source_feature_id, "feature_class": row.feature_class, } ) return { "type": "Feature", "id": str(row.id), "geometry": mapping(geometry), "properties": properties, } @staticmethod def select_features_by_bbox( db, dataset_id: UUID, bbox: dict[str, Any], limit: int = 100, dataset: Dataset | None = None, selection_geometry: Any | None = None, selection_area_id: UUID | None = None, full_dataset_area: bool = False, preclipped_partition_filter: tuple[str, str] | None = None, dataset_ids: list[UUID] | None = None, deduplicate_source_features: bool = False, ) -> dict[str, Any]: normalized_bbox = VectorFeatureService._normalize_selection_bbox(bbox) safe_limit = max(1, min(int(limit), 1000)) selection_shape = selection_geometry if selection_shape is None: selection_shape = ST_MakeEnvelope( normalized_bbox["min_x"], normalized_bbox["min_y"], normalized_bbox["max_x"], normalized_bbox["max_y"], 4326, ) selected_dataset_ids = dataset_ids or [dataset_id] query = db.query(VectorFeature).filter(VectorFeature.dataset_id.in_(selected_dataset_ids)) if preclipped_partition_filter is not None: partition_property, partition_value = preclipped_partition_filter query = query.filter(VectorFeature.properties_json.op("->>")(partition_property) == partition_value) if not full_dataset_area: query = query.filter(ST_Intersects(VectorFeature.geometry, selection_shape)) if deduplicate_source_features: identity = func.coalesce(VectorFeature.source_feature_id, cast(VectorFeature.id, String)) total_feature_count = int( query.with_entities(func.count(func.distinct(identity))).scalar() or 0 ) elif hasattr(query, "count"): total_feature_count = int(query.count()) else: # Lightweight unit-test sessions do not always implement Query.count(). total_feature_count = len(query.all()) rows = ( query.order_by(VectorFeature.created_at.asc()) .limit(safe_limit + 1) .all() ) if deduplicate_source_features: # ``total_feature_count`` is already distinct; without this the map # would draw a boundary feature once per partition and the returned # count would exceed the headline number beside it. rows = VectorFeatureService.deduplicate_rows(rows) truncated = total_feature_count > safe_limit selected_rows = rows[:safe_limit] features = [VectorFeatureService._row_to_geojson_feature(row) for row in selected_rows] summary = None if dataset and VectorFeatureService.supports_selection_summary(dataset): summary = VectorFeatureService.summarize_features_by_bbox( db, dataset=dataset, dataset_ids=selected_dataset_ids, bbox=normalized_bbox, total_feature_count=total_feature_count, selection_geometry=selection_geometry, full_dataset_area=full_dataset_area, preclipped_partition_filter=preclipped_partition_filter, ) result = { "selection_bbox": normalized_bbox, "feature_count": len(features), "total_feature_count": total_feature_count, "limit": safe_limit, "truncated": truncated, "geojson": { "type": "FeatureCollection", "features": features, }, "summary": summary, } if selection_area_id is not None: result["selection_area_id"] = str(selection_area_id) return result @staticmethod def summarize_features_by_bbox( db, *, dataset: Dataset, bbox: dict[str, Any], dataset_ids: list[UUID] | None = None, total_feature_count: int | None = None, selection_geometry: Any | None = None, full_dataset_area: bool = False, preclipped_partition_filter: tuple[str, str] | None = None, disclose_selection_edge: bool = True, ) -> dict[str, Any]: normalized_bbox = VectorFeatureService._normalize_selection_bbox(bbox) selection_shape = selection_geometry if selection_shape is None: selection_shape = ST_MakeEnvelope( normalized_bbox["min_x"], normalized_bbox["min_y"], normalized_bbox["max_x"], normalized_bbox["max_y"], 4326, ) selection_filter = ( (VectorFeature.dataset_id.in_(dataset_ids),) if dataset_ids is not None else (VectorFeature.dataset_id == dataset.id,) ) if preclipped_partition_filter is not None: partition_property, partition_value = preclipped_partition_filter selection_filter += ( VectorFeature.properties_json.op("->>")(partition_property) == partition_value, ) if not full_dataset_area: selection_filter += (ST_Intersects(VectorFeature.geometry, selection_shape),) selection_is_preclipped = full_dataset_area feature_count = total_feature_count if feature_count is None: feature_count = int(db.query(func.count(VectorFeature.id)).filter(*selection_filter).scalar() or 0) # How many of the counted features the selection edge cuts. Skipped for # a pre-clipped whole-area selection, which has no edge to cut against. fully_covered_feature_count: int | None = None if disclose_selection_edge and not full_dataset_area and feature_count: try: fully_covered_feature_count = int( db.query(func.count(VectorFeature.id)) .filter(*selection_filter) .filter(func.ST_CoveredBy(VectorFeature.geometry, selection_shape)) .scalar() or 0 ) except Exception: # Lightweight unit-test sessions do not implement every spatial # predicate; the count then simply carries no edge disclosure. fully_covered_feature_count = None count_disclosure = VectorFeatureService.count_disclosure( total_feature_count=feature_count, fully_covered_feature_count=fully_covered_feature_count, ) source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {} config = source_metadata.get("selection_aggregation") if not isinstance(config, dict): config = {} theme = VectorFeatureService._dataset_theme(dataset) configured_metric = { "metric_key": str(config.get("metric_key") or config.get("method") or "feature_count"), "method": str(config.get("method") or "feature_count"), "label": str(config.get("label") or SEMANTIC_COUNT_LABELS.get(theme or "", "Objecten")), "unit": str(config.get("unit") or "objecten"), "warning": str(config["warning"]) if config.get("warning") else None, "is_estimate": bool(config.get("is_estimate", False)), **({"property": config.get("property")} if config.get("property") else {}), } provenance = dataset.provenance_metadata if isinstance(dataset.provenance_metadata, dict) else {} semantic_metrics_disabled = ( source_metadata.get("semantic_metrics") is False or provenance.get("operator_tool") in SEMANTIC_METRICS_DISABLED_OPERATOR_TOOLS ) semantic_metrics = ( [] if semantic_metrics_disabled else [dict(metric) for metric in SEMANTIC_SELECTION_METRICS.get(theme or "", ())] ) primary_config = configured_metric if configured_metric["method"] == "feature_count" and semantic_metrics: primary_config = semantic_metrics[0] metric_configs = [primary_config] configured_metrics = source_metadata.get("selection_metrics") if isinstance(configured_metrics, list): existing_metric_keys = {str(primary_config.get("metric_key") or "")} for configured_item in configured_metrics: if not isinstance(configured_item, dict): continue metric_key = str(configured_item.get("metric_key") or "").strip() if not metric_key or metric_key in existing_metric_keys: continue metric_configs.append(dict(configured_item)) existing_metric_keys.add(metric_key) for semantic_metric in semantic_metrics: signature = (semantic_metric["method"], semantic_metric["unit"]) existing = { (item["method"], item["unit"]) for item in metric_configs } if signature not in existing: metric_configs.append(semantic_metric) if not any(item["method"] == "feature_count" for item in metric_configs): metric_configs.append( { "metric_key": "feature_count", "method": "feature_count", "label": SEMANTIC_COUNT_LABELS.get(theme or "", "Objecten"), "unit": "objecten", } ) metrics = [ VectorFeatureService._calculate_selection_metric( db, dataset=dataset, config=metric_config, selection_filter=selection_filter, selection_shape=selection_shape, feature_count=feature_count, full_dataset_area=selection_is_preclipped, partially_covered_feature_count=count_disclosure["partially_covered_feature_count"], ) for metric_config in metric_configs ] # The whole-feature count carries the edge disclosure; a metric that # already clips to the selection (area, length) does not need it. for computed_metric in metrics: if computed_metric["aggregation_method"] == "feature_count" and count_disclosure["is_estimate"]: computed_metric["is_estimate"] = True computed_metric["warning"] = computed_metric.get("warning") or count_disclosure["warning"] primary_metric = metrics[0] return { "metric_label": primary_metric["metric_label"], "metric_value": primary_metric["metric_value"], "metric_unit": primary_metric["metric_unit"], "aggregation_method": primary_metric["aggregation_method"], "primary_metric_key": primary_metric["metric_key"], "feature_count": feature_count, "fully_covered_feature_count": fully_covered_feature_count, "partially_covered_feature_count": count_disclosure["partially_covered_feature_count"], "selection_edge_warning": count_disclosure["warning"], "is_estimate": primary_metric["is_estimate"], "warning": primary_metric.get("warning"), "metrics": metrics, } @staticmethod def _calculate_selection_metric( db, *, dataset: Dataset, config: dict[str, Any], selection_filter: tuple[Any, ...], selection_shape: Any, feature_count: int, full_dataset_area: bool, partially_covered_feature_count: int | None = None, ) -> dict[str, Any]: method = str(config.get("method") or "feature_count") unit = str(config.get("unit") or "objecten") warning = str(config["warning"]) if config.get("warning") else None is_estimate = bool(config.get("is_estimate", False)) metric_value = float(feature_count) dimension = config.get("geometry_dimension") metric_filter = selection_filter if dimension in {1, 2}: metric_filter += (func.ST_Dimension(VectorFeature.geometry) == int(dimension),) filter_property = str(config.get("filter_property") or "").strip() filter_values = config.get("filter_values") if filter_property: if not isinstance(filter_values, list) or not filter_values: raise AppError( code="INVALID_SELECTION_AGGREGATION", message="Dataset selection metric filter requires one or more values", details={"dataset_id": str(dataset.id), "filter_property": filter_property}, status_code=500, ) normalized_filter_values = VectorFeatureService._expanded_selection_filter_values( filter_property, filter_values, ) metric_filter += ( VectorFeature.properties_json.op("->>")(filter_property).in_(normalized_filter_values), ) if method == "intersection_area": source_area = func.ST_Area(func.ST_Transform(VectorFeature.geometry, 31370)) if full_dataset_area: area_expression = source_area else: covered_by_selection = func.ST_CoveredBy(VectorFeature.geometry, selection_shape) intersection_area = func.ST_Area( func.ST_Transform(func.ST_Intersection(VectorFeature.geometry, selection_shape), 31370) ) area_expression = case( (covered_by_selection, source_area), else_=intersection_area, ) area_m2 = db.query(func.coalesce(func.sum(area_expression), 0.0)).filter(*metric_filter).scalar() divisor = 10_000.0 if unit == "ha" else 1.0 metric_value = float(area_m2 or 0.0) / divisor elif method == "intersection_length": measured_geometry = ( VectorFeature.geometry if full_dataset_area else func.ST_Intersection(VectorFeature.geometry, selection_shape) ) length_expression = func.ST_Length(func.ST_Transform(measured_geometry, 31370)) length_m = db.query(func.coalesce(func.sum(length_expression), 0.0)).filter(*metric_filter).scalar() divisor = 1_000.0 if unit == "km" else 1.0 metric_value = float(length_m or 0.0) / divisor elif method in PROPERTY_AGGREGATION_METHODS | PROPERTY_EXTREMA_METHODS: property_name = str(config.get("property") or "").strip() if not property_name: raise AppError( code="INVALID_SELECTION_AGGREGATION", message="Dataset selection aggregation requires a numeric property", details={"dataset_id": str(dataset.id), "method": method}, status_code=500, ) numeric_value = cast(VectorFeature.properties_json.op("->>")(property_name), Float) value_expression = numeric_value covered_by_selection = None if method == "area_weighted_sum" and not full_dataset_area: source_area = func.ST_Area(func.ST_Transform(VectorFeature.geometry, 31370)) intersection_area = func.ST_Area( func.ST_Transform(func.ST_Intersection(VectorFeature.geometry, selection_shape), 31370) ) covered_by_selection = func.ST_CoveredBy(VectorFeature.geometry, selection_shape) coverage_ratio = case( (covered_by_selection, 1.0), else_=intersection_area / func.nullif(source_area, 0.0), ) value_expression = numeric_value * coverage_ratio aggregate_function = { "mean": func.avg, "min": func.min, "max": func.max, }.get(method, func.sum) aggregate_value = ( db.query(func.coalesce(aggregate_function(value_expression), 0.0)) .filter(*metric_filter) .filter(VectorFeature.properties_json.op("->>")(property_name).isnot(None)) .scalar() ) metric_value = float(aggregate_value or 0.0) if method == "area_weighted_sum" and not full_dataset_area: # The selection-edge count was already established for the # feature count; a second query would ask the same question. partial_feature_count = partially_covered_feature_count if partial_feature_count is None: partial_feature_count = ( db.query(func.count(VectorFeature.id)) .filter(*metric_filter) .filter(~covered_by_selection) .scalar() ) is_estimate = bool(config.get("is_estimate", False)) or bool(partial_feature_count) if not is_estimate and config.get("warning_only_when_estimate", True): warning = None elif method == "area_weighted_sum": is_estimate = bool(config.get("is_estimate", False)) if not is_estimate and config.get("warning_only_when_estimate", True): warning = None elif method == "feature_count" and (filter_property or dimension in {1, 2}): metric_value = float( db.query(func.count(VectorFeature.id)).filter(*metric_filter).scalar() or 0 ) elif method != "feature_count": raise AppError( code="INVALID_SELECTION_AGGREGATION", message="Unsupported dataset selection aggregation", details={"dataset_id": str(dataset.id), "method": method}, status_code=500, ) return { "metric_key": str(config.get("metric_key") or method), "metric_label": str(config.get("label") or "Objecten"), "metric_value": metric_value, "metric_unit": unit, "aggregation_method": method, "is_estimate": is_estimate, "warning": warning, } @staticmethod def persist_geojson_features( db, dataset_id: UUID, payload: dict[str, Any], feature_class: str | None = None, *, commit: bool = True, source_crs: str = "EPSG:4326", ) -> list[VectorFeature]: features = payload.get("features") if payload.get("type") != "FeatureCollection" or not isinstance(features, list): raise AppError(code="INVALID_GEOJSON", message="GeoJSON payload must be a FeatureCollection", status_code=400) persisted: list[VectorFeature] = [] for index, feature in enumerate(features): if not isinstance(feature, dict): raise AppError(code="INVALID_GEOJSON", message=f"Feature {index} must be an object", status_code=400) row = VectorFeatureService._feature_row( dataset_id, feature, index, feature_class, source_crs=source_crs, ) if row is None: continue db.add(row) persisted.append(row) if commit: db.flush() db.commit() return persisted @staticmethod def persist_geojson_partitions( db, dataset_id: UUID, partition_paths: Iterable[str | Path], feature_class: str | None = None, *, batch_size: int = 1000, source_crs: str = "EPSG:4326", ) -> int: if batch_size <= 0: raise ValueError("batch_size must be positive") persisted_count = 0 source_feature_ids: set[str] = set() for partition_path in partition_paths: path = Path(partition_path) try: payload = json.loads(path.read_text(encoding="utf-8")) except (OSError, json.JSONDecodeError) as exc: raise AppError( code="INVALID_GEOJSON_PARTITION", message=f"Could not read GeoJSON partition {path.name}", status_code=400, ) from exc features = payload.get("features") if payload.get("type") != "FeatureCollection" or not isinstance(features, list): raise AppError( code="INVALID_GEOJSON_PARTITION", message=f"GeoJSON partition {path.name} must be a FeatureCollection", status_code=400, ) batch: list[VectorFeature] = [] for index, feature in enumerate(features): if not isinstance(feature, dict): raise AppError( code="INVALID_GEOJSON_PARTITION", message=f"Feature {index} in {path.name} must be an object", status_code=400, ) row = VectorFeatureService._feature_row( dataset_id, feature, index, feature_class, source_crs=source_crs, ) if row is None: continue if row.source_feature_id: if row.source_feature_id in source_feature_ids: raise AppError( code="DUPLICATE_SOURCE_FEATURE", message=f"Duplicate source feature {row.source_feature_id} across regional partitions", status_code=400, ) source_feature_ids.add(row.source_feature_id) db.add(row) batch.append(row) persisted_count += 1 if len(batch) >= batch_size: db.flush() for persisted in batch: db.expunge(persisted) batch.clear() if batch: db.flush() for persisted in batch: db.expunge(persisted) return persisted_count