249 lines
10 KiB
Python
249 lines
10 KiB
Python
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 GeometryCollection
|
|
from shapely.geometry.base import BaseGeometry
|
|
from shapely.ops import unary_union
|
|
from shapely.validation import make_valid
|
|
from shapely.geometry import shape
|
|
from app.core.errors import AppError
|
|
from app.models import Area, Dataset
|
|
from app.schemas.qa import QaProviderComparisonResult
|
|
from app.services.vector_operations_service import VectorOperationsService
|
|
|
|
|
|
def _extract_crs_warnings(source_dataset: Dataset, reference_dataset: Dataset) -> list[str]:
|
|
warnings: list[str] = []
|
|
for dataset, label in ((source_dataset, "candidate"), (reference_dataset, "reference")):
|
|
metadata = dataset.metadata_json
|
|
crs_assumed = None
|
|
if isinstance(metadata, dict):
|
|
crs_assumed = metadata.get("crs_assumed")
|
|
if crs_assumed:
|
|
warnings.append(f"CRS assumption is weak for {label} dataset ({dataset.id}); geometry metrics are approximate")
|
|
if dataset.crs is None:
|
|
warnings.append(f"Missing CRS on {label} dataset ({dataset.id})")
|
|
return warnings
|
|
|
|
|
|
class QaService:
|
|
SUPPORTED_GEOMETRY_TYPES = {"Polygon", "MultiPolygon"}
|
|
|
|
@staticmethod
|
|
def _load_dataset_payload(db, dataset_id: UUID, *, expected_project_id: UUID | None = None) -> tuple[Dataset, dict[str, Any], list[tuple[dict[str, Any], BaseGeometry]]]:
|
|
dataset = db.get(Dataset, dataset_id)
|
|
if not dataset:
|
|
raise AppError(code="DATASET_NOT_FOUND", message="Dataset not found", status_code=404)
|
|
if expected_project_id is not None and dataset.project_id != expected_project_id:
|
|
raise AppError(code="INVALID_DATASET_SCOPE", message="Dataset does not belong to this project", status_code=400)
|
|
if dataset.dataset_type not in {"vector", "geojson"}:
|
|
raise AppError(code="INVALID_DATASET_TYPE", message="Dataset is not a vector dataset", status_code=400)
|
|
|
|
payload, raw_features = VectorOperationsService._load_dataset_payload(dataset)
|
|
geometries = VectorOperationsService._extract_geometries(raw_features)
|
|
return dataset, payload, geometries
|
|
|
|
@staticmethod
|
|
def _apply_area_filter(
|
|
geometries: list[tuple[dict[str, Any], BaseGeometry]],
|
|
area_geometry: BaseGeometry,
|
|
*,
|
|
dataset_id: UUID,
|
|
) -> list[tuple[dict[str, Any], BaseGeometry]]:
|
|
area_geom = area_geometry
|
|
if isinstance(area_geom, GeometryCollection):
|
|
area_geom = unary_union(area_geom.geoms)
|
|
|
|
filtered: list[tuple[dict[str, Any], BaseGeometry]] = []
|
|
for feature, feature_geometry in geometries:
|
|
clipped = feature_geometry.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=f"Area filtering produced invalid geometry for feature in dataset {dataset_id}",
|
|
status_code=400,
|
|
)
|
|
filtered.append((feature, clipped))
|
|
return filtered
|
|
|
|
@staticmethod
|
|
def _validate_area(db, area_id: UUID | None, project_id: UUID, *, dataset_ids: tuple[UUID, UUID]) -> BaseGeometry | None:
|
|
if not area_id:
|
|
return None
|
|
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 != project_id:
|
|
raise AppError(code="INVALID_DATASET_SCOPE", message="Area does not belong to this project", status_code=400)
|
|
if area.id in dataset_ids:
|
|
raise AppError(code="INVALID_PARAMETERS", message="area_id must reference an area, not a dataset", status_code=400)
|
|
|
|
area_geometry = to_shape(area.geometry)
|
|
if area_geometry.is_empty:
|
|
raise AppError(code="INVALID_GEOMETRY", message="Area geometry is empty", status_code=400)
|
|
return area_geometry
|
|
|
|
@staticmethod
|
|
def _match_io_u_metrics(
|
|
source_geometries: list[tuple[dict[str, Any], BaseGeometry]],
|
|
reference_geometries: list[tuple[dict[str, Any], BaseGeometry]],
|
|
iou_threshold: float,
|
|
) -> tuple[int, int, int, list[float], list[str], bool]:
|
|
source_supported = [
|
|
(feature, geom) for feature, geom in source_geometries if geom.geom_type in QaService.SUPPORTED_GEOMETRY_TYPES
|
|
]
|
|
reference_supported = [
|
|
(feature, geom)
|
|
for feature, geom in reference_geometries
|
|
if geom.geom_type in QaService.SUPPORTED_GEOMETRY_TYPES
|
|
]
|
|
|
|
unsupported = sorted(
|
|
{
|
|
geom.geom_type
|
|
for _, geom in source_geometries + reference_geometries
|
|
if geom.geom_type not in QaService.SUPPORTED_GEOMETRY_TYPES
|
|
}
|
|
)
|
|
if not source_supported or not reference_supported:
|
|
return (
|
|
0,
|
|
len(source_supported),
|
|
len(reference_supported),
|
|
[],
|
|
[f"Unsupported geometry types: {unsupported}"] if unsupported else [],
|
|
True,
|
|
)
|
|
|
|
unmatched_reference_indices = set(range(len(reference_supported)))
|
|
matches = 0
|
|
match_iou_values: list[float] = []
|
|
false_positives = 0
|
|
|
|
for _, source_geom in source_supported:
|
|
if source_geom.area <= 0:
|
|
false_positives += 1
|
|
continue
|
|
|
|
best_iou = 0.0
|
|
best_index = None
|
|
for reference_index in list(unmatched_reference_indices):
|
|
_, reference_geom = reference_supported[reference_index]
|
|
if reference_geom.area <= 0:
|
|
unmatched_reference_indices.discard(reference_index)
|
|
continue
|
|
try:
|
|
intersection = source_geom.intersection(reference_geom)
|
|
except Exception as exc: # pragma: no cover - robustness path
|
|
raise AppError(code="GEOMETRY_OPERATION_UNSUPPORTED", message="Geometry operations failed", details={"reason": str(exc)}, status_code=422)
|
|
|
|
if intersection.is_empty:
|
|
continue
|
|
|
|
intersection_area = intersection.area
|
|
if intersection_area < 0:
|
|
intersection_area = 0.0
|
|
union_area = source_geom.area + reference_geom.area - intersection_area
|
|
if union_area <= 0:
|
|
continue
|
|
|
|
candidate_iou = intersection_area / union_area
|
|
if candidate_iou > best_iou:
|
|
best_iou = candidate_iou
|
|
best_index = reference_index
|
|
|
|
if best_index is not None and best_iou >= iou_threshold:
|
|
matches += 1
|
|
match_iou_values.append(best_iou)
|
|
unmatched_reference_indices.discard(best_index)
|
|
else:
|
|
false_positives += 1
|
|
|
|
false_negatives = len(unmatched_reference_indices)
|
|
warnings: list[str] = [f"Unsupported geometry types: {unsupported}"] if unsupported else []
|
|
|
|
return matches, false_positives, false_negatives, match_iou_values, warnings, bool(unsupported)
|
|
|
|
@staticmethod
|
|
def compare_candidate_with_reference(
|
|
db,
|
|
project_id: UUID,
|
|
candidate_dataset_id: UUID,
|
|
reference_dataset_id: UUID,
|
|
iou_threshold: float = 0.5,
|
|
area_id: UUID | None = None,
|
|
) -> QaProviderComparisonResult:
|
|
if candidate_dataset_id == reference_dataset_id:
|
|
raise AppError(code="INVALID_PARAMETERS", message="Candidate and reference dataset must differ", status_code=400)
|
|
|
|
candidate_dataset, candidate_payload, candidate_geometries = QaService._load_dataset_payload(
|
|
db,
|
|
candidate_dataset_id,
|
|
expected_project_id=project_id,
|
|
)
|
|
reference_dataset, reference_payload, reference_geometries = QaService._load_dataset_payload(
|
|
db,
|
|
reference_dataset_id,
|
|
expected_project_id=project_id,
|
|
)
|
|
|
|
area_geometry = QaService._validate_area(
|
|
db,
|
|
area_id=area_id,
|
|
project_id=project_id,
|
|
dataset_ids=(candidate_dataset_id, reference_dataset_id),
|
|
)
|
|
|
|
if area_geometry is not None:
|
|
candidate_geometries = QaService._apply_area_filter(candidate_geometries, area_geometry, dataset_id=candidate_dataset.id)
|
|
reference_geometries = QaService._apply_area_filter(reference_geometries, area_geometry, dataset_id=reference_dataset.id)
|
|
|
|
matches, false_positives, false_negatives, match_iou_values, warnings, unsupported = QaService._match_io_u_metrics(
|
|
candidate_geometries,
|
|
reference_geometries,
|
|
iou_threshold,
|
|
)
|
|
|
|
candidate_feature_count = len(candidate_payload.get("features", [])) if isinstance(candidate_payload, dict) else 0
|
|
reference_feature_count = len(reference_payload.get("features", [])) if isinstance(reference_payload, dict) else 0
|
|
mean_iou = None if not match_iou_values else sum(match_iou_values) / len(match_iou_values)
|
|
|
|
precision = None
|
|
if matches + false_positives > 0:
|
|
precision = matches / (matches + false_positives)
|
|
|
|
recall = None
|
|
if matches + false_negatives > 0:
|
|
recall = matches / (matches + false_negatives)
|
|
|
|
f1_score = None
|
|
if precision is not None and recall is not None and precision + recall > 0:
|
|
f1_score = (2 * precision * recall) / (precision + recall)
|
|
|
|
status = "unsupported" if unsupported else "ok"
|
|
return QaProviderComparisonResult(
|
|
status=status,
|
|
warnings=_extract_crs_warnings(candidate_dataset, reference_dataset) + warnings,
|
|
candidate_feature_count=candidate_feature_count,
|
|
reference_feature_count=reference_feature_count,
|
|
matches=matches,
|
|
false_positives=false_positives,
|
|
false_negatives=false_negatives,
|
|
precision=precision,
|
|
recall=recall,
|
|
f1_score=f1_score,
|
|
mean_iou=mean_iou,
|
|
iou_threshold=iou_threshold,
|
|
unsupported_geometry=unsupported,
|
|
unsupported_geometries=warnings,
|
|
generated_at=datetime.now(timezone.utc),
|
|
)
|