Small adjustments for xcube-eopf Sentinel-1 Implementation

CHANGES.md

@@ -1,4 +1,4 @@
-## Changes in 0.3.0 (in development)
+## Changes in 0.3.0 (from 2026-07-22)
 
 - Sentinel-1 GRD analysis mode is now fully lazy, enabling seamless execution on 
   local and distributed Dask clusters.

tests/amodes/test_sentinel1.py

@@ -144,6 +144,10 @@ def test_get_applicable_params(self: TestCase):
                 cache_uri="file:///tmp/cache",
             ),
         )
+        with pytest.raises(
+            TypeError, match="resolution argument must contain exactly two"
+        ):
+            self.mode.get_applicable_params(resolution=(1, "x"))
         with pytest.raises(TypeError, match="interp_methods"):
             self.mode.get_applicable_params(interp_methods="cubic")
         with pytest.raises(TypeError, match="footprint_scale_factor"):
@@ -413,9 +417,12 @@ def test_is_not_valid_source(self):
 
     def test_get_groups(self):
         groups = self.mode._get_groups(self.dt)
-        self.assertEqual(("mode", "swath", "burst"), groups.dims)
-        self.assertEqual((2, 2, 1), groups.shape)
-        self.assertEqual("S1A_IW_SLC_TEST_VV_IW1_0", groups.sel(mode="VV").item(0))
+        self.assertEqual(2, len(groups))
+        self.assertEqual(2, len(groups[0]))
+        self.assertEqual(3, len(groups[0][0]))
+        self.assertEqual(1, len(groups[0][0][0]))
+        self.assertEqual("S1A_IW_SLC_TEST_VV_IW1_0", groups[0][0][0][0])
+        self.assertEqual(["VH", "VV"], [str(mode) for mode in groups[1]])
 
     def test_get_grid_parameters(self):
         params = sen1._get_grid_parameters(
@@ -441,7 +448,7 @@ def test_open_data(self):
         self.assertCountEqual(["beta0_vv", "beta0_vh"], out.data_vars)
         self.assertNotIn("line", out.coords)
         self.assertNotIn("pixel", out.coords)
-        self.assertEqual({"azimuth_time": 2, "slant_range_time": 3}, out.sizes)
+        self.assertEqual({"azimuth_time": 2, "slant_range_time": 4}, out.sizes)
 
     def test_open_data_accepts_string_include(self):
         out = self.mode._open_data(self.dt, includes="vv")
@@ -670,6 +677,10 @@ def test_get_applicable_params(self: TestCase):
                 agg_methods="nearest",
             ),
         )
+        with pytest.raises(
+            TypeError, match="resolution argument must contain exactly two"
+        ):
+            self.mode.get_applicable_params(resolution=(1, "x"))
         with pytest.raises(TypeError):
             self.mode.get_applicable_params(interp_methods="cubic")
 
@@ -823,7 +834,7 @@ def test_get_dem_requires_credentials(self):
             with pytest.raises(ValueError, match="Missing AWS credentials"):
                 sen1.get_dem([0, 50, 1, 51])
 
-    def test_get_dem_resolution_required_if_crs_given(self):
+    def test_get_dem_with_projected_crs_uses_inferred_resolution(self):
         with patch.dict(
             os.environ,
             {"AWS_ACCESS_KEY_ID": "k", "AWS_SECRET_ACCESS_KEY": "s"},
@@ -832,6 +843,18 @@ def test_get_dem_resolution_required_if_crs_given(self):
             with (
                 patch.object(sen1.pystac_client.Client, "open") as client_open,
                 patch.object(sen1.rioxarray, "open_rasterio") as open_rasterio,
+                patch.object(sen1, "transform_resolution", return_value=30.0) as tr,
+                patch.object(
+                    sen1,
+                    "resample_in_space",
+                    return_value=SimpleNamespace(
+                        dem=xr.DataArray(
+                            np.ones((2, 2), dtype="float32"),
+                            dims=("lat", "lon"),
+                            coords={"lat": [1.0, 0.0], "lon": [0.0, 1.0]},
+                        )
+                    ),
+                ) as resample,
             ):
                 fake_item = SimpleNamespace(assets={"data": SimpleNamespace(href="x")})
                 search = SimpleNamespace(items=lambda: [fake_item])
@@ -842,10 +865,12 @@ def test_get_dem_resolution_required_if_crs_given(self):
                     coords={"band": [1], "y": [3, 2, 1, 0], "x": [0, 1, 2, 3]},
                 )
 
-                with pytest.raises(
-                    ValueError, match="Resolution must be provided if CRS is not None"
-                ):
-                    sen1.get_dem([0, 0, 1, 1], crs=pyproj.CRS.from_epsg(32632))
+                out = sen1.get_dem([0, 0, 900, 900], crs=pyproj.CRS.from_epsg(32632))
+
+        tr.assert_called_once()
+        resample.assert_called_once()
+        self.assertIsInstance(out, xr.DataArray)
+        self.assertEqual((2, 2), out.shape)
 
     def test_get_dem_reprojects_bbox_and_resamples(self):
         with patch.dict(

tests/amodes/test_sentinel3.py

@@ -49,6 +49,12 @@ def test_get_applicable_params(self: TestCase):
                 agg_methods={"scl": "mode"},
             ),
         )
+        with pytest.raises(
+            TypeError, match="resolution argument must contain exactly two"
+        ):
+            self.mode.get_applicable_params(resolution=(1, "x"))
+        with pytest.raises(TypeError):
+            self.mode.get_applicable_params(interp_methods=["nearest"])
 
     def test_process_metadata(self: TestCase):
         self.assertEqual({}, self.mode.process_metadata(xr.DataTree()))

tests/helpers/sentinel1.py

@@ -128,14 +128,23 @@ def make_s1_slc_datatree() -> xr.DataTree:
             ],
             dtype="datetime64[ns]",
         ),
+        "IW3": np.array(
+            [
+                "2024-01-01T00:00:00",
+                "2024-01-01T00:00:01",
+                "2024-01-01T00:00:02",
+            ],
+            dtype="datetime64[ns]",
+        ),
     }
     slant_range_time_by_swath = {
         "IW1": np.array([0.0, 1.0, 2.0, 3.0], dtype="float32"),
         "IW2": np.array([1.0, 2.0, 3.0, 4.0], dtype="float32"),
+        "IW3": np.array([2.0, 3.0, 4.0, 5.0], dtype="float32"),
     }
     axis = np.array(["x", "y", "z"])
     polarizations = ("VV", "VH")
-    swaths = ("IW1", "IW2")
+    swaths = ("IW1", "IW2", "IW3")
 
     beta_nought = xr.Dataset(
         {

xarray_eopf/amodes/sentinel1.py

@@ -35,6 +35,7 @@
 )
 
 from xarray_eopf.amode import AnalysisMode, AnalysisModeRegistry
+from xarray_eopf.constants import FloatInt
 from xarray_eopf.source import get_source_path
 from xarray_eopf.utils import NameFilter, assert_arg_has_length, assert_arg_is_instance
 
@@ -45,6 +46,8 @@
 _CRS_WGS84 = pyproj.CRS.from_string("EPSG:4326")
 _DEM_CHUNKSIZE = dict(lat=1800, lon=1800)
 _CHUNKSIZE = (2048, 2048)
+_SLC_SWATHS = ["IW1", "IW2", "IW3"]
+_SENTINEL1_POLARZIATION_MODES = ["VV", "VH", "HV", "HH"]
 _REGISTERED_CACHE_URIS = []
 
 
@@ -111,7 +114,14 @@ def get_applicable_params(self, **kwargs) -> dict[str, Any]:
 
         resolution = kwargs.get("resolution")
         if resolution is not None:
-            assert_arg_is_instance(resolution, "resolution", (float, int))
+            assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+            if isinstance(resolution, tuple):
+                assert_arg_has_length(resolution, "resolution", 2)
+                if not all(isinstance(v, (float, int)) for v in resolution):
+                    raise TypeError(
+                        "resolution argument must contain exactly "
+                        "two float or int values."
+                    )
             params.update(resolution=resolution)
 
         bbox = kwargs.get("bbox")
@@ -171,7 +181,7 @@ def convert_datatree(
         datatree: xr.DataTree,
         includes: str | Iterable[str] | None = None,
         excludes: str | Iterable[str] | None = None,
-        resolution: float = None,
+        resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
         bbox: Sequence[float | int] | None = None,
         crs: pyproj.CRS | None = None,
         interp_methods: Literal["nearest", "bilinear"] = "bilinear",
@@ -231,7 +241,7 @@ def _open_data(
 
         dataset = None
         measurement_group = ""
-        for mode in ["VV", "VH", "HV", "HH"]:
+        for mode in _SENTINEL1_POLARZIATION_MODES:
             children = [x for x in datatree.children if mode in x]
             if children:
                 if len(children) != 1:
@@ -359,26 +369,28 @@ def _open_data(
     ) -> xr.Dataset:
         """Load, calibrate, and merge SLC bursts into a single dataset."""
 
-        children = self._get_groups(datatree)
+        children, modes = self._get_groups(datatree)
 
         # First combine all polarizations within each burst.
-        dss_burst = np.empty(children.shape[1:], dtype=object)
-        for swath_i, swath in enumerate(children.swath.values):
-            for burst_idx in children.burst.values:
+        dss_all = []
+        for swath_i, _ in enumerate(_SLC_SWATHS):
+            dss_swath = []
+            for burst_i, _ in enumerate(children[0][swath_i]):
                 burst_dataset = xr.Dataset()
-                for mode_i, mode in enumerate(children.mode.values):
-                    child = children.sel(mode=mode, swath=swath, burst=burst_idx).item()
+                for mode_i, mode in enumerate(modes):
+                    child = children[mode_i][swath_i][burst_i]
                     burst = datatree[child]
                     beta0 = self._calibrate_burst(burst)
                     beta0 = self._extract_valid_region(beta0, burst)
                     beta0 = beta0.drop_vars(["line", "pixel"])
                     beta0 = beta0.rename({"slc": f"beta0_{mode.lower()}"})
                     burst_dataset.update(beta0)
-                dss_burst[swath_i, burst_idx] = burst_dataset
+                dss_swath.append(burst_dataset)
+            dss_all.append(dss_swath)
 
-        dss_swaths = np.empty(children.shape[1], dtype=object)
-        for swath_i, swath in enumerate(children.swath.values):
-            dss_swaths[swath_i] = self._merge_bursts(dss_burst[swath_i, :])
+        dss_swaths = []
+        for swath_i, swath in enumerate(_SLC_SWATHS):
+            dss_swaths.append(self._merge_bursts(dss_all[swath_i]))
         dss_swaths = self._align_azimuth(dss_swaths)
         merged = self._merge_swaths(dss_swaths)
 
@@ -407,32 +419,20 @@ def _expand_names(
         return merged
 
     @staticmethod
-    def _get_groups(datatree):
+    def _get_groups(datatree: xr.DataTree) -> tuple[list[list[list[str]]], list[str]]:
         """Return child group names organized by polarization, swath, and burst."""
-        swaths = ["IW1", "IW2", "IW3"]
-        modes = ["VV", "VH", "HV", "HH"]
+        modes_sel = []
         children = []
-        for mode in modes:
+        for mode in _SENTINEL1_POLARZIATION_MODES:
             mode_children = []
-            for swath_i, swath in enumerate(swaths):
+            for swath_i, swath in enumerate(_SLC_SWATHS):
                 bursts = [x for x in datatree.children if f"_{mode}_{swath}" in x]
                 if bursts:
                     mode_children.append(bursts)
+                    modes_sel.append(mode)
             if mode_children:
                 children.append(mode_children)
-        children = np.array(children, dtype=str)
-        return xr.DataArray(
-            children,
-            dims=("mode", "swath", "burst"),
-            coords={
-                "mode": [
-                    m for m in modes if any(f"_{m}_" in x for x in datatree.children)
-                ],
-                "swath": swaths[: children.shape[1]],
-                "burst": np.arange(children.shape[2]),
-            },
-            name="burst_groups",
-        )
+        return children, list(np.unique(modes_sel))
 
     @staticmethod
     def _calibrate_burst(burst: xr.DataTree) -> xr.Dataset:
@@ -460,7 +460,7 @@ def _extract_valid_region(dataset: xr.Dataset, burst: xr.DataTree) -> xr.Dataset
         )
 
     @staticmethod
-    def _merge_bursts(dss: np.ndarray, tolerance: float = 0.01) -> xr.Dataset:
+    def _merge_bursts(dss: list[xr.Dataset], tolerance: float = 0.01) -> xr.Dataset:
         """Merge overlapping bursts along azimuth time."""
         idxs = np.zeros((len(dss), 2), dtype=int)
         tol = 0.01 * np.diff(dss[0].azimuth_time.values[:2])[0]
@@ -517,7 +517,9 @@ def _merge_bursts(dss: np.ndarray, tolerance: float = 0.01) -> xr.Dataset:
         return out
 
     @staticmethod
-    def _align_azimuth(dss: np.ndarray, tolerance: float = 0.01) -> list[xr.Dataset]:
+    def _align_azimuth(
+        dss: list[xr.Dataset], tolerance: float = 0.01
+    ) -> list[xr.Dataset]:
         """Align all swaths to a shared azimuth-time axis."""
         tol = 0.1 * np.diff(dss[0].azimuth_time.values[:2])[0]
 
@@ -613,7 +615,14 @@ def get_applicable_params(self, **kwargs) -> dict[str, Any]:
 
         resolution = kwargs.get("resolution")
         if resolution is not None:
-            assert_arg_is_instance(resolution, "resolution", (float, int))
+            assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+            if isinstance(resolution, tuple):
+                assert_arg_has_length(resolution, "resolution", 2)
+                if not all(isinstance(v, (float, int)) for v in resolution):
+                    raise TypeError(
+                        "resolution argument must contain exactly "
+                        "two float or int values."
+                    )
             params.update(resolution=resolution)
 
         bbox = kwargs.get("bbox")
@@ -648,7 +657,7 @@ def convert_datatree(
         datatree: xr.DataTree,
         includes: str | Iterable[str] | None = None,
         excludes: str | Iterable[str] | None = None,
-        resolution: float = None,
+        resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
         bbox: Sequence[float | int] | None = None,
         crs: pyproj.CRS | None = None,
         interp_methods: SpatialInterpMethods | None = None,
@@ -708,12 +717,12 @@ def _apply_valid_range(array, *, dtype=None, fill_value=None):
                 bbox = source_gm.xy_bbox
         if resolution is None:
             if crs and not crs.is_geographic:
-                center_lat = (
+                ref_point = (
                     (source_gm.xy_bbox[0] + source_gm.xy_bbox[2]) / 2,
                     (source_gm.xy_bbox[1] + source_gm.xy_bbox[3]) / 2,
                 )
                 resolution = transform_resolution(
-                    center_lat, source_gm.xy_res, source_gm.crs, crs
+                    ref_point, source_gm.xy_res, source_gm.crs, crs
                 )
             else:
                 resolution = source_gm.xy_res
@@ -753,7 +762,7 @@ def _cleanup_registered_cache_uris() -> None:
 
 def get_dem(
     bbox: Sequence[float | int],
-    resolution: float | None = None,
+    resolution: FloatInt | tuple[FloatInt, FloatInt] | None = None,
     crs: pyproj.CRS | None = None,
 ):
     """Fetch and prepare a DEM for the given area of interest.
@@ -817,17 +826,25 @@ def get_dem(
             lon=slice(bbox_wgs84[0], bbox_wgs84[2]),
         ).chunk(_DEM_CHUNKSIZE)
     else:
-        if resolution is None:
-            raise ValueError("Resolution must be provided if CRS is not None.")
+        dem = dem.to_dataset(name="dem")
         if crs is None:
             crs = _CRS_WGS84
+        if resolution is None:
+            source_gm = GridMapping.from_dataset(dem)
+            ref_point = (
+                (bbox_wgs84[0] + bbox_wgs84[2]) / 2,
+                (bbox_wgs84[1] + bbox_wgs84[3]) / 2,
+            )
+            resolution = transform_resolution(
+                ref_point, source_gm.xy_res, source_gm.crs, crs
+            )
         target_gm = GridMapping.regular_from_bbox(
             bbox,
             resolution,
             crs,
             tile_size=(_DEM_CHUNKSIZE["lat"], _DEM_CHUNKSIZE["lon"]),
         )
-        dem = resample_in_space(dem.to_dataset(name="dem"), target_gm=target_gm).dem
+        dem = resample_in_space(dem, target_gm=target_gm).dem
 
     return dem
 

xarray_eopf/amodes/sentinel3.py

@@ -56,7 +56,14 @@ def get_applicable_params(self, **kwargs) -> dict[str, Any]:
 
         resolution = kwargs.get("resolution")
         if resolution is not None:
-            assert_arg_is_instance(resolution, "resolution", (int, float))
+            assert_arg_is_instance(resolution, "resolution", (float, int, tuple))
+            if isinstance(resolution, tuple):
+                assert_arg_has_length(resolution, "resolution", 2)
+                if not all(isinstance(v, (float, int)) for v in resolution):
+                    raise TypeError(
+                        "resolution argument must contain exactly "
+                        "two float or int values."
+                    )
             params.update(resolution=resolution)
 
         bbox = kwargs.get("bbox")