from __future__ import annotations import json import uuid from datetime import datetime, timezone from hashlib import sha256 from pathlib import Path import re from typing import Any from geoalchemy2.shape import to_shape from pyproj import CRS, Transformer from shapely.geometry import GeometryCollection, MultiPolygon, shape from shapely.geometry.base import BaseGeometry from shapely.geometry import mapping from shapely.ops import transform as shapely_transform from shapely.ops import unary_union from shapely.validation import make_valid from sqlalchemy.orm import Session from app.core.errors import AppError from app.models import Area, Dataset, DatasetVersion from app.schemas.dataset import DatasetCreateResponse from app.schemas.operations import VectorOperationResult from app.services.derived_dataset_governance_service import DerivedDatasetGovernanceService from app.services.geojson_service import parse_geojson_payload from app.services.storage_service import StorageService from app.services.vector_feature_service import VectorFeatureService class VectorOperationsService: CANONICAL_VECTOR_CRS = "EPSG:4326" _CHECKSUM_SHA256 = re.compile(r"^[0-9a-f]{64}$", re.IGNORECASE) @staticmethod def _require_vector_dataset(dataset: Dataset) -> None: if dataset.dataset_type not in {"vector", "geojson"}: raise AppError(code="INVALID_DATASET_TYPE", message="Dataset is not a vector dataset", status_code=400) @staticmethod def _load_dataset_payload(dataset: Dataset) -> tuple[dict[str, Any], list[dict[str, Any]]]: if not dataset.storage_path: raise AppError(code="DATASET_FILE_MISSING", message="Stored file missing", status_code=404) path = Path(dataset.storage_path) if not path.exists(): raise AppError(code="DATASET_FILE_MISSING", message="Stored file missing", status_code=404) try: stored_bytes = path.read_bytes() payload = json.loads(stored_bytes.decode("utf-8")) except Exception as exc: raise AppError(code="INVALID_GEOJSON", message="Stored dataset is not valid JSON", status_code=400) from exc if not isinstance(payload, dict) or payload.get("type") != "FeatureCollection": raise AppError(code="INVALID_GEOJSON", message="Dataset payload is not a FeatureCollection", status_code=400) features = payload.get("features") if not isinstance(features, list): raise AppError(code="INVALID_GEOJSON", message="Dataset payload is missing features", status_code=400) # The normal Dataset storage path is a consumption artifact, not a # provenance source archive. Refuse projected/original source bytes # here rather than letting a spatial operation interpret them as # canonical map coordinates. raw_crs = payload.get("crs") if isinstance(raw_crs, dict): crs_properties = raw_crs.get("properties") raw_crs = crs_properties.get("name") if isinstance(crs_properties, dict) else None stored_crs = str(raw_crs or VectorOperationsService.CANONICAL_VECTOR_CRS).strip().upper() dataset_crs = str(dataset.crs or "").strip().upper() if stored_crs != VectorOperationsService.CANONICAL_VECTOR_CRS or ( dataset_crs and dataset_crs != VectorOperationsService.CANONICAL_VECTOR_CRS ): raise AppError( code="DATASET_STORAGE_CRS_MISMATCH", message="Vector operations require canonical EPSG:4326 dataset storage.", details={ "stored_crs": raw_crs or VectorOperationsService.CANONICAL_VECTOR_CRS, "dataset_crs": dataset.crs, "expected_crs": VectorOperationsService.CANONICAL_VECTOR_CRS, }, status_code=409, ) expected_checksum = str(dataset.checksum_sha256 or "").strip().lower() actual_checksum = sha256(stored_bytes).hexdigest() governed_artifact = bool( getattr(dataset, "data_contract_key", None) or ( isinstance(getattr(dataset, "metadata_json", None), dict) and dataset.metadata_json.get("canonical_storage_crs") ) ) if expected_checksum and VectorOperationsService._CHECKSUM_SHA256.fullmatch(expected_checksum): if expected_checksum != actual_checksum: raise AppError( code="DATASET_STORAGE_CHECKSUM_MISMATCH", message="Vector dataset storage no longer matches its validated checksum.", details={"expected_checksum_sha256": expected_checksum, "actual_checksum_sha256": actual_checksum}, status_code=409, ) elif governed_artifact: raise AppError( code="DATASET_STORAGE_CHECKSUM_UNVERIFIABLE", message="Governed vector storage requires a valid SHA-256 checksum before use.", details={"checksum_sha256": dataset.checksum_sha256}, status_code=409, ) return payload, [feature for feature in features if isinstance(feature, dict)] @staticmethod def _extract_geometries(features: list[dict[str, Any]]) -> list[tuple[dict[str, Any], BaseGeometry]]: geometries: list[tuple[dict[str, Any], BaseGeometry]] = [] for feature in features: if not isinstance(feature, dict): continue geometry = feature.get("geometry") if not geometry: continue try: shapely_geom = shape(geometry) except Exception as exc: raise AppError(code="INVALID_GEOMETRY", message="Feature geometry invalid", status_code=400) from exc if not shapely_geom.is_valid: shapely_geom = make_valid(shapely_geom) if not shapely_geom.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Feature geometry cannot be repaired", status_code=400) geometries.append((feature, shapely_geom)) if not geometries: raise AppError(code="EMPTY_VECTOR_DATASET", message="Dataset has no parseable geometries", status_code=422) return geometries @staticmethod def _buffer_in_metres(geometry: BaseGeometry, distance_m: float, source_crs: str) -> BaseGeometry: """Buffer in a Belgian projected CRS, never in angular degrees.""" try: input_crs = CRS.from_user_input(source_crs) metric_crs = CRS.from_epsg(31370) if input_crs == metric_crs: return geometry.buffer(distance_m) forward = Transformer.from_crs(input_crs, metric_crs, always_xy=True) backward = Transformer.from_crs(metric_crs, input_crs, always_xy=True) return shapely_transform(backward.transform, shapely_transform(forward.transform, geometry).buffer(distance_m)) except Exception as exc: raise AppError( code="INVALID_CRS", message="A valid explicit CRS is required for metre-based vector buffering.", status_code=400, ) from exc @staticmethod def inspect(db: Session, dataset_id: uuid.UUID) -> VectorOperationResult: dataset = db.get(Dataset, dataset_id) if not dataset: raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(dataset) payload, features = VectorOperationsService._load_dataset_payload(dataset) geometries = VectorOperationsService._extract_geometries(features) geometry_type_summary: dict[str, int] = {} for _, geometry in geometries: geometry_type_summary[geometry.geom_type] = geometry_type_summary.get(geometry.geom_type, 0) + 1 unioned = unary_union([geometry for _, geometry in geometries]) bounds = unioned.bounds return VectorOperationResult( source_dataset_id=str(dataset_id), feature_count=len(geometries), geometry_type_summary=geometry_type_summary, bounds_json={"min_x": float(bounds[0]), "min_y": float(bounds[1]), "max_x": float(bounds[2]), "max_y": float(bounds[3])}, crs=VectorOperationsService.CANONICAL_VECTOR_CRS, ) @staticmethod def bbox(db: Session, dataset_id: uuid.UUID) -> dict[str, Any]: summary = VectorOperationsService.inspect(db, dataset_id) return { "dataset_id": str(dataset_id), "bounds_json": summary.bounds_json, "feature_count": summary.feature_count, "crs": summary.crs, } @staticmethod def stats(db: Session, dataset_id: uuid.UUID) -> dict[str, Any]: summary = VectorOperationsService.inspect(db, dataset_id) return { "dataset_id": str(dataset_id), "feature_count": summary.feature_count, "geometry_type_summary": summary.geometry_type_summary, "bounds_json": summary.bounds_json, "crs": summary.crs, } @staticmethod def clip_by_area(db: Session, dataset_id: uuid.UUID, area_id: uuid.UUID, output_name: str | None) -> uuid.UUID: source_dataset = db.get(Dataset, dataset_id) if not source_dataset: raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(source_dataset) area = db.get(Area, area_id) if not area: raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404) if area.project_id != source_dataset.project_id: raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to dataset project", status_code=400) payload, features = VectorOperationsService._load_dataset_payload(source_dataset) geometries = VectorOperationsService._extract_geometries(features) area_geom = to_shape(area.geometry) if area_geom.is_empty: raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry is empty", status_code=400) if isinstance(area_geom, GeometryCollection): area_geom = unary_union(area_geom.geoms) if area_geom.geom_type == "MultiPolygon": area_geom = MultiPolygon(area_geom.geoms) if not area_geom.is_valid: area_geom = make_valid(area_geom) if not area_geom.is_valid: raise AppError(code="INVALID_AREA_GEOMETRY", message="Area geometry cannot be repaired", status_code=400) output_features: list[dict[str, Any]] = [] for feature, source_geom in geometries: clipped = source_geom.intersection(area_geom) if clipped.is_empty: continue if not clipped.is_valid: clipped = make_valid(clipped) if not clipped.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Clipped geometry became invalid", status_code=400) output_features.append({ "type": "Feature", "geometry": mapping(clipped), "properties": feature.get("properties", {}) or {}, }) if not output_features: raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Clip operation produced no output features", status_code=422) return VectorOperationsService._persist_derived_dataset( db=db, source_dataset=source_dataset, source_id=dataset_id, operation="clip", feature_collection={"type": "FeatureCollection", "features": output_features}, output_name=output_name, default_name="vector_clipped", ) @staticmethod def buffer(db: Session, dataset_id: uuid.UUID, distance_m: float, dissolve: bool, output_name: str | None) -> uuid.UUID: source_dataset = db.get(Dataset, dataset_id) if not source_dataset: raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(source_dataset) if distance_m <= 0: raise AppError(code="INVALID_PARAMETERS", message="distance_m must be greater than 0", status_code=400) payload, features = VectorOperationsService._load_dataset_payload(source_dataset) source_crs = VectorOperationsService.CANONICAL_VECTOR_CRS geometries = VectorOperationsService._extract_geometries(features) buffered_features = [ (feature, VectorOperationsService._buffer_in_metres(geometry, distance_m, source_crs)) for feature, geometry in geometries ] output_features: list[dict[str, Any]] = [] for feature, geometry in buffered_features: if geometry.is_empty: continue if not geometry.is_valid: geometry = make_valid(geometry) if not geometry.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Buffer geometry became invalid", status_code=400) output_features.append({ "type": "Feature", "geometry": mapping(geometry), "properties": feature.get("properties", {}) or {}, }) if dissolve: dissolved = unary_union([shape(feature["geometry"]) for feature in output_features]) output_features = [{ "type": "Feature", "geometry": mapping(dissolved), "properties": {"operation": "vector_buffer", "distance_m": distance_m, "dissolve": True}, }] if not output_features: raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Buffer operation produced no output features", status_code=422) return VectorOperationsService._persist_derived_dataset( db=db, source_dataset=source_dataset, source_id=dataset_id, operation="buffer", feature_collection={"type": "FeatureCollection", "features": output_features}, output_name=output_name, default_name="vector_buffered", ) @staticmethod def intersect( db: Session, source_dataset_id: uuid.UUID, target_dataset_id: uuid.UUID, output_name: str | None, ) -> uuid.UUID: if source_dataset_id == target_dataset_id: raise AppError(code="INVALID_PARAMETERS", message="other_dataset_id must be different from source dataset", status_code=400) source_dataset = db.get(Dataset, source_dataset_id) if not source_dataset: raise AppError(code="DATASET_NOT_FOUND", message="Source dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(source_dataset) target_dataset = db.get(Dataset, target_dataset_id) if not target_dataset: raise AppError(code="DATASET_NOT_FOUND", message="Target dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(target_dataset) if target_dataset.project_id != source_dataset.project_id: raise AppError(code="INVALID_DATASET_SCOPE", message="Datasets must belong to same project", status_code=400) source_payload, source_features = VectorOperationsService._load_dataset_payload(source_dataset) target_payload, _ = VectorOperationsService._load_dataset_payload(target_dataset) source_geometries = VectorOperationsService._extract_geometries(source_features) target_geometries = VectorOperationsService._extract_geometries(target_payload.get("features", [])) target_union = unary_union([geometry for _, geometry in target_geometries]) output_features: list[dict[str, Any]] = [] for source_feature, source_geometry in source_geometries: intersection = source_geometry.intersection(target_union) if intersection.is_empty: continue if not intersection.is_valid: intersection = make_valid(intersection) if not intersection.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Intersection geometry became invalid", status_code=400) output_features.append({ "type": "Feature", "geometry": mapping(intersection), "properties": source_feature.get("properties", {}) or {}, }) if not output_features: raise AppError(code="VECTOR_OPERATION_EMPTY_RESULT", message="Intersection operation produced no output features", status_code=422) return VectorOperationsService._persist_derived_dataset( db=db, source_dataset=source_dataset, source_id=source_dataset_id, operation="intersect", feature_collection={"type": "FeatureCollection", "features": output_features}, output_name=output_name, default_name="vector_intersect", ) @staticmethod def derive_selection_dataset( db: Session, dataset_id: uuid.UUID, bbox: dict[str, Any], selection_geometry: Any | None = None, selection_area_id: uuid.UUID | None = None, limit: int = 250, output_name: str | None = None, ) -> DatasetCreateResponse: source_dataset = db.get(Dataset, dataset_id) if not source_dataset: raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404) VectorOperationsService._require_vector_dataset(source_dataset) selection = VectorFeatureService.select_features_by_bbox( db, dataset_id=dataset_id, bbox=bbox, selection_geometry=selection_geometry, selection_area_id=selection_area_id, limit=limit, ) if selection["feature_count"] <= 0: raise AppError( code="VECTOR_OPERATION_EMPTY_RESULT", message="Selection produced no output features", status_code=422, ) feature_collection = VectorOperationsService._selection_geojson_for_derived_dataset( selection["geojson"], source_dataset_id=dataset_id, ) derived_id = VectorOperationsService._persist_derived_dataset( db=db, source_dataset=source_dataset, source_id=dataset_id, operation="selection", feature_collection=feature_collection, output_name=output_name, default_name="map_selection", dataset_role="derived", source_name="map_selection", source_metadata={ "selection_bbox": selection["selection_bbox"], "selection_area_id": selection.get("selection_area_id"), "feature_count": selection["feature_count"], "limit": selection["limit"], "truncated": selection["truncated"], "source_table": "vector_features", }, provenance_metadata={ "operation": "map_bbox_selection", "source_dataset_id": str(dataset_id), "source_table": "vector_features", "selection_bbox": selection["selection_bbox"], "selection_area_id": selection.get("selection_area_id"), }, metadata_extra={ "selection_bbox": selection["selection_bbox"], "selection_area_id": selection.get("selection_area_id"), "source_feature_count": selection["feature_count"], "selection_limit": selection["limit"], "selection_truncated": selection["truncated"], "source_dataset_id": str(dataset_id), "source_table": "vector_features", }, persist_vector_features=True, ) derived = db.get(Dataset, derived_id) if not derived: raise AppError(code="DATASET_NOT_FOUND", message="Derived dataset was not persisted", status_code=500) metadata = derived.metadata_json or {} return DatasetCreateResponse( id=derived.id, name=derived.name, dataset_type=derived.dataset_type, source=derived.source, dataset_role=derived.dataset_role, source_name=derived.source_name, reference_layer_name=derived.reference_layer_name, source_metadata=derived.source_metadata, provenance_metadata=derived.provenance_metadata, imported_at=derived.imported_at, project_id=derived.project_id, area_id=derived.area_id, storage_path=derived.storage_path, original_filename=derived.original_filename, stored_filename=derived.stored_filename, content_type=derived.content_type, size_bytes=derived.size_bytes, checksum_sha256=derived.checksum_sha256, crs=derived.crs, bounds_json=derived.bounds_json, resolution_json=derived.resolution_json, bands_json=derived.bands_json, metadata_json=derived.metadata_json, vector_summary=None, status=derived.status, derived_from_dataset_id=derived.derived_from_dataset_id, created_at=derived.created_at, feature_count=metadata.get("feature_count") if isinstance(metadata, dict) else None, ) @staticmethod def _selection_geojson_for_derived_dataset(payload: dict[str, Any], source_dataset_id: uuid.UUID) -> dict[str, Any]: features = payload.get("features") if payload.get("type") != "FeatureCollection" or not isinstance(features, list): raise AppError(code="INVALID_GEOJSON", message="Selection payload must be a FeatureCollection", status_code=500) output_features: list[dict[str, Any]] = [] for feature in features: if not isinstance(feature, dict): continue properties = dict(feature.get("properties") or {}) source_vector_feature_id = properties.pop("vector_feature_id", feature.get("id")) properties.pop("dataset_id", None) properties["source_dataset_id"] = str(source_dataset_id) if source_vector_feature_id is not None: properties["source_vector_feature_id"] = str(source_vector_feature_id) output_features.append( { "type": "Feature", "geometry": feature.get("geometry"), "properties": properties, } ) return {"type": "FeatureCollection", "features": output_features} @staticmethod def _persist_derived_dataset( db: Session, source_dataset: Dataset, source_id: uuid.UUID, operation: str, feature_collection: dict[str, Any], output_name: str | None, default_name: str, dataset_role: str = "derived", source_name: str | None = None, source_metadata: dict[str, Any] | None = None, provenance_metadata: dict[str, Any] | None = None, metadata_extra: dict[str, Any] | None = None, persist_vector_features: bool = False, ) -> uuid.UUID: derived_id = uuid.uuid4() output_name_value = f"{(output_name or default_name)}.geojson" if not output_name_value.strip(): output_name_value = f"{default_name}.geojson" # All current governed vector storage is EPSG:4326. A legacy source # with another CRS is not relabelled here: the derived contract will # quarantine the result instead of placing non-WGS84 coordinates on # the map as if they were WGS84. output_crs = VectorOperationsService.CANONICAL_VECTOR_CRS output_feature_collection = dict(feature_collection) output_feature_collection["crs"] = output_crs stored = json.dumps(output_feature_collection, ensure_ascii=False, separators=(",", ":")).encode("utf-8") storage_info = StorageService.persist_dataset_file( project_id=str(source_dataset.project_id), dataset_id=str(derived_id), dataset_type="vector", original_filename=output_name_value, content=stored, content_type="application/geo+json", ) metadata = parse_geojson_payload(json.dumps(output_feature_collection, ensure_ascii=False, separators=(",", ":"))) if metadata_extra: metadata.update(metadata_extra) derived_dataset = Dataset( id=derived_id, project_id=source_dataset.project_id, area_id=source_dataset.area_id, name=output_name_value, dataset_type="vector", source=f"operation:{operation}", dataset_role=dataset_role, source_name=source_name or "derived", source_metadata=source_metadata, provenance_metadata=provenance_metadata, imported_at=datetime.now(timezone.utc), temporal_series_key=( f"{source_dataset.temporal_series_key}:{operation}" if source_dataset.temporal_series_key else None ), observed_at=source_dataset.observed_at, valid_from=source_dataset.valid_from, valid_to=source_dataset.valid_to, temporal_granularity=source_dataset.temporal_granularity, source_version=source_dataset.source_version, storage_path=storage_info["storage_path"], original_filename=storage_info["original_filename"], stored_filename=storage_info["stored_filename"], content_type=storage_info["content_type"], size_bytes=storage_info["size_bytes"], checksum_sha256=storage_info["checksum_sha256"], derived_from_dataset_id=source_id, crs=metadata.get("crs"), bounds_json=metadata.get("bounds_json"), resolution_json=metadata.get("resolution_json"), bands_json=metadata.get("bands_json"), metadata_json=metadata, status="validating", ) db.add(derived_dataset) dataset_version = DatasetVersion( dataset_id=derived_dataset.id, version=1, storage_path=derived_dataset.storage_path, source_version=derived_dataset.source_version, observed_at=derived_dataset.observed_at, valid_from=derived_dataset.valid_from, valid_to=derived_dataset.valid_to, checksum_sha256=derived_dataset.checksum_sha256, source_metadata=derived_dataset.source_metadata, provenance_metadata=derived_dataset.provenance_metadata, ) db.add(dataset_version) derived_source_key = "map_selection" if source_name == "map_selection" else "derived" is_ready = DerivedDatasetGovernanceService.govern_vector( db, dataset=derived_dataset, dataset_version=dataset_version, feature_collection=output_feature_collection, source_key=derived_source_key, operation=f"vector.{operation}", parent_dataset=source_dataset, operation_parameters={ "operation": operation, "output_name": output_name_value, **(metadata_extra or {}), }, ) if persist_vector_features and is_ready: VectorFeatureService.persist_geojson_features( db=db, dataset_id=derived_dataset.id, payload=output_feature_collection, commit=False, ) db.commit() db.refresh(derived_dataset) return derived_id