from __future__ import annotations import json 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 shapely.geometry import mapping from shapely.geometry import shape from shapely.ops import transform as transform_geometry from shapely.validation import make_valid from sqlalchemy import Float, 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", } 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", } PRECLIPPED_MUNICIPALITY_PARTITION_OPERATOR_TOOLS = { "provision_regional_bwk_natura2000.py", } SEMANTIC_SELECTION_METRICS: dict[str, tuple[dict[str, Any], ...]] = { "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": (), } SEMANTIC_COUNT_LABELS = { "buildings": "Gebouwen", "population": "Statistische sectoren", "forest": "Bosvlakken", "water": "Waterobjecten", "roads": "Wegsegmenten", "parcels": "Percelen", "nature_value": "BWK-kaartvlakken", "agriculture": "Landbouwgebruikspercelen", } # 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: @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 = { "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", } 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 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) -> 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) if geometry.has_z: geometry = transform_geometry(lambda x, y, z=None: (x, y), geometry) 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 _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 _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, ) -> 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, ) query = db.query(VectorFeature).filter(VectorFeature.dataset_id == dataset_id) 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) elif not full_dataset_area: query = query.filter(ST_Intersects(VectorFeature.geometry, selection_shape)) if 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() ) 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, 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], total_feature_count: int | None = None, selection_geometry: Any | None = None, full_dataset_area: bool = False, preclipped_partition_filter: tuple[str, str] | None = None, ) -> 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 == 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, ) elif not full_dataset_area: selection_filter += (ST_Intersects(VectorFeature.geometry, selection_shape),) selection_is_preclipped = full_dataset_area or preclipped_partition_filter is not None 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) 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, ) for metric_config in metric_configs ] 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, "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, ) -> 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": measured_geometry = ( VectorFeature.geometry if full_dataset_area else func.ST_Intersection(VectorFeature.geometry, selection_shape) ) area_expression = func.ST_Area(func.ST_Transform(measured_geometry, 31370)) 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 {"sum", "mean", "area_weighted_sum"}: 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 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) ) coverage_ratio = intersection_area / func.nullif(source_area, 0.0) value_expression = numeric_value * coverage_ratio aggregate_function = func.avg if method == "mean" else 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: partial_feature_count = ( db.query(func.count(VectorFeature.id)) .filter(*metric_filter) .filter(coverage_ratio < 0.999999) .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, ) -> 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) 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, ) -> 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) 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