Merge branch 'issue1293-load-stac-start-end-datetime'

soxofaan · soxofaan · commit 2543fd39fd94 · 2025-07-25T14:36:26.000+02:00
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -18,6 +18,7 @@ without compromising stable operations.
 - Configurable usage of `async_task` ([eu-cdse/openeo-cdse-infra#387](https://github.com/eu-cdse/openeo-cdse-infra/issues/387))
 - Add job option "omit-derived-from-links" to omit "derived_from" links in batch job results metadata ([ESA-WEED-project/eo_processing#175](https://github.com/ESA-WEED-project/eo_processing/issues/175))
 - Better freeIPA configurability for proxy user lookup ([#1261](https://github.com/Open-EO/openeo-geopyspark-driver/issues/1261), eu-cdse/openeo-cdse-infra#626)
+- `load_stac`: temporal filtering of STAC Items: prefer more descriptive `start_datetime`+`end_datetime` over nominal `datetime` if possible ([#1293](https://github.com/Open-EO/openeo-geopyspark-driver/issues/1293))
 
 
 ## 0.66.0
diff --git a/openeogeotrellis/load_stac.py b/openeogeotrellis/load_stac.py
@@ -1,6 +1,7 @@
 import re
 
 import datetime as dt
+import datetime
 import json
 import time
 from functools import partial
@@ -48,6 +49,7 @@
 from openeogeotrellis.constants import EVAL_ENV_KEY
 from openeogeotrellis.geopysparkcubemetadata import GeopysparkCubeMetadata
 from openeogeotrellis.geopysparkdatacube import GeopysparkDataCube
+from openeogeotrellis.util.datetime import to_datetime_utc_unless_none
 from openeogeotrellis.utils import normalize_temporal_extent, get_jvm, to_projected_polygons, map_optional, unzip
 from openeogeotrellis.integrations.stac import ResilientStacIO
 
@@ -92,22 +94,7 @@ def load_stac(
     to_date = (dt.datetime.combine(until_date, dt.time.max, until_date.tzinfo) if from_date == until_date
                else until_date - dt.timedelta(milliseconds=1))
 
-    def intersects_spatiotemporally(itm: pystac.Item) -> bool:
-        def intersects_temporally() -> bool:
-            nominal_date = itm.datetime or dateutil.parser.parse(itm.properties["start_datetime"])
-            return from_date <= nominal_date <= to_date
-
-        def intersects_spatially() -> bool:
-            if not requested_bbox or itm.bbox is None:
-                return True
-
-            requested_bbox_lonlat = requested_bbox.reproject("EPSG:4326")
-            return requested_bbox_lonlat.as_polygon().intersects(
-                Polygon.from_bounds(*itm.bbox)
-            )
-
-        return intersects_temporally() and intersects_spatially()
-
+    spatiotemporal_extent = _SpatioTemporalExtent(bbox=requested_bbox, from_date=from_date, to_date=to_date)
 
     def get_pixel_value_offset(itm: pystac.Item, asst: pystac.Asset) -> float:
         raster_scale = asst.extra_fields.get("raster:scale", itm.properties.get("raster:scale", 1.0))
@@ -200,7 +187,7 @@ def operator_value(criterion: Dict[str, object]) -> (str, object):
                                               "proj:shape": asset.get("proj:shape"),
                                           }))
 
-                if intersects_spatiotemporally(pystac_item) and "data" in asset.get("roles", []):
+                if spatiotemporal_extent.item_intersects(pystac_item) and "data" in asset.get("roles", []):
                     pystac_asset = pystac.Asset(
                         href=asset["href"],
                         extra_fields={
@@ -229,7 +216,7 @@ def operator_value(criterion: Dict[str, object]) -> (str, object):
 
                 item = stac_object
                 band_names = stac_metadata_parser.bands_from_stac_item(item=item).band_names()
-                intersecting_items = [item] if intersects_spatiotemporally(item) else []
+                intersecting_items = [item] if spatiotemporal_extent.item_intersects(item) else []
             elif isinstance(stac_object, pystac.Collection) and _supports_item_search(stac_object):
                 collection = stac_object
                 netcdf_with_time_dimension = contains_netcdf_with_time_dimension(collection)
@@ -364,7 +351,7 @@ def intersects_temporally(interval) -> bool:
                     itm
                     for intersecting_catalog in intersecting_catalogs(root=catalog)
                     for itm in intersecting_catalog.get_items()
-                    if intersects_spatiotemporally(itm)
+                    if spatiotemporal_extent.item_intersects(itm)
                 )
 
         jvm = get_jvm()
@@ -723,6 +710,101 @@ def intersects_temporally(interval) -> bool:
     return GeopysparkDataCube(pyramid=gps.Pyramid(levels), metadata=metadata)
 
 
+class _TemporalExtent:
+    """
+    Helper to represent a load_collection/load_stac-style temporal extent
+    with a from_date (inclusive) and to_date (exclusive)
+    and calculate intersection with STAC entities
+    based on nominal datetime or start_datetime+end_datetime
+
+    refs:
+    - https://github.com/radiantearth/stac-spec/blob/master/item-spec/item-spec.md#datetime
+    - https://github.com/radiantearth/stac-spec/blob/master/commons/common-metadata.md#date-and-time-range
+    """
+    # TODO: move this to a more generic location for better reuse
+
+    __slots__ = ("from_date", "to_date")
+
+    def __init__(
+        self,
+        from_date: Union[str, datetime.datetime, datetime.date, None],
+        to_date: Union[str, datetime.datetime, datetime.date, None],
+    ):
+        self.from_date: Union[datetime.datetime, None] = to_datetime_utc_unless_none(from_date)
+        self.to_date: Union[datetime.datetime, None] = to_datetime_utc_unless_none(to_date)
+
+    def intersects(
+        self,
+        nominal: Union[str, datetime.datetime, datetime.date, None] = None,
+        start_datetime: Union[str, datetime.datetime, datetime.date, None] = None,
+        end_datetime: Union[str, datetime.datetime, datetime.date, None] = None,
+    ) -> bool:
+        """
+        Check if the given datetime/interval intersects with the spatiotemporal extent.
+
+        :param nominal: nominal datetime (e.g. typically the "datetime" property of a STAC Item)
+        :param start_datetime: start of the interval (e.g. "start_datetime" property of a STAC Item)
+        :param end_datetime: end of the interval (e.g. "end_datetime" property of a STAC Item)
+        """
+        start_datetime = to_datetime_utc_unless_none(start_datetime)
+        end_datetime = to_datetime_utc_unless_none(end_datetime)
+        nominal = to_datetime_utc_unless_none(nominal)
+
+        # If available, start+end are preferred (cleanly defined interval)
+        # fall back on nominal otherwise
+        if start_datetime and end_datetime and start_datetime <= end_datetime:
+            pass
+        elif nominal:
+            start_datetime = end_datetime = nominal
+        else:
+            raise ValueError(f"Ill-defined instant/interval {nominal=} {start_datetime=} {end_datetime=}")
+
+        return (self.from_date is None or self.from_date <= end_datetime) and (
+            self.to_date is None or start_datetime < self.to_date
+        )
+
+
+class _SpatialExtent:
+    """
+    Helper to represent a spatial extent with a bounding box
+    and calculate intersection with STAC entities (e.g. bbox of a STAC Item).
+    """
+    # TODO: move this to a more generic location for better reuse
+
+    __slots__ = ("bbox", "_bbox_lonlat_shape")
+
+    def __init__(self, *, bbox: Union[BoundingBox, None]):
+        # TODO: support more bbox representations as input
+        self.bbox = bbox
+        self._bbox_lonlat_shape = self.bbox.reproject("EPSG:4326").as_polygon() if self.bbox else None
+
+    def intersects(self, bbox: Union[List[float], Tuple[float, float, float, float], None]):
+        # TODO: this assumes bbox is in lon/lat coordinates, also support other CRSes?
+        if not self.bbox or bbox is None:
+            return True
+        return self._bbox_lonlat_shape.intersects(Polygon.from_bounds(*bbox))
+
+
+class _SpatioTemporalExtent:
+    # TODO: move this to a more generic location for better reuse
+    def __init__(
+        self,
+        *,
+        bbox: Union[BoundingBox, None],
+        from_date: Union[str, datetime.datetime, datetime.date, None],
+        to_date: Union[str, datetime.datetime, datetime.date, None],
+    ):
+        self._spatial_extent = _SpatialExtent(bbox=bbox)
+        self._temporal_extent = _TemporalExtent(from_date=from_date, to_date=to_date)
+
+    def item_intersects(self, item: pystac.Item) -> bool:
+        return self._temporal_extent.intersects(
+            nominal=item.datetime,
+            start_datetime=item.properties.get("start_datetime"),
+            end_datetime=item.properties.get("end_datetime"),
+        ) and self._spatial_extent.intersects(item.bbox)
+
+
 def _is_supported_raster_mime_type(mime_type: str) -> bool:
     mime_type = mime_type.lower()
     # https://github.com/radiantearth/stac-spec/blob/master/best-practices.md#common-media-types-in-stac
diff --git a/openeogeotrellis/util/datetime.py b/openeogeotrellis/util/datetime.py
@@ -0,0 +1,28 @@
+from typing import Union
+import datetime
+
+import dateutil.parser
+
+
+def to_datetime_utc(d: Union[str, datetime.datetime, datetime.date]) -> datetime.datetime:
+    """Parse/convert to datetime in UTC."""
+    if isinstance(d, str):
+        d = dateutil.parser.parse(d)
+    elif isinstance(d, datetime.datetime):
+        pass
+    elif isinstance(d, datetime.date):
+        d = datetime.datetime.combine(d, datetime.time.min)
+    else:
+        raise ValueError(f"Expected str/datetime, but got {type(d)}")
+    if d.tzinfo is None:
+        d = d.replace(tzinfo=datetime.timezone.utc)
+    else:
+        d = d.astimezone(datetime.timezone.utc)
+    return d
+
+
+def to_datetime_utc_unless_none(
+    d: Union[str, datetime.datetime, datetime.date, None]
+) -> Union[datetime.datetime, None]:
+    """Parse/convert to datetime in UTC, but preserve None."""
+    return None if d is None else to_datetime_utc(d)
diff --git a/tests/test_load_stac.py b/tests/test_load_stac.py
@@ -1,6 +1,7 @@
 import dirty_equals
 import pystac
 from contextlib import nullcontext
+import datetime
 
 import mock
 import pytest
@@ -12,6 +13,7 @@
 from openeo_driver.backend import BatchJobMetadata, BatchJobs, LoadParameters
 from openeo_driver.errors import OpenEOApiException
 from openeo_driver.util.date_math import now_utc
+from openeo_driver.util.geometry import BoundingBox
 from openeo_driver.utils import EvalEnv
 
 from openeogeotrellis.load_stac import (
@@ -22,6 +24,9 @@
     _is_band_asset,
     _supports_item_search,
     _get_proj_metadata,
+    _TemporalExtent,
+    _SpatioTemporalExtent,
+    _SpatialExtent,
 )
 
 
@@ -718,3 +723,166 @@ def test_get_proj_metadata_from_asset():
     )
     item = pystac.Item.from_dict(StacDummyBuilder.item())
     assert _get_proj_metadata(asset, item=item) == (32631, (12.0, 34.0, 56.0, 78.0), (12, 34))
+
+
+class TestTemporalExtent:
+    def test_empty(self):
+        extent = _TemporalExtent(None, None)
+        assert extent.intersects("1789-07-14") == True
+        assert extent.intersects(nominal="1789-07-14") == True
+        assert extent.intersects(start_datetime="1914-07-28", end_datetime="1918-11-11") == True
+        assert extent.intersects(nominal="2025-07-24") == True
+
+    def test_nominal_basic(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11", "2025-05-06T22:22:22")
+        assert extent.intersects(nominal="2022-10-11") == False
+        assert extent.intersects(nominal="2025-03-03T12:13:14") == False
+        assert extent.intersects(nominal="2025-03-05T05:05:05") == True
+        assert extent.intersects(nominal="2025-07-07T07:07:07") == False
+
+        assert extent.intersects(nominal=datetime.date(2025, 4, 10)) == True
+        assert extent.intersects(nominal=datetime.datetime(2025, 4, 10, 12)) == True
+
+    def test_nominal_edges(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11", "2025-05-06T22:22:22")
+        assert extent.intersects(nominal="2025-03-04T11:11:10") == False
+        assert extent.intersects(nominal="2025-03-04T11:11:11") == True
+        assert extent.intersects(nominal="2025-03-05T05:05:05") == True
+        assert extent.intersects(nominal="2025-05-06T22:22:21") == True
+        assert extent.intersects(nominal="2025-05-06T22:22:22") == False
+
+    def test_nominal_timezones(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11Z", "2025-05-06T22:22:22-03")
+        assert extent.intersects(nominal="2025-03-04T11:11:10") == False
+        assert extent.intersects(nominal="2025-03-04T11:11:10-02") == True
+        assert extent.intersects(nominal="2025-03-04T11:11:10-04:00") == True
+        assert extent.intersects(nominal="2025-03-04T13:11:11") == True
+        assert extent.intersects(nominal="2025-03-04T13:11:11+02") == True
+        assert extent.intersects(nominal="2025-03-04T13:11:10+02") == False
+
+        assert extent.intersects(nominal="2025-05-06T22:22:22") == True
+        assert extent.intersects(nominal="2025-05-06T22:22:22-02") == True
+        assert extent.intersects(nominal="2025-05-06T22:22:22-03") == False
+        assert extent.intersects(nominal="2025-05-07T01:22:21Z") == True
+        assert extent.intersects(nominal="2025-05-07T01:22:22Z") == False
+
+    def test_nominal_half_open(self):
+        extent = _TemporalExtent(None, "2025-05-06")
+        assert extent.intersects(nominal="1789-07-14") == True
+        assert extent.intersects(nominal="2025-05-05") == True
+        assert extent.intersects(nominal="2025-05-06") == False
+        assert extent.intersects(nominal="2025-11-11") == False
+
+        extent = _TemporalExtent("2025-05-06", None)
+        assert extent.intersects(nominal="2025-05-05") == False
+        assert extent.intersects(nominal="2025-05-06") == True
+        assert extent.intersects(nominal="2099-11-11") == True
+
+    def test_start_end_basic(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11", "2025-05-06T22:22:22")
+        assert extent.intersects(start_datetime="2022-02-02", end_datetime="2022-02-03") == False
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-04-04") == True
+        assert extent.intersects(start_datetime="2025-03-10", end_datetime="2025-04-04") == True
+        assert extent.intersects(start_datetime="2025-03-10", end_datetime="2025-08-08") == True
+        assert extent.intersects(start_datetime="2025-06-10", end_datetime="2025-08-08") == False
+
+    def test_start_end_edges(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11", "2025-05-06T22:22:22")
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T11:11:10") == False
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T11:11:11") == True
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T11:11:12") == True
+
+        assert extent.intersects(start_datetime="2025-05-06T22:22:21", end_datetime="2025-08-08") == True
+        assert extent.intersects(start_datetime="2025-05-06T22:22:22", end_datetime="2025-08-08") == False
+        assert extent.intersects(start_datetime="2025-05-06T22:22:23", end_datetime="2025-08-08") == False
+
+    def test_start_end_timezones(self):
+        extent = _TemporalExtent("2025-03-04T11:11:11Z", "2025-05-06T22:22:22-03")
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T12:12:12") == True
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T12:12:12Z") == True
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T12:12:12+06") == False
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T10:10:10") == False
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-03-04T10:10:10-03") == True
+
+    def test_start_end_half_open(self):
+        extent = _TemporalExtent(None, "2025-05-06")
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-05-05") == True
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-08-08") == True
+        assert extent.intersects(start_datetime="2025-06-06", end_datetime="2025-08-08") == False
+
+        extent = _TemporalExtent("2025-05-06", None)
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-05-05") == False
+        assert extent.intersects(start_datetime="2025-02-02", end_datetime="2025-08-08") == True
+        assert extent.intersects(start_datetime="2025-06-06", end_datetime="2025-08-08") == True
+
+    def test_nominal_vs_start_end(self):
+        """https://github.com/Open-EO/openeo-geopyspark-driver/issues/1293"""
+        extent = _TemporalExtent("2024-02-01", "2024-02-10")
+        assert extent.intersects(nominal="2024-01-01", start_datetime="2024-01-01", end_datetime="2024-12-31") == True
+
+
+class TestSpatialExtent:
+    def test_empty(self):
+        extent = _SpatialExtent(bbox=None)
+        assert extent.intersects(None) is True
+        assert extent.intersects((1, 2, 3, 4)) == True
+
+    def test_basic(self):
+        extent = _SpatialExtent(bbox=BoundingBox(west=3, south=51, east=4, north=52, crs=4326))
+        assert extent.intersects((1, 2, 3, 4)) == False
+        assert extent.intersects((2, 50, 3.1, 51.1)) == True
+        assert extent.intersects((3.3, 51.1, 3.5, 51.5)) == True
+        assert extent.intersects((3.9, 51.9, 4.4, 52.2)) == True
+        assert extent.intersects((5, 51.1, 6, 52.2)) == False
+
+
+class TestSpatioTemporalExtent:
+    @pytest.mark.parametrize(
+        ["bbox", "properties", "expected"],
+        [
+            (
+                [20, 34, 26, 40],
+                {
+                    "datetime": "2024-01-01T00:00:00Z",
+                    "start_datetime": "2024-01-01T00:00:00Z",
+                    "end_datetime": "2024-12-31T23:59:59Z",
+                },
+                True,
+            ),
+            (
+                [20, 34, 26, 40],
+                {"datetime": "2024-01-01T00:00:00Z"},
+                False,
+            ),
+            (
+                [20, 34, 26, 40],
+                {"datetime": "2024-02-02T00:00:00Z"},
+                True,
+            ),
+            (
+                [60, 34, 66, 40],
+                {
+                    "datetime": "2024-01-01T00:00:00Z",
+                    "start_datetime": "2024-01-01T00:00:00Z",
+                    "end_datetime": "2024-12-31T23:59:59Z",
+                },
+                False,
+            ),
+        ],
+    )
+    def test_item_intersects(self, bbox, properties, expected):
+        extent = _SpatioTemporalExtent(
+            bbox=BoundingBox(west=21, south=35, east=25, north=38, crs=4326),
+            from_date="2024-02-01",
+            to_date="2024-02-10",
+        )
+        item = pystac.Item.from_dict(
+            {
+                "type": "Feature",
+                "stac_version": "1.0.0",
+                "id": "2024_GRC_V00",
+                "bbox": bbox,
+                "properties": properties,
+            }
+        )
+        assert extent.item_intersects(item) == expected
diff --git a/tests/util/test_datetime.py b/tests/util/test_datetime.py
