geointel/backend/app/services/change_detection_service.py

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,
            },
        }