Constrain selections to active work areas

backend/app/services/vector_feature_service.py

@@ -8,8 +8,7 @@ from uuid import UUID
 from geoalchemy2.functions import ST_Intersects, ST_MakeEnvelope
 from geoalchemy2.shape import from_shape
 from geoalchemy2.shape import to_shape
-from shapely.geometry import mapping
-from shapely.geometry import shape
+from shapely.geometry import box, mapping, shape
 from shapely.ops import transform as transform_geometry
 from shapely.validation import make_valid
 from sqlalchemy import Float, cast, func
@@ -189,6 +188,27 @@ class VectorFeatureService:
         source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
         return isinstance(source_metadata.get("selection_aggregation"), dict) or VectorFeatureService._dataset_theme(dataset) is not None
 
+    @staticmethod
+    def constrain_bbox_to_area(
+        bbox: dict[str, Any],
+        area_geometry: Any,
+    ) -> tuple[Any, bool]:
+        bbox_geometry = box(
+            float(bbox["min_x"]),
+            float(bbox["min_y"]),
+            float(bbox["max_x"]),
+            float(bbox["max_y"]),
+        )
+        area_shape = to_shape(area_geometry)
+        constrained_geometry = bbox_geometry.intersection(area_shape)
+        if constrained_geometry.is_empty or constrained_geometry.area <= 0:
+            raise AppError(
+                code="VECTOR_SELECTION_OUTSIDE_AREA",
+                message="Selection does not overlap the selected work area",
+                status_code=422,
+            )
+        return from_shape(constrained_geometry, srid=4326), constrained_geometry.equals(area_shape)
+
     @staticmethod
     def can_use_full_area_fast_path(dataset: Dataset, selection_area_id: UUID | None) -> bool:
         if selection_area_id is None or dataset.area_id != selection_area_id:
@@ -343,7 +363,7 @@ class VectorFeatureService:
         if preclipped_partition_filter is not None:
             partition_property, partition_value = preclipped_partition_filter
             query = query.filter(VectorFeature.properties_json.op("->>")(partition_property) == partition_value)
-        elif not full_dataset_area:
+        if not full_dataset_area:
             query = query.filter(ST_Intersects(VectorFeature.geometry, selection_shape))
         if hasattr(query, "count"):
             total_feature_count = int(query.count())
@@ -413,9 +433,9 @@ class VectorFeatureService:
             selection_filter += (
                 VectorFeature.properties_json.op("->>")(partition_property) == partition_value,
             )
-        elif not full_dataset_area:
+        if not full_dataset_area:
             selection_filter += (ST_Intersects(VectorFeature.geometry, selection_shape),)
-        selection_is_preclipped = full_dataset_area or preclipped_partition_filter is not None
+        selection_is_preclipped = full_dataset_area
         feature_count = total_feature_count
         if feature_count is None:
             feature_count = int(db.query(func.count(VectorFeature.id)).filter(*selection_filter).scalar() or 0)