diff --git a/CHANGES.md b/CHANGES.md index 4209517..1002f83 100644 --- a/CHANGES.md +++ b/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 diff --git a/tests/amodes/test_sentinel2.py b/tests/amodes/test_sentinel2.py index 2a52ee1..120b92d 100644 --- a/tests/amodes/test_sentinel2.py +++ b/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")) diff --git a/tests/amodes/test_sentinel3.py b/tests/amodes/test_sentinel3.py index 56c0799..e4e735c 100644 --- a/tests/amodes/test_sentinel3.py +++ b/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( diff --git a/xarray_eopf/amode.py b/xarray_eopf/amode.py index 7f4a555..1e742e0 100644 --- a/xarray_eopf/amode.py +++ b/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()`. diff --git a/xarray_eopf/amodes/sentinel2.py b/xarray_eopf/amodes/sentinel2.py index 3c06789..b120808 100644 --- a/xarray_eopf/amodes/sentinel2.py +++ b/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,7 +283,7 @@ 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") @@ -287,14 +291,22 @@ def assign_grid_mapping(self, dataset: xr.Dataset) -> xr.Dataset: 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( diff --git a/xarray_eopf/amodes/sentinel3.py b/xarray_eopf/amodes/sentinel3.py index 350c34b..ecb13b7 100644 --- a/xarray_eopf/amodes/sentinel3.py +++ b/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( diff --git a/xarray_eopf/backend.py b/xarray_eopf/backend.py index 3e9923d..2fefbac 100644 --- a/xarray_eopf/backend.py +++ b/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( diff --git a/xarray_eopf/version.py b/xarray_eopf/version.py index ee5c7a0..ff3b410 100644 --- a/xarray_eopf/version.py +++ b/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"