Subsetting Sentinel-3 products based on footprints given in the STAC metadata

CHANGES.md

@@ -4,6 +4,9 @@
   (and its subpackages) is included in the PyPI wheel.
 - Remove the `coarsen.py` module, as it has been moved to [xcube-resampling](https://github.com/xcube-dev/xcube-resampling) 
   and is no longer used internally.
+- Added footprint-based subsetting for Sentinel-3 OLCI and SLSTR LST using STAC 
+  metadata, improving performance by avoiding full latitude/longitude grid downloads 
+  during subsetting.
 
 
 ## Changes in 0.2.7 (from 2026-03-27)

integration/test_sen2_analysis.py

@@ -7,7 +7,6 @@
 import xarray as xr
 
 from integration.helpers import assert_dataset_is_chunked
-from xarray_eopf.constants import DEFAULT_ENDPOINT_URL
 from xarray_eopf.utils import timeit
 
 allowed_open_time = 1000  # seconds

tests/amodes/test_sentinel3.py

@@ -14,7 +14,12 @@
 
 from tests.helpers import make_s3_olci_efr, make_s3_slstr_lst, make_s3_slstr_rbt
 from xarray_eopf.amode import AnalysisModeRegistry
-from xarray_eopf.amodes.sentinel3 import Sen3Ol1Efr, Sen3Sl1Rbt, Sen3Sl2Lst, register
+from xarray_eopf.amodes.sentinel3 import (
+    Sen3Ol1Efr,
+    Sen3Sl1Rbt,
+    Sen3Sl2Lst,
+    register,
+)
 from xarray_eopf.constants import FloatInt
 
 
@@ -91,7 +96,7 @@ def assert_convert_datatree_ok(
         self,
         original_dt: xr.DataTree,
         expected_var_names: list[str],
-        expected_size: (int, int),
+        expected_size: tuple[int, int],
         resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
         bbox: Sequence[float | int] | None = None,
     ):
@@ -118,6 +123,12 @@ def assert_convert_datatree_fail(self, original_dt: xr.DataTree):
         with pytest.raises(ValueError, match="No variables selected"):
             self.mode.convert_datatree(original_dt, includes="bibo")
 
+    def assert_convert_datatree_fail_with_include_exclude(
+        self, original_dt: xr.DataTree
+    ):
+        with pytest.raises(ValueError, match="No variables selected"):
+            self.mode.convert_datatree(original_dt, includes=".+", excludes=".+")
+
 
 class OlciEfrTest(Sen3TestMixin, TestCase):
     mode = Sen3Ol1Efr()
@@ -163,7 +174,7 @@ def test_convert_datatree_bbox(self):
                 "oa02_radiance",
                 "oa03_radiance",
             ],
-            expected_size=(372, 421),
+            expected_size=(372, 454),
             bbox=[1, 55, 3, 56],
         )
 
@@ -187,6 +198,21 @@ def test_convert_datatree_raise_warning(self):
     def test_convert_datatree_fail(self):
         self.assert_convert_datatree_fail(make_s3_olci_efr(size=48))
 
+    def test_convert_datatree_fail_include_exclude_overlap(self):
+        self.assert_convert_datatree_fail_with_include_exclude(
+            make_s3_olci_efr(size=48)
+        )
+
+    def test_convert_datatree_sets_other_metadata_as_attrs(self):
+        dt = make_s3_olci_efr(size=100)
+        dt.attrs["other_metadata"] = {"test_key": "test_val"}
+        ds = self.mode.convert_datatree(
+            dt,
+            includes=["oa01_radiance"],
+            resolution=0.1,
+        )
+        self.assertEqual({"test_key": "test_val"}, ds.attrs)
+
 
 class SlstrRbtTest(Sen3TestMixin, TestCase):
     mode = Sen3Sl1Rbt()
@@ -260,6 +286,11 @@ def test_convert_datatree_raise_warning(self):
     def test_convert_datatree_fail(self):
         self.assert_convert_datatree_fail(make_s3_slstr_rbt(size=48))
 
+    def test_convert_datatree_fail_include_exclude_overlap(self):
+        self.assert_convert_datatree_fail_with_include_exclude(
+            make_s3_slstr_rbt(size=48)
+        )
+
     def test_get_outer_bbox(self):
         bboxs = np.array([[-2, 10, 8, 20], [2, 12, 13, 25]])
         expected = [-2, 10, 13, 25]
@@ -307,9 +338,14 @@ def test_convert_datatree_bbox(self):
         self.assert_convert_datatree_ok(
             make_s3_slstr_lst(size=1000),
             expected_var_names=["lst"],
-            expected_size=(112, 127),
+            expected_size=(112, 148),
             bbox=[1, 55, 3, 56],
         )
 
     def test_convert_datatree_fail(self):
         self.assert_convert_datatree_fail(make_s3_slstr_lst(size=48))
+
+    def test_convert_datatree_fail_include_exclude_overlap(self):
+        self.assert_convert_datatree_fail_with_include_exclude(
+            make_s3_slstr_lst(size=48)
+        )

tests/helpers/sentinel3.py

@@ -11,25 +11,39 @@
 
 
 def make_s3_olci_efr(size: int = 48) -> xr.DataTree:
+    ds = make_s3_meas(size, bands=[f"oa{i:02}_radiance" for i in range(1, 22)])
 
-    return create_datatree(
-        {
-            "measurements": make_s3_meas(
-                size, bands=[f"oa{i:02}_radiance" for i in range(1, 22)]
-            ),
+    footprint = derive_footprint(ds)
+
+    dt = create_datatree(
+        {"measurements": ds},
+        attrs={
+            "stac_discovery": {
+                "geometry": {"type": "Polygon", "coordinates": [footprint]},
+                "properties": {"sat:orbit_state": "descending"},
+            }
         },
     )
+    return dt
 
 
 def make_s3_slstr_lst(size: int = 48) -> xr.DataTree:
+    ds = make_s3_meas(size, bands=["lst"])
+    footprint = derive_footprint(ds)
     return create_datatree(
         {
             "conditions/auxiliary": make_s3_meas(size, bands=["elevation"]),
             "conditions/meteorology": make_s3_meas((size, size // 10), bands=["s2m"]),
             "conditions/geometry": make_s3_meas(
                 (size, size // 10), bands=["sat_azimuth_tn", "sat_zenith_tn"]
             ),
-            "measurements": make_s3_meas(size, bands=["lst"]),
+            "measurements": ds,
+        },
+        attrs={
+            "stac_discovery": {
+                "geometry": {"type": "Polygon", "coordinates": [footprint]},
+                "properties": {"sat:orbit_state": "ascending"},
+            }
         },
     )
 
@@ -116,12 +130,30 @@ def make_coords(w: int, h: int, oblique_view=False) -> dict[str, xr.DataArray]:
     return {
         "latitude": xr.DataArray(lat_final, dims=("rows", "columns")),
         "longitude": xr.DataArray(lon_final, dims=("rows", "columns")),
-        "time_stamps": xr.DataArray(
-            np.arange(h).astype("datetime64[ns]"), dims=("rows")
-        ),
+        "time_stamps": xr.DataArray(np.arange(h).astype("datetime64[ns]"), dims="rows"),
     }
 
 
+def derive_footprint(ds: xr.Dataset) -> list[list[float]]:
+    lon = ds["longitude"]
+    lat = ds["latitude"]
+    corners = [
+        (0, 0),
+        (0, ds.sizes["columns"] - 1),
+        (ds.sizes["rows"] - 1, ds.sizes["columns"] - 1),
+        (ds.sizes["rows"] - 1, 0),
+    ]
+    footprint = [
+        [
+            float(lon.isel(rows=i, columns=j).values),
+            float(lat.isel(rows=i, columns=j).values),
+        ]
+        for i, j in corners
+    ]
+    footprint.append(footprint[0])
+    return footprint
+
+
 def create_datatree(
     datasets: dict[str, xr.Dataset], attrs: dict[str, Any] | None = None
 ) -> xr.DataTree:

tests/test_utils.py

@@ -3,16 +3,20 @@
 #  https://opensource.org/license/apache-2-0.
 
 from unittest import TestCase
+from unittest.mock import patch
 
+import numpy as np
 import pytest
 import xarray as xr
 
 from tests.helpers import make_s2_msi
 from xarray_eopf.utils import (
     NameFilter,
+    _find_relative_bbox,
     assert_arg_has_length,
     assert_arg_is_instance,
     assert_arg_is_one_of,
+    build_footprint_uv_mapping,
     get_data_tree_item,
     timeit,
 )
@@ -113,3 +117,49 @@ def test_filter(self):
         self.assertEqual(
             ["ernie", "emmie"], list(f.filter(["bibo", "ernie", "bert", "emmie"]))
         )
+
+
+class BuildFootprintUvMappingTest(TestCase):
+    def test_accepts_closed_ring_points(self):
+        open_ring = np.array(
+            [[10.0, 50.0], [12.0, 50.0], [12.0, 52.0], [10.0, 52.0]],
+            dtype=float,
+        )
+        closed_ring = np.vstack([open_ring, open_ring[0]])
+
+        open_xy, open_uv = build_footprint_uv_mapping(open_ring)
+        closed_xy, closed_uv = build_footprint_uv_mapping(closed_ring)
+
+        self.assertTrue(np.allclose(open_xy, closed_xy))
+        self.assertTrue(np.allclose(open_uv, closed_uv))
+
+    def test_find_relative_bbox_uses_southern_utm_epsg(self):
+        stac_meta = {
+            "geometry": {
+                "coordinates": [
+                    [
+                        [10.0, -11.0],
+                        [11.0, -11.0],
+                        [11.0, -10.0],
+                        [10.0, -10.0],
+                        [10.0, -11.0],
+                    ]
+                ]
+            },
+            "properties": {"sat:orbit_state": "descending"},
+        }
+        bbox = [10.2, -10.8, 10.8, -10.2]
+
+        with patch("xarray_eopf.utils.pyproj.Transformer.from_crs") as from_crs:
+            transformer = from_crs.return_value
+            transformer.transform.return_value = (
+                np.array([0.0, 1.0, 1.0, 0.0, 0.0]),
+                np.array([0.0, 0.0, 1.0, 1.0, 0.0]),
+            )
+            transformer.transform_bounds.return_value = (0.2, 0.2, 0.8, 0.8)
+
+            rel_bbox = _find_relative_bbox(stac_meta, bbox)
+
+        _, utm_epsg = from_crs.call_args.args[:2]
+        self.assertEqual("EPSG:32732", utm_epsg)
+        self.assertEqual(4, len(rel_bbox))