reading routine added for angles

environment.yml

@@ -16,6 +16,7 @@ dependencies:
   - pystac-client
   - shapely
   - xarray
+  - xmltodict
   - xcube >= 1.7.0
   - zarr >=2.11,<3  # until we can ensure zarr 3 compatibility; see Issue xcube #1102
   # for testing

pyproject.toml

@@ -29,6 +29,7 @@ dependencies = [
     "rioxarray",
     "requests",
     "shapely",
+    "xmltodict",
     "xarray",
     "xcube"
 ]

test/accessor/__init__.py

@@ -0,0 +1,20 @@
+# The MIT License (MIT)
+# Copyright (c) 2024 by the xcube development team and contributors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NON INFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.

test/test_accessor.py → test/accessor/test_sen2.py

@@ -19,17 +19,21 @@
 # OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
 # SOFTWARE.
 
+import json
 import unittest
 from unittest.mock import patch
 from unittest.mock import MagicMock
 
 import dask
 import dask.array as da
+import requests
+import pystac
 import xarray as xr
 import rasterio
 from xcube.core.mldataset import MultiLevelDataset
 
-from xcube_stac.accessor import S3Sentinel2DataAccessor
+from xcube_stac.accessor.sen2 import S3Sentinel2DataAccessor
+from ..sampledata import sentinel_2_band_data
 
 
 class TestS3Sentinel2DataAccessor(unittest.TestCase):
@@ -119,3 +123,18 @@ def test_open_data(self, mock_rioxarray_open, mock_rasterio_open):
         self.assertCountEqual(
             [1024, 1024], [ds.chunksizes["x"][0], ds.chunksizes["y"][0]]
         )
+
+    # @pytest.mark.vcr()
+    def test_add_sen2_angles(self):
+        ds = sentinel_2_band_data()
+        item_url = (
+            "https://stac.dataspace.copernicus.eu/v1/collections/sentinel-2-l2a/"
+            "items/S2A_MSIL2A_20200301T090901_N0500_R050_T35UPU_20230630T033416"
+        )
+        response = requests.request(method="GET", url=item_url)
+        item = pystac.Item.from_dict(
+            json.loads(response.text), href=item_url, preserve_dict=False
+        )
+        ds_test = self.accessor.add_sen2_angles(item, ds)
+        self.assertIn("solar_angle", ds_test)
+        self.assertIn("viewing_angle", ds_test)

test/sampledata.py

@@ -0,0 +1,84 @@
+# The MIT License (MIT)
+# Copyright (c) 2024 by the xcube development team and contributors
+#
+# Permission is hereby granted, free of charge, to any person obtaining a
+# copy of this software and associated documentation files (the "Software"),
+# to deal in the Software without restriction, including without limitation
+# the rights to use, copy, modify, merge, publish, distribute, sublicense,
+# and/or sell copies of the Software, and to permit persons to whom the
+# Software is furnished to do so, subject to the following conditions:
+#
+# The above copyright notice and this permission notice shall be included in all
+# copies or substantial portions of the Software.
+#
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+# FITNESS FOR A PARTICULAR PURPOSE AND NON INFRINGEMENT. IN NO EVENT SHALL THE
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+# SOFTWARE.
+
+import dask.array as da
+import numpy as np
+import xarray as xr
+
+
+def sentinel_2_band_data():
+    mock_data = {
+        "band_1": (
+            ("y", "x"),
+            da.ones((10980, 10980), chunks=(1024, 1024), dtype=np.uint16),
+        ),
+    }
+    spatial_ref = xr.DataArray(
+        np.array(0),
+        attrs={
+            "crs_wkt": (
+                'PROJCS["WGS 84 / UTM zone 35N",GEOGCS["WGS 84",DATUM["WGS_1984",'
+                'SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],'
+                'AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG"'
+                ',"8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",'
+                '"9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_'
+                'Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central'
+                '_meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false'
+                '_easting",500000],PARAMETER["false_northing",0],UNIT["metre"'
+                ',1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing"'
+                ',NORTH],AUTHORITY["EPSG","32635"]]'
+            ),
+            "semi_major_axis": 6378137.0,
+            "semi_minor_axis": 6356752.314245179,
+            "inverse_flattening": 298.257223563,
+            "reference_ellipsoid_name": "WGS 84",
+            "longitude_of_prime_meridian": 0.0,
+            "prime_meridian_name": "Greenwich",
+            "geographic_crs_name": "WGS 84",
+            "horizontal_datum_name": "World Geodetic System 1984",
+            "projected_crs_name": "WGS 84 / UTM zone 35N",
+            "grid_mapping_name": "transverse_mercator",
+            "latitude_of_projection_origin": 0.0,
+            "longitude_of_central_meridian": 27.0,
+            "false_easting": 500000.0,
+            "false_northing": 0.0,
+            "scale_factor_at_central_meridian": 0.9996,
+            "spatial_ref": (
+                'PROJCS["WGS 84 / UTM zone 35N",GEOGCS["WGS 84",DATUM["WGS_1984",'
+                'SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]'
+                ',AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG"'
+                ',"8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",'
+                '"9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_'
+                'Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_'
+                'meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false_'
+                'easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,'
+                'AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",'
+                'NORTH],AUTHORITY["EPSG","32635"]]'
+            ),
+            "GeoTransform": "600000.0 10.0 0.0 5900040.0 0.0 -10.0",
+        },
+    )
+    coords = {
+        "x": np.arange(600005.0, 709796.0, 10.0),
+        "y": np.arange(5900035.0, 5790244.0, -10.0),
+        "spatial_ref": spatial_ref,
+    }
+    return xr.Dataset(mock_data, coords=coords)

test/test_stack.py

@@ -25,7 +25,7 @@
 import pystac
 import xarray as xr
 
-from xcube_stac.stack import mosaic_3d_take_first
+from xcube_stac.stack import mosaic_spatial_along_time_take_first
 from xcube_stac.stack import groupby_solar_day
 
 
@@ -119,14 +119,14 @@ def test_mosaic_3d_take_first(self):
 
         # test only one tile
         dts = np.array(["2025-01-01", "2025-01-02", "2025-01-03"], dtype="datetime64")
-        ds_test = mosaic_3d_take_first(list_ds[:1], dts)
+        ds_test = mosaic_spatial_along_time_take_first(list_ds[:1], dts)
         xr.testing.assert_allclose(ds_test, list_ds[0])
 
         # test two tiles
         dts = np.array(
             ["2025-01-01", "2025-01-02", "2025-01-03", "2025-01-04"], dtype="datetime64"
         )
-        ds_test = mosaic_3d_take_first(list_ds, dts)
+        ds_test = mosaic_spatial_along_time_take_first(list_ds, dts)
         data = np.array(
             [
                 [[1, 2, 3], [4, 5, 6], [7, 8, 9]],
@@ -155,5 +155,5 @@ def test_mosaic_3d_take_first(self):
             list_ds[i] = ds
         ds_expected = xr.Dataset({"B01": da})
         ds_expected = ds_expected.assign_coords({"spatial_ref": spatial_ref})
-        ds_test = mosaic_3d_take_first(list_ds, dts)
+        ds_test = mosaic_spatial_along_time_take_first(list_ds, dts)
         xr.testing.assert_allclose(ds_test, ds_expected)

test/test_store.py

@@ -41,9 +41,10 @@
 from xcube_stac.constants import DATA_STORE_ID
 from xcube_stac.constants import DATA_STORE_ID_XCUBE
 from xcube_stac.constants import DATA_STORE_ID_CDSE
-from xcube_stac.constants import SENITNEL_2_L2A_BANDS
-from xcube_stac.accessor import HttpsDataAccessor
-from xcube_stac.accessor import S3DataAccessor
+from xcube_stac.constants import SENITNEL2_L2A_BANDS
+from xcube_stac.accessor.https import HttpsDataAccessor
+from xcube_stac.accessor.s3 import S3DataAccessor
+import sampledata
 
 SKIP_HELP = (
     "Skipped, because server is not running:"
@@ -646,19 +647,19 @@ def test_open_data_stack_mode(self):
         store = new_data_store(DATA_STORE_ID, url=self.url_searchable, stack_mode=True)
 
         # open data as dataset
-        bbox_utm = [659574, 5892990, 659724, 5893140]
+        bbox_utm = [5599905, 3511735, 5600064, 3511894]
         ds = store.open_data(
             data_id="sentinel-2-l2a",
             bbox=bbox_utm,
             time_range=["2023-11-01", "2023-11-10"],
             spatial_res=10,
-            crs="EPSG:32635",
+            crs="EPSG:3035",
             asset_names=["red", "green", "blue"],
             apply_scaling=True,
             open_params_dataset_geotiff=dict(tile_size=(512, 512)),
         )
         self.assertIsInstance(ds, xr.Dataset)
-        self.assertCountEqual(["red", "green", "blue"], list(ds.data_vars))
+        self.assertCountEqual(["red", "green", "blue", "crs"], list(ds.data_vars))
         self.assertCountEqual(
             [4, 16, 16],
             [ds.sizes["time"], ds.sizes["y"], ds.sizes["x"]],
@@ -699,64 +700,7 @@ def test_open_data_stack_mode_no_items_found(self):
     @pytest.mark.vcr()
     @patch("rioxarray.open_rasterio")
     def test_open_data_stack_mode_cdse_sen2(self, mock_rioxarray_open):
-        mock_data = {
-            "band_1": (
-                ("y", "x"),
-                da.ones((10980, 10980), chunks=(1024, 1024), dtype=np.uint16),
-            ),
-        }
-        spatial_ref = xr.DataArray(
-            np.array(0),
-            attrs={
-                "crs_wkt": (
-                    'PROJCS["WGS 84 / UTM zone 35N",GEOGCS["WGS 84",DATUM["WGS_1984",'
-                    'SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],'
-                    'AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG"'
-                    ',"8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",'
-                    '"9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_'
-                    'Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central'
-                    '_meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false'
-                    '_easting",500000],PARAMETER["false_northing",0],UNIT["metre"'
-                    ',1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing"'
-                    ',NORTH],AUTHORITY["EPSG","32635"]]'
-                ),
-                "semi_major_axis": 6378137.0,
-                "semi_minor_axis": 6356752.314245179,
-                "inverse_flattening": 298.257223563,
-                "reference_ellipsoid_name": "WGS 84",
-                "longitude_of_prime_meridian": 0.0,
-                "prime_meridian_name": "Greenwich",
-                "geographic_crs_name": "WGS 84",
-                "horizontal_datum_name": "World Geodetic System 1984",
-                "projected_crs_name": "WGS 84 / UTM zone 35N",
-                "grid_mapping_name": "transverse_mercator",
-                "latitude_of_projection_origin": 0.0,
-                "longitude_of_central_meridian": 27.0,
-                "false_easting": 500000.0,
-                "false_northing": 0.0,
-                "scale_factor_at_central_meridian": 0.9996,
-                "spatial_ref": (
-                    'PROJCS["WGS 84 / UTM zone 35N",GEOGCS["WGS 84",DATUM["WGS_1984",'
-                    'SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]'
-                    ',AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG"'
-                    ',"8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG",'
-                    '"9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Transverse_'
-                    'Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_'
-                    'meridian",27],PARAMETER["scale_factor",0.9996],PARAMETER["false_'
-                    'easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,'
-                    'AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",'
-                    'NORTH],AUTHORITY["EPSG","32635"]]'
-                ),
-                "GeoTransform": "600000.0 10.0 0.0 5900040.0 0.0 -10.0",
-            },
-        )
-        coords = {
-            "x": np.arange(600005.0, 709796.0, 10.0),
-            "y": np.arange(5900035.0, 5790244.0, -10.0),
-            "spatial_ref": spatial_ref,
-        }
-        mock_ds = xr.Dataset(mock_data, coords=coords)
-        mock_rioxarray_open.return_value = mock_ds
+        mock_rioxarray_open.return_value = sampledata.sentinel_2_band_data()
 
         store = new_data_store(
             DATA_STORE_ID_CDSE,
@@ -766,25 +710,45 @@ def test_open_data_stack_mode_cdse_sen2(self, mock_rioxarray_open):
         )
 
         # open data as dataset
-        bbox_utm = [659574, 5892990, 659724, 5893140]
+        bbox_utm = [620000, 5800000, 680000, 5860000]
         ds = store.open_data(
             data_id="sentinel-2-l2a",
             bbox=bbox_utm,
             time_range=["2023-11-01", "2023-11-10"],
-            spatial_res=10,
+            spatial_res=20,
             crs="EPSG:32635",
             apply_scaling=True,
+            angles_sentinel2=True,
         )
         self.assertIsInstance(ds, xr.Dataset)
 
-        self.assertCountEqual(SENITNEL_2_L2A_BANDS + ["crs"], list(ds.data_vars))
         self.assertCountEqual(
-            [4, 16, 16],
-            [ds.sizes["time"], ds.sizes["y"], ds.sizes["x"]],
+            SENITNEL2_L2A_BANDS + ["crs", "solar_angle", "viewing_angle"],
+            list(ds.data_vars),
         )
         self.assertCountEqual(
-            [1, 16, 16],
-            [ds.chunksizes["time"][0], ds.chunksizes["y"][0], ds.chunksizes["x"][0]],
+            [4, 3001, 3001, 13, 13, 2, 12],
+            [
+                ds.sizes["time"],
+                ds.sizes["y"],
+                ds.sizes["x"],
+                ds.sizes["angle_y"],
+                ds.sizes["angle_x"],
+                ds.sizes["angle"],
+                ds.sizes["band"],
+            ],
+        )
+        self.assertCountEqual(
+            [1, 1024, 1024, 13, 13, 2, 12],
+            [
+                ds.chunksizes["time"][0],
+                ds.chunksizes["y"][0],
+                ds.chunksizes["x"][0],
+                ds.chunksizes["angle_y"][0],
+                ds.chunksizes["angle_x"][0],
+                ds.chunksizes["angle"][0],
+                ds.chunksizes["band"][0],
+            ],
         )
 
     @pytest.mark.vcr()

test/test_utils.py

@@ -312,8 +312,8 @@ def test_reproject_bbox(self):
             5768595.563692021,
         ]
         np.testing.assert_almost_equal(
-            reproject_bbox(bbox_utm, crs_utm, crs_wgs84),
-            [178.8245008, 49.9483787, -178.9158425, 52.0509362],
+            reproject_bbox(bbox_utm, crs_utm, crs_wgs84, buffer=0.02),
+            [178.77930769, 49.90632759, -178.87064939, 52.09298731],
         )
 
     def test_normalize_crs(self):