JensandClaude Opus 5
12aaf1bb4d
Change detection was the one analysis that ignored the selection entirely. It compared two datasets in full, loaded every feature of both into Python with no spatial predicate, and — with include_unchanged defaulting to true — returned a FeatureCollection holding both datasets. For a regional building layer that is the wrong answer to "what changed here" and a response no browser should be asked to hold. It now accepts bbox and area_id, resolved the way every other analysis resolves them, and loads through an indexed ST_Intersects predicate. Features are deliberately not clipped to the selection. A change class describes a whole object: comparing a clipped earlier footprint against an unclipped later one would report the selection edge itself as a change. Objects the edge crosses are compared in full and counted in a warning. The returned geometry is capped by preview_limit, spending that budget on modified, added and removed before unchanged, while every count still describes the whole selection. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
413 lines
17 KiB
Python
413 lines
17 KiB
Python
from __future__ import annotations
|
|
|
|
from datetime import datetime, timezone
|
|
from typing import Any
|
|
from uuid import UUID
|
|
|
|
from geoalchemy2.shape import from_shape, 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 import func
|
|
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,
|
|
bbox: dict[str, Any] | None = None,
|
|
area_id: UUID | None = None,
|
|
preview_limit: int = 2_000,
|
|
) -> 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")
|
|
|
|
selection_geometry = ChangeDetectionService._selection_geometry(db, project_id, bbox=bbox, area_id=area_id)
|
|
|
|
source_features, source_warnings = ChangeDetectionService._load_features(db, source_dataset, selection_geometry)
|
|
target_features, target_warnings = ChangeDetectionService._load_features(db, target_dataset, selection_geometry)
|
|
|
|
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)
|
|
|
|
source_features = ChangeDetectionService.restrict_to_selection(source_features, selection_geometry, label="Source")
|
|
target_features = ChangeDetectionService.restrict_to_selection(target_features, selection_geometry, label="Target")
|
|
|
|
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, preview_truncated = ChangeDetectionService.limit_preview(
|
|
buckets["added"] + buckets["removed"] + buckets["modified"] + buckets["unchanged"],
|
|
limit=preview_limit,
|
|
)
|
|
warnings = source_warnings + target_warnings
|
|
edge_count = sum(
|
|
1 for feature in source_features + target_features if feature.get("partially_covered")
|
|
)
|
|
if edge_count:
|
|
warnings.append(
|
|
f"{edge_count} objecten liggen deels buiten de selectie. Ze zijn volledig vergeleken, zodat de "
|
|
"selectierand zelf geen wijziging veroorzaakt."
|
|
)
|
|
if preview_truncated:
|
|
warnings.append(
|
|
f"De tellingen gelden voor de volledige selectie; de kaart toont maximaal {preview_limit} objecten, "
|
|
"wijzigingen eerst."
|
|
)
|
|
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,
|
|
selection_area_id=area_id,
|
|
preview_limit=preview_limit,
|
|
preview_truncated=preview_truncated,
|
|
warnings=warnings,
|
|
generated_at=datetime.now(timezone.utc),
|
|
geojson={"type": "FeatureCollection", "features": geojson_features},
|
|
)
|
|
|
|
@staticmethod
|
|
def _selection_geometry(
|
|
db: Session,
|
|
project_id: UUID,
|
|
*,
|
|
bbox: dict[str, Any] | None,
|
|
area_id: UUID | None,
|
|
) -> BaseGeometry | None:
|
|
"""Resolve the drawn rectangle against the named work area, if any."""
|
|
|
|
from app.models import Area
|
|
from shapely.geometry import box as shapely_box
|
|
|
|
selection = None
|
|
if bbox:
|
|
selection = shapely_box(
|
|
float(bbox["min_x"]), float(bbox["min_y"]), float(bbox["max_x"]), float(bbox["max_y"])
|
|
)
|
|
if area_id is None:
|
|
return selection
|
|
|
|
area = db.get(Area, area_id)
|
|
if area is None or area.project_id != project_id:
|
|
raise AppError(code="AREA_NOT_FOUND", message="Area not found", status_code=404)
|
|
area_geometry = to_shape(area.geometry)
|
|
if selection is None:
|
|
return area_geometry
|
|
intersection = selection.intersection(area_geometry)
|
|
if intersection.is_empty or intersection.area <= 0:
|
|
raise AppError(
|
|
code="CHANGE_DETECTION_SELECTION_OUTSIDE_AREA",
|
|
message="Selection does not overlap the selected work area",
|
|
status_code=422,
|
|
)
|
|
return intersection
|
|
|
|
# Order the preview spends its budget in. An operator asking what changed
|
|
# is not helped by a cap filled with unchanged footprints.
|
|
PREVIEW_PRIORITY = {"modified": 0, "added": 1, "removed": 2, "unchanged": 3}
|
|
|
|
@staticmethod
|
|
def restrict_to_selection(
|
|
features: list[dict[str, Any]],
|
|
selection_geometry: BaseGeometry | None,
|
|
*,
|
|
label: str = "Dataset",
|
|
) -> list[dict[str, Any]]:
|
|
"""Keep the features a drawn selection reaches, and say which it cuts.
|
|
|
|
Geometry is deliberately *not* clipped. A change class describes a whole
|
|
object: comparing a clipped 2020 footprint against an unclipped 2024 one
|
|
would manufacture "modified" along the selection edge. Clipping is right
|
|
for an area metric and wrong for an identity comparison.
|
|
"""
|
|
|
|
if selection_geometry is None:
|
|
return features
|
|
|
|
kept: list[dict[str, Any]] = []
|
|
for feature in features:
|
|
geometry = feature["geometry"]
|
|
if not geometry.intersects(selection_geometry):
|
|
continue
|
|
kept.append({**feature, "partially_covered": not selection_geometry.covers(geometry)})
|
|
|
|
if not kept:
|
|
raise AppError(
|
|
code="CHANGE_DETECTION_SELECTION_EMPTY",
|
|
message=f"{label} dataset has no features inside this selection",
|
|
status_code=422,
|
|
)
|
|
return kept
|
|
|
|
@staticmethod
|
|
def limit_preview(
|
|
features: list[dict[str, Any]],
|
|
*,
|
|
limit: int,
|
|
) -> tuple[list[dict[str, Any]], bool]:
|
|
"""Cap the returned geometry without capping the counts.
|
|
|
|
``include_unchanged`` defaulted to true and nothing bounded the result,
|
|
so a regional comparison returned a FeatureCollection holding both
|
|
datasets in full. The counts describe the whole selection; the preview
|
|
describes what a map can usefully draw.
|
|
"""
|
|
|
|
if limit <= 0 or len(features) <= limit:
|
|
return features, False
|
|
ordered = sorted(
|
|
features,
|
|
key=lambda item: ChangeDetectionService.PREVIEW_PRIORITY.get(item["change_type"], 9),
|
|
)
|
|
return ordered[:limit], True
|
|
|
|
@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,
|
|
selection_geometry: BaseGeometry | None = None,
|
|
) -> tuple[list[dict[str, Any]], list[str]]:
|
|
query = db.query(VectorFeature).filter(VectorFeature.dataset_id == dataset.id)
|
|
if selection_geometry is not None and hasattr(query, "filter"):
|
|
# Bound the load in the database. Pulling a regional building layer
|
|
# into Python to then discard most of it costs memory and time for
|
|
# nothing, and the fallback below has no such option.
|
|
try:
|
|
query = query.filter(
|
|
func.ST_Intersects(VectorFeature.geometry, from_shape(selection_geometry, srid=4326))
|
|
)
|
|
except Exception:
|
|
# Lightweight unit-test sessions do not implement every spatial
|
|
# predicate; restrict_to_selection still bounds the population.
|
|
pass
|
|
rows = query.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,
|
|
},
|
|
}
|