Add crs for Sen2 subgroups from STAC metadata

CHANGES.md

@@ -1,3 +1,9 @@
+## Changes in 0.2.7 (under development)
+
+- Add the CRS information from the STAC metadata stored in the datatree's attributes;
+  Temporally fixes the issue https://gitlab.eopf.copernicus.eu/cpm/eopf-cpm/-/issues/932
+
+
 ## Changes in 0.2.6 (from 2026-03-20)
 
 - Fixed an issue in `xr.open_dataset` (native mode) where selecting variables did not

tests/amodes/test_sentinel2.py

@@ -58,6 +58,7 @@ def test_assign_grid_mapping(self: TestCase):
         def make_band():
             return xr.DataArray(np.zeros((10, 10)), dims=("y", "x"))
 
+        # in dataset attrs with key "horizontal_CRS_code"
         dataset = self.mode.assign_grid_mapping(
             xr.Dataset(
                 dict(
@@ -66,7 +67,45 @@ def make_band():
                     b03=make_band(),
                 ),
                 attrs={"horizontal_CRS_code": "ESPG:32632"},
-            )
+            ),
+            {"test": "dict"},
+        )
+        self.assertIn("spatial_ref", dataset)
+        self.assertEqual(
+            "transverse_mercator", dataset.spatial_ref.attrs.get("grid_mapping_name")
+        )
+        self.assertEqual("spatial_ref", dataset.b01.attrs.get("grid_mapping"))
+        self.assertEqual("spatial_ref", dataset.b02.attrs.get("grid_mapping"))
+        self.assertEqual("spatial_ref", dataset.b03.attrs.get("grid_mapping"))
+
+        # in band attrs with key "proj:epsg"
+        ds = xr.Dataset(
+            dict(
+                b01=make_band(),
+                b02=make_band(),
+                b03=make_band(),
+            ),
+        )
+        ds["b01"].attrs = {"proj:epsg": 32632}
+        dataset = self.mode.assign_grid_mapping(ds, {"test": "dict"})
+        self.assertIn("spatial_ref", dataset)
+        self.assertEqual(
+            "transverse_mercator", dataset.spatial_ref.attrs.get("grid_mapping_name")
+        )
+        self.assertEqual("spatial_ref", dataset.b01.attrs.get("grid_mapping"))
+        self.assertEqual("spatial_ref", dataset.b02.attrs.get("grid_mapping"))
+        self.assertEqual("spatial_ref", dataset.b03.attrs.get("grid_mapping"))
+
+        # in band attrs with key "proj:epsg"
+        dataset = self.mode.assign_grid_mapping(
+            xr.Dataset(
+                dict(
+                    b01=make_band(),
+                    b02=make_band(),
+                    b03=make_band(),
+                ),
+            ),
+            {"properties": {"proj:code": "EPSG:32632"}},
         )
         self.assertIn("spatial_ref", dataset)
         self.assertEqual(
@@ -88,7 +127,8 @@ def make_band():
                     b03=make_band(),
                 ),
                 attrs={"horizontal_CRS_code": "ESPG:-1"},
-            )
+            ),
+            {},
         )
         self.assertNotIn("spatial_ref", dataset)
         self.assertEqual(None, dataset.b01.attrs.get("grid_mapping"))
@@ -110,7 +150,8 @@ def make_band():
                     b02=make_band(),
                     b03=make_band(),
                 )
-            )
+            ),
+            {},
         )
         self.assertNotIn("spatial_ref", dataset)
         self.assertEqual(None, dataset.b01.attrs.get("grid_mapping"))

tests/amodes/test_sentinel3.py

@@ -75,7 +75,7 @@ def test_assign_grid_mapping(self: TestCase):
 
     def test_transform_dataset(self: TestCase):
         dataset = self.mode.transform_dataset(
-            self.create_simple_dataset(), interp_method=0
+            self.create_simple_dataset(), {"stac_meta": "test"}, interp_method=0
         )
         self.assertIn("spatial_ref", dataset)
         self.assertEqual(

xarray_eopf/amode.py

@@ -95,7 +95,9 @@ def transform_datatree(self, datatree: xr.DataTree, **params) -> xr.DataTree:
         """
 
     @abstractmethod
-    def transform_dataset(self, dataset: xr.Dataset, **params) -> xr.Dataset:
+    def transform_dataset(
+        self, dataset: xr.Dataset, stac_meta: dict, **params
+    ) -> xr.Dataset:
         """Transform `dataset` into an analysis-ready form.
         Called from the backend's `open_dataset()` implementation to transform
         the given `xr.Dataset` object.
@@ -105,6 +107,7 @@ def transform_dataset(self, dataset: xr.Dataset, **params) -> xr.Dataset:
 
         Args:
             dataset: The dataset to be transformed.
+            stac_meta: The STAC metadata.
             params: Product type specific parameters.
                 See `get_applicable_params()`.
 

xarray_eopf/amodes/sentinel2.py

@@ -147,8 +147,10 @@ def transform_datatree(self, datatree: xr.DataTree, **params) -> xr.DataTree:
         )
         return datatree
 
-    def transform_dataset(self, dataset: xr.Dataset, **params) -> xr.Dataset:
-        return self.assign_grid_mapping(dataset)
+    def transform_dataset(
+        self, dataset: xr.Dataset, stac_meta: dict, **params
+    ) -> xr.Dataset:
+        return self.assign_grid_mapping(dataset, stac_meta)
 
     def convert_datatree(
         self,
@@ -221,7 +223,9 @@ def convert_datatree(
             if da_mapping:
                 ds = xr.Dataset(da_mapping)
                 ds.attrs.update(self.process_metadata(datatree))
-                datasets[res] = self.assign_grid_mapping(ds)
+                datasets[res] = self.assign_grid_mapping(
+                    ds, datatree.attrs.get("stac_discovery")
+                )
 
         # resample dataset
         if resolution in datasets and crs is None and bbox is None:
@@ -279,22 +283,30 @@ def convert_datatree(
         return rescaled_ds
 
     # noinspection PyMethodMayBeStatic
-    def assign_grid_mapping(self, dataset: xr.Dataset) -> xr.Dataset:
+    def assign_grid_mapping(self, dataset: xr.Dataset, stac_meta: dict) -> xr.Dataset:
         crs = None
         try:
             crs_code = dataset.attrs.get("horizontal_CRS_code", "EPSG:-1")
             epsg_int = int(crs_code.split(":")[1])
             crs = pyproj.CRS.from_epsg(epsg_int)
         except pyproj.exceptions.CRSError:
             pass
-        try:
-            for var_name, var in dataset.data_vars.items():
-                epsg_int = var.attrs.get("proj:epsg")
-                if isinstance(epsg_int, int):
-                    crs = pyproj.CRS.from_epsg(epsg_int)
-                    break
-        except pyproj.exceptions.CRSError:
-            pass
+        if crs is None:
+            try:
+                for var_name, var in dataset.data_vars.items():
+                    epsg_int = var.attrs.get("proj:epsg")
+                    if isinstance(epsg_int, int):
+                        crs = pyproj.CRS.from_epsg(epsg_int)
+                        break
+            except pyproj.exceptions.CRSError:
+                pass
+        if crs is None:
+            try:
+                crs_code = stac_meta.get("properties", {}).get("proj:code", "EPSG:-1")
+                epsg_int = int(crs_code.split(":")[1])
+                crs = pyproj.CRS.from_epsg(epsg_int)
+            except pyproj.exceptions.CRSError:
+                pass
 
         if crs:
             dataset = dataset.assign_coords(

xarray_eopf/amodes/sentinel3.py

@@ -85,7 +85,9 @@ def transform_datatree(self, datatree: xr.DataTree, **params) -> xr.DataTree:
         )
         return datatree
 
-    def transform_dataset(self, dataset: xr.Dataset, **params) -> xr.Dataset:
+    def transform_dataset(
+        self, dataset: xr.Dataset, stac_meta: dict, **params
+    ) -> xr.Dataset:
         return self.assign_grid_mapping(dataset)
 
     def convert_datatree(

xarray_eopf/backend.py

@@ -222,7 +222,9 @@ def open_dataset(
             if datatree.has_data:
                 # subgroup level, so we transform the dataset
                 dataset = datatree.to_dataset()
-                dataset = analysis_mode.transform_dataset(dataset)
+                dataset = analysis_mode.transform_dataset(
+                    dataset, datatree.attrs.get("stac_discovery")
+                )
             else:
                 # product level, so we convert the tree into a dataset
                 params = analysis_mode.get_applicable_params(

xarray_eopf/version.py

@@ -2,4 +2,4 @@
 #  Permissions are hereby granted under the terms of the Apache 2.0 License:
 #  https://opensource.org/license/apache-2-0.
 
-version = "0.2.6"
+version = "0.2.7dev0"