924 lines
38 KiB
Python
924 lines
38 KiB
Python
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 box, mapping, 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",
|
|
"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], ...]] = {
|
|
"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.",
|
|
},
|
|
),
|
|
}
|
|
|
|
SEMANTIC_COUNT_LABELS = {
|
|
"buildings": "Gebouwen",
|
|
"population": "Statistische sectoren",
|
|
"forest": "Bosvlakken",
|
|
"water": "Waterobjecten",
|
|
"roads": "Wegsegmenten",
|
|
"parcels": "Percelen",
|
|
"nature_value": "BWK-kaartvlakken",
|
|
"agriculture": "Landbouwgebruikspercelen",
|
|
"soil": "Bodemkaartvlakken",
|
|
}
|
|
|
|
# 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 = {
|
|
"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 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) -> 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 _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,
|
|
) -> 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)
|
|
if 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],
|
|
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,
|
|
) -> 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)
|
|
|
|
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 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
|
|
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 = {
|
|
"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:
|
|
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
|