from __future__ import annotations from datetime import datetime, timezone from typing import Any from uuid import UUID from geoalchemy2.shape import to_shape from shapely.geometry import mapping from shapely.geometry.base import BaseGeometry from shapely.strtree import STRtree from shapely.validation import make_valid from sqlalchemy.orm import Session from app.core.errors import AppError from app.models import Dataset, VectorFeature from app.schemas.analysis import ChangeDetectionSummary from app.services.vector_operations_service import VectorOperationsService class ChangeDetectionService: SUPPORTED_GEOMETRY_TYPES = {"Polygon", "MultiPolygon"} @staticmethod def compare_vector_datasets( db: Session, *, project_id: UUID, source_dataset_id: UUID, target_dataset_id: UUID, iou_threshold: float = 0.8, include_unchanged: bool = True, modified_threshold: float = 0.3, ) -> ChangeDetectionSummary: if source_dataset_id == target_dataset_id: raise AppError(code="INVALID_PARAMETERS", message="Source and target datasets must differ", status_code=400) if iou_threshold < 0 or iou_threshold > 1: raise AppError(code="INVALID_PARAMETERS", message="iou_threshold must be between 0 and 1", status_code=400) if modified_threshold < 0 or modified_threshold > iou_threshold: raise AppError( code="INVALID_PARAMETERS", message="modified_threshold must be between 0 and iou_threshold", status_code=400, ) source_dataset = ChangeDetectionService._get_project_vector_dataset(db, source_dataset_id, project_id, "Source") target_dataset = ChangeDetectionService._get_project_vector_dataset(db, target_dataset_id, project_id, "Target") source_features, source_warnings = ChangeDetectionService._load_features(db, source_dataset) target_features, target_warnings = ChangeDetectionService._load_features(db, target_dataset) if not source_features: raise AppError(code="EMPTY_VECTOR_DATASET", message="Source dataset has no comparable vector features", status_code=422) if not target_features: raise AppError(code="EMPTY_VECTOR_DATASET", message="Target dataset has no comparable vector features", status_code=422) classified = ChangeDetectionService._classify_features( source_features, target_features, iou_threshold=iou_threshold, modified_threshold=modified_threshold, ) buckets: dict[str, list[dict[str, Any]]] = {"added": [], "removed": [], "modified": [], "unchanged": []} for item in classified: buckets[item["change_type"]].append( ChangeDetectionService._feature( geometry=item["geometry"], change_type=item["change_type"], source_dataset_id=source_dataset_id, target_dataset_id=target_dataset_id, source_feature_id=item["source_feature_id"], target_feature_id=item["target_feature_id"], iou=item["iou"], properties=item["properties"], ) ) unchanged_count = len(buckets["unchanged"]) if not include_unchanged: buckets["unchanged"] = [] geojson_features = buckets["added"] + buckets["removed"] + buckets["modified"] + buckets["unchanged"] return ChangeDetectionSummary( source_dataset_id=source_dataset_id, target_dataset_id=target_dataset_id, source_feature_count=len(source_features), target_feature_count=len(target_features), added_count=len(buckets["added"]), removed_count=len(buckets["removed"]), modified_count=len(buckets["modified"]), unchanged_count=unchanged_count, iou_threshold=iou_threshold, modified_iou_threshold=modified_threshold, warnings=source_warnings + target_warnings, generated_at=datetime.now(timezone.utc), geojson={"type": "FeatureCollection", "features": geojson_features}, ) @staticmethod def _classify_features( source_features: list[dict[str, Any]], target_features: list[dict[str, Any]], *, iou_threshold: float, modified_threshold: float, ) -> list[dict[str, Any]]: """Pair source with target footprints and label how each one changed. Matching is indexed rather than a full cross product: comparing two municipal building layers is otherwise hundreds of millions of geometry intersections. Sources are considered largest first so a big footprint is not left over after a small neighbour claimed its counterpart. """ target_geometries = [feature["geometry"] for feature in target_features] tree = STRtree(target_geometries) if target_geometries else None claimed: set[int] = set() classified: list[dict[str, Any]] = [] order = sorted( range(len(source_features)), key=lambda index: (-source_features[index]["geometry"].area, str(source_features[index]["feature_id"])), ) for source_index in order: source_feature = source_features[source_index] geometry = source_feature["geometry"] best_iou = 0.0 best_index: int | None = None candidates = [] if tree is None else sorted(int(value) for value in tree.query(geometry)) for target_index in candidates: if target_index in claimed: continue candidate_iou = ChangeDetectionService._iou(geometry, target_geometries[target_index]) if candidate_iou > best_iou: best_iou = candidate_iou best_index = target_index if best_index is not None and best_iou >= iou_threshold: claimed.add(best_index) change_type = "unchanged" elif best_index is not None and best_iou >= modified_threshold: # The same object, redrawn: an annexe, a demolition of one wing, # or a resurvey. Reporting it as removed + added would hide it. claimed.add(best_index) change_type = "modified" else: change_type = "removed" classified.append( { "change_type": change_type, "geometry": geometry if change_type != "modified" else target_geometries[best_index], "source_feature_id": source_feature["feature_id"], "target_feature_id": target_features[best_index]["feature_id"] if change_type != "removed" else None, "iou": best_iou if best_iou > 0 else None, "properties": source_feature["properties"], } ) classified.extend( { "change_type": "added", "geometry": target_feature["geometry"], "source_feature_id": None, "target_feature_id": target_feature["feature_id"], "iou": None, "properties": target_feature["properties"], } for target_index, target_feature in enumerate(target_features) if target_index not in claimed ) return classified @staticmethod def _get_project_vector_dataset(db: Session, dataset_id: UUID, project_id: UUID, label: str) -> Dataset: dataset = db.get(Dataset, dataset_id) if not dataset: raise AppError(code="DATASET_NOT_FOUND", message=f"{label} dataset not found", status_code=404) if dataset.project_id != project_id: raise AppError(code="INVALID_DATASET_SCOPE", message=f"{label} dataset does not belong to this project", status_code=400) VectorOperationsService._require_vector_dataset(dataset) return dataset @staticmethod def _load_features(db: Session, dataset: Dataset) -> tuple[list[dict[str, Any]], list[str]]: rows = db.query(VectorFeature).filter(VectorFeature.dataset_id == dataset.id).all() warnings: list[str] = [] if rows: return [ChangeDetectionService._row_to_feature(row) for row in rows], warnings warnings.append(f"Dataset {dataset.id} has no persisted vector_features; falling back to stored GeoJSON artifact") _payload, raw_features = VectorOperationsService._load_dataset_payload(dataset) extracted = VectorOperationsService._extract_geometries(raw_features) return [ ChangeDetectionService._raw_feature_to_feature(index, raw_feature, geometry) for index, (raw_feature, geometry) in enumerate(extracted) ], warnings @staticmethod def _row_to_feature(row: VectorFeature) -> dict[str, Any]: geometry = ChangeDetectionService._valid_comparable_geometry(to_shape(row.geometry)) return { "feature_id": str(row.source_feature_id or row.id), "properties": dict(row.properties_json or {}), "geometry": geometry, } @staticmethod def _raw_feature_to_feature(index: int, raw_feature: dict[str, Any], geometry: BaseGeometry) -> dict[str, Any]: properties = raw_feature.get("properties") if isinstance(raw_feature.get("properties"), dict) else {} source_id = raw_feature.get("id") or properties.get("id") or properties.get("source_feature_id") or str(index) return { "feature_id": str(source_id), "properties": dict(properties), "geometry": ChangeDetectionService._valid_comparable_geometry(geometry), } @staticmethod def _valid_comparable_geometry(geometry: BaseGeometry) -> BaseGeometry: if geometry.is_empty: raise AppError(code="INVALID_GEOMETRY", message="Empty geometry cannot be compared", status_code=400) if not geometry.is_valid: geometry = make_valid(geometry) if geometry.is_empty or not geometry.is_valid: raise AppError(code="INVALID_GEOMETRY", message="Geometry cannot be repaired for comparison", status_code=400) if geometry.geom_type not in ChangeDetectionService.SUPPORTED_GEOMETRY_TYPES: raise AppError( code="UNSUPPORTED_GEOMETRY", message="Change detection supports Polygon and MultiPolygon geometries only", details={"geometry_type": geometry.geom_type}, status_code=422, ) return geometry @staticmethod def _iou(left: BaseGeometry, right: BaseGeometry) -> float: if left.area <= 0 or right.area <= 0: return 0.0 intersection = left.intersection(right) if intersection.is_empty: return 0.0 union_area = left.area + right.area - intersection.area if union_area <= 0: return 0.0 return float(intersection.area / union_area) @staticmethod def _feature( *, geometry: BaseGeometry, change_type: str, source_dataset_id: UUID, target_dataset_id: UUID, source_feature_id: str | None, target_feature_id: str | None, iou: float | None, properties: dict[str, Any], ) -> dict[str, Any]: return { "type": "Feature", "geometry": mapping(geometry), "properties": { **properties, "change_type": change_type, "source_dataset_id": str(source_dataset_id), "target_dataset_id": str(target_dataset_id), "source_feature_id": source_feature_id, "target_feature_id": target_feature_id, "iou": iou, }, }