perf: optimize trusted full-area analysis

CHANGELOG.md

@@ -15,6 +15,7 @@
 - Kept every fetch operator-triggered, idempotent and behind the canonical DatasetService upload path; no startup fetch, migration or API contract change was introduced.
 - Preserved separate Mol/regional series keys and honest partial-sector population and 10 m forest-area limitations.
 - Replaced internal provider identifiers with readable source labels in the primary map.
+- Added a provenance-gated full-Area query path for official pre-clipped datasets, avoiding redundant intersection of every feature against the same detailed municipal/regional boundary while preserving exact rectangle selection behavior.
 - Added focused scope filtering, command construction, boundary resolution, multi-NIS provenance, packaging and frontend-label tests.
 
 ## Sprint 193 End-user regional workbench simplification (2026-07-14)

backend/README.md

@@ -1043,6 +1043,14 @@ stay within upstream response limits, then builds one retained 10 m mosaic and
 one normal regional vector Dataset. A failed source request leaves completed
 partition artifacts reusable and never lowers source resolution silently.
 
+Official operator datasets record that their geometries were clipped to the
+persisted Area. When that exact Area is selected, vector totals and aggregate
+metrics use the already clipped geometries directly rather than intersecting
+every row with the same detailed boundary again. This optimization is allowed
+only for matching Dataset/Area ids with explicit clipping metadata or a known
+clipping operator; drawn rectangles and ordinary uploads keep the normal exact
+PostGIS intersection path.
+
 ## Helpful repository scripts
 
 - `bash scripts/backend_install.sh`

backend/app/api/routes/datasets.py

@@ -248,10 +248,13 @@ def select_vector_features(
         "bbox": payload.bbox.model_dump(),
         "limit": payload.limit,
     }
+    full_dataset_area = False
     if selection_area is not None:
+        full_dataset_area = VectorFeatureService.can_use_full_area_fast_path(dataset, selection_area.id)
         selection_kwargs.update(
             selection_geometry=selection_area.geometry,
             selection_area_id=selection_area.id,
+            full_dataset_area=full_dataset_area,
         )
     result = VectorFeatureService.select_features_by_bbox(db, **selection_kwargs)
     if isinstance(dataset.source_metadata, dict) and dataset.source_metadata.get("selection_aggregation"):
@@ -262,6 +265,7 @@ def select_vector_features(
         }
         if selection_area is not None:
             summary_kwargs["selection_geometry"] = selection_area.geometry
+            summary_kwargs["full_dataset_area"] = full_dataset_area
         result["summary"] = VectorFeatureService.summarize_features_by_bbox(db, **summary_kwargs)
     return envelope(VectorSelectionResponse(**result).model_dump(exclude_none=True))
 

backend/app/services/vector_feature_service.py

@@ -18,7 +18,25 @@ from app.core.errors import AppError
 from app.models import Dataset, VectorFeature
 
 
+FULL_AREA_CLIPPED_OPERATOR_TOOLS = {
+    "provision_mol_population_history.py",
+    "provision_official_landuse_timeseries.py",
+    "provision_regional_grb_buildings.py",
+    "provision_regional_grb_context.py",
+}
+
+
 class VectorFeatureService:
+    @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:
+            return False
+        source_metadata = dataset.source_metadata if isinstance(dataset.source_metadata, dict) else {}
+        if source_metadata.get("geometry_clipped_to_area") is True:
+            return True
+        provenance = dataset.provenance_metadata if isinstance(dataset.provenance_metadata, dict) else {}
+        return provenance.get("operator_tool") in FULL_AREA_CLIPPED_OPERATOR_TOOLS
+
     @staticmethod
     def _feature_row(dataset_id: UUID, feature: dict[str, Any], index: int, feature_class: str | None) -> VectorFeature | None:
         geometry_payload = feature.get("geometry")
@@ -126,6 +144,7 @@ class VectorFeatureService:
         dataset: Dataset | None = None,
         selection_geometry: Any | None = None,
         selection_area_id: UUID | None = None,
+        full_dataset_area: bool = False,
     ) -> dict[str, Any]:
         normalized_bbox = VectorFeatureService._normalize_selection_bbox(bbox)
         safe_limit = max(1, min(int(limit), 1000))
@@ -139,11 +158,9 @@ class VectorFeatureService:
                 4326,
             )
 
-        query = (
-            db.query(VectorFeature)
-            .filter(VectorFeature.dataset_id == dataset_id)
-            .filter(ST_Intersects(VectorFeature.geometry, selection_shape))
-        )
+        query = db.query(VectorFeature).filter(VectorFeature.dataset_id == dataset_id)
+        if not full_dataset_area:
+            query = query.filter(ST_Intersects(VectorFeature.geometry, selection_shape))
         if hasattr(query, "count"):
             total_feature_count = int(query.count())
         else:  # Lightweight unit-test sessions do not always implement Query.count().
@@ -165,6 +182,7 @@ class VectorFeatureService:
                 bbox=normalized_bbox,
                 total_feature_count=total_feature_count,
                 selection_geometry=selection_geometry,
+                full_dataset_area=full_dataset_area,
             )
 
         result = {
@@ -191,6 +209,7 @@ class VectorFeatureService:
         bbox: dict[str, Any],
         total_feature_count: int | None = None,
         selection_geometry: Any | None = None,
+        full_dataset_area: bool = False,
     ) -> dict[str, Any]:
         normalized_bbox = VectorFeatureService._normalize_selection_bbox(bbox)
         selection_shape = selection_geometry
@@ -202,10 +221,9 @@ class VectorFeatureService:
                 normalized_bbox["max_y"],
                 4326,
             )
-        selection_filter = (
-            VectorFeature.dataset_id == dataset.id,
-            ST_Intersects(VectorFeature.geometry, selection_shape),
-        )
+        selection_filter = (VectorFeature.dataset_id == dataset.id,)
+        if not full_dataset_area:
+            selection_filter += (ST_Intersects(VectorFeature.geometry, selection_shape),)
         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)
@@ -222,14 +240,22 @@ class VectorFeatureService:
 
         metric_value = float(feature_count)
         if method == "intersection_area":
-            intersection = func.ST_Intersection(VectorFeature.geometry, selection_shape)
-            area_expression = func.ST_Area(func.ST_Transform(intersection, 31370))
+            measured_geometry = (
+                VectorFeature.geometry
+                if full_dataset_area
+                else func.ST_Intersection(VectorFeature.geometry, selection_shape)
+            )
+            area_expression = func.ST_Area(func.ST_Transform(measured_geometry, 31370))
             area_m2 = db.query(func.coalesce(func.sum(area_expression), 0.0)).filter(*selection_filter).scalar()
             divisor = 10_000.0 if unit == "ha" else 1.0
             metric_value = float(area_m2 or 0.0) / divisor
         elif method == "intersection_length":
-            intersection = func.ST_Intersection(VectorFeature.geometry, selection_shape)
-            length_expression = func.ST_Length(func.ST_Transform(intersection, 31370))
+            measured_geometry = (
+                VectorFeature.geometry
+                if full_dataset_area
+                else func.ST_Intersection(VectorFeature.geometry, selection_shape)
+            )
+            length_expression = func.ST_Length(func.ST_Transform(measured_geometry, 31370))
             length_m = db.query(func.coalesce(func.sum(length_expression), 0.0)).filter(*selection_filter).scalar()
             divisor = 1_000.0 if unit == "km" else 1.0
             metric_value = float(length_m or 0.0) / divisor
@@ -244,7 +270,7 @@ class VectorFeatureService:
                 )
             numeric_value = cast(VectorFeature.properties_json.op("->>")(property_name), Float)
             value_expression = numeric_value
-            if method == "area_weighted_sum":
+            if method == "area_weighted_sum" and not full_dataset_area:
                 source_area = func.ST_Area(func.ST_Transform(VectorFeature.geometry, 31370))
                 intersection_area = func.ST_Area(
                     func.ST_Transform(func.ST_Intersection(VectorFeature.geometry, selection_shape), 31370)
@@ -258,7 +284,7 @@ class VectorFeatureService:
                 .scalar()
             )
             metric_value = float(aggregate_value or 0.0)
-            if method == "area_weighted_sum":
+            if method == "area_weighted_sum" and not full_dataset_area:
                 partial_feature_count = (
                     db.query(func.count(VectorFeature.id))
                     .filter(*selection_filter)
@@ -268,6 +294,10 @@ class VectorFeatureService:
                 is_estimate = bool(partial_feature_count)
                 if not is_estimate and config.get("warning_only_when_estimate", True):
                     warning = None
+            elif method == "area_weighted_sum":
+                is_estimate = False
+                if config.get("warning_only_when_estimate", True):
+                    warning = None
         elif method != "feature_count":
             raise AppError(
                 code="INVALID_SELECTION_AGGREGATION",

backend/tests/test_sprint194_regional_timeseries.py

@@ -7,11 +7,16 @@ import json
 from pathlib import Path
 import sys
 import zipfile
+from types import SimpleNamespace
+from uuid import uuid4
 
 import numpy as np
 import rasterio
 from rasterio.transform import from_origin
 
+from app.models import Dataset
+from app.services.vector_feature_service import VectorFeatureService
+
 
 ROOT = Path(__file__).resolve().parents[2]
 SCRIPTS = ROOT / "scripts"
@@ -184,3 +189,94 @@ def test_end_user_dataset_sources_are_human_readable() -> None:
     assert "statbel: 'Statbel'" in display
     assert "getDatasetSourceDisplayName(activeThemeDataset)" in workspace
     assert "dataset ? getDatasetSourceDisplayName(dataset)" in workspace
+
+
+def test_full_area_fast_path_requires_matching_area_and_clipped_operator_provenance() -> None:
+    project_id = uuid4()
+    area_id = uuid4()
+    trusted = Dataset(
+        id=uuid4(),
+        project_id=project_id,
+        area_id=area_id,
+        name="Regional forest",
+        dataset_type="vector",
+        source="operator_official_import",
+        provenance_metadata={"operator_tool": "provision_official_landuse_timeseries.py"},
+    )
+    explicit = Dataset(
+        id=uuid4(),
+        project_id=project_id,
+        area_id=area_id,
+        name="Clipped vector",
+        dataset_type="vector",
+        source="manual",
+        source_metadata={"geometry_clipped_to_area": True},
+    )
+    untrusted = Dataset(
+        id=uuid4(),
+        project_id=project_id,
+        area_id=area_id,
+        name="Assigned only",
+        dataset_type="vector",
+        source="manual",
+    )
+
+    assert VectorFeatureService.can_use_full_area_fast_path(trusted, area_id) is True
+    assert VectorFeatureService.can_use_full_area_fast_path(explicit, area_id) is True
+    assert VectorFeatureService.can_use_full_area_fast_path(untrusted, area_id) is False
+    assert VectorFeatureService.can_use_full_area_fast_path(trusted, uuid4()) is False
+    assert VectorFeatureService.can_use_full_area_fast_path(trusted, None) is False
+
+
+def test_full_area_summary_uses_exact_stored_values_without_partial_intersection() -> None:
+    class ScalarQuery:
+        def __init__(self, value: float) -> None:
+            self.value = value
+
+        def filter(self, *_args):
+            return self
+
+        def scalar(self):
+            return self.value
+
+    class ScalarSession:
+        def __init__(self, value: float) -> None:
+            self.value = value
+            self.query_count = 0
+
+        def query(self, *_args):
+            self.query_count += 1
+            return ScalarQuery(self.value)
+
+    dataset = Dataset(
+        id=uuid4(),
+        project_id=uuid4(),
+        name="Population",
+        dataset_type="vector",
+        source="operator_official_import",
+        source_metadata={
+            "selection_aggregation": {
+                "method": "area_weighted_sum",
+                "property": "population_total",
+                "label": "Inwoners",
+                "unit": "inwoners",
+                "warning_only_when_estimate": True,
+                "warning": "Partial-sector estimate",
+            }
+        },
+    )
+    session = ScalarSession(506_473.0)
+
+    result = VectorFeatureService.summarize_features_by_bbox(
+        session,
+        dataset=dataset,
+        bbox={"min_x": 4.5, "min_y": 51.0, "max_x": 5.3, "max_y": 51.6, "crs": "EPSG:4326"},
+        total_feature_count=733,
+        selection_geometry=SimpleNamespace(),
+        full_dataset_area=True,
+    )
+
+    assert result["metric_value"] == 506_473.0
+    assert result["is_estimate"] is False
+    assert result["warning"] is None
+    assert session.query_count == 1

docs/CODEX_EXECUTION_LOG.md

@@ -7982,6 +7982,7 @@ Live source finding:
 - The regional Statbel synchronization imported all five requested snapshots successfully.
 - The first complete-region forest request was rejected by MercatorNet before download because the estimated 97.68 MB response exceeded its 78.12 MB service limit.
 - The operator now uses the 28 retained official municipality boundaries as resumable WCS partitions, merges them locally at the unchanged 10 m grid and applies the exact region-union clip before polygonization. It does not lower resolution or truncate the source.
+- Live full-region analysis returned the correct six-theme result but initially took about 92 seconds because every already-clipped feature was intersected again with a 947 kB regional geometry. Added a provenance-gated full-Area path that skips only this redundant intersection; arbitrary rectangles and untrusted uploads retain the existing exact query.
 
 ## Sprint 190 Regional Kempen GRB buildings (2026-07-14)
 

scripts/provision_mol_population_history.py

@@ -255,6 +255,7 @@ def build_snapshot(
         "scope_type": scope.scope_type,
         "member_count": len(scope.members),
         "member_nis_codes": list(scope.nis_codes),
+        "geometry_clipped_to_area": True,
         "observation_year": year,
         "missing_population_sector_count": missing_population,
         "attribution": ATTRIBUTION,
@@ -301,6 +302,7 @@ def upload_snapshot(
         "scope_display_name": scope.display_name,
         "member_count": len(scope.members),
         "member_nis_codes": list(scope.nis_codes),
+        "geometry_clipped_to_area": True,
         "attribution": ATTRIBUTION,
         "license": "CC BY 4.0",
         "temporal_series_label": "Officiële bevolkingscijfers per statistische sector",
@@ -320,6 +322,7 @@ def upload_snapshot(
         "operator_tool": "provision_mol_population_history.py",
         "operator_explicit_fetch": True,
         "scope_key": scope.key,
+        "geometry_clipped_to_area": True,
         "sector_geometry_url": SECTOR_URL.format(year=year),
         "population_url": POPULATION_URLS[year],
         "generated_at": datetime.now(timezone.utc).isoformat(),

scripts/provision_official_landuse_timeseries.py

@@ -788,6 +788,7 @@ def build_source_metadata(args: argparse.Namespace, snapshot: PreparedSnapshot)
         "authority_level": "authoritative",
         "coverage_scope": args.scope_key,
         **identity,
+        "geometry_clipped_to_area": True,
         "attribution": ATTRIBUTION,
         "license_note": "Publieke Vlaamse overheidsdata; raadpleeg de toegangs- en gebruiksvoorwaarden in de bronmetadata.",
         "methodology_version": "3",
@@ -815,6 +816,7 @@ def build_provenance_metadata(args: argparse.Namespace, snapshot: PreparedSnapsh
     return {
         "operator_tool": "provision_official_landuse_timeseries.py",
         "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
         "wcs_url": WCS_URL,
         "wcs_version": WCS_VERSION,
         "coverage_id": coverage_id(snapshot.year),

scripts/provision_regional_grb_buildings.py

@@ -657,6 +657,7 @@ def provision_dataset(
         "scope_limitation": scope.limitation_message,
         "member_count": len(scope.members),
         "member_nis_codes": list(scope.nis_codes),
+        "geometry_clipped_to_area": True,
         "feature_count": manifest["feature_count"],
         "partition_count": len(partition_paths),
         "partition_strategy": manifest["partition_strategy"],
@@ -671,6 +672,7 @@ def provision_dataset(
     provenance_metadata = {
         "operator_tool": "provision_regional_grb_buildings.py",
         "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
         "manifest_path": str(manifest_path),
         "source_url": GRB_GBG_ITEMS_URL,
         "artifact_sha256": manifest["artifact_sha256"],

scripts/provision_regional_grb_context.py

@@ -716,6 +716,7 @@ def provision_dataset(
         "layer_limitation": definition.limitation_message,
         "member_count": len(scope.members),
         "member_nis_codes": list(scope.nis_codes),
+        "geometry_clipped_to_area": True,
         "feature_count": manifest["feature_count"],
         "partition_count": len(partition_paths),
         "partition_strategy": manifest["partition_strategy"],
@@ -730,6 +731,7 @@ def provision_dataset(
     provenance_metadata = {
         "operator_tool": "provision_regional_grb_context.py",
         "operator_explicit_fetch": True,
+        "geometry_clipped_to_area": True,
         "manifest_path": str(manifest_path),
         "source_urls": manifest["grb_source_urls"],
         "artifact_sha256": manifest["artifact_sha256"],