Fix wl_bounds and unstack_dates (#359)

<!-- Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #360 
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features).
- [x] (If applicable) Tests have been added.
- [x] This PR does not seem to break the templates.
- [x] CHANGES.rst has been updated (with summary of main changes).
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added.

### What kind of change does this PR introduce?

* Fix `unstack_dates` where a "A" freq base was still to be found.
* Fix `get_warming_levels` when used with a dataset that already has
bounds (like for the grid). The `warminglevel_bounds` variable now has
the dimension "wl_bounds" instead of "bounds" to avoid name conflicts.

### Does this PR introduce a breaking change?
No. Unless someone was using the `bounds` dimension ?

### Other information:

CHANGES.rst

@@ -4,7 +4,7 @@ Changelog
 
 v0.9.0 (unreleased)
 -------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Gabriel Rondeau-Genesse (:user:`RondeauG`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^
@@ -13,6 +13,11 @@ Internal changes
     * Added a few free `grep`-based hooks for finding unwanted artifacts in the code base.
     * Updated `ruff` to v0.2.0 and `black` to v24.2.0.
 
+Bug fixes
+^^^^^^^^^
+* Fix ``unstack_dates`` for the new frequency syntax introduced by pandas v2.2. (:pull:`359`).
+* ``subset_warming_level`` will not return partial subsets if the warming level is reached at the end of the timeseries. (:issue:`360`, :pull:`359`).
+
 v0.8.3 (2024-02-28)
 -------------------
 Contributors to this version: Juliette Lavoie (:user:`juliettelavoie`), Trevor James Smith (:user:`Zeitsperre`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Pascal Bourgault (:user:`aulemahal`).

tests/test_extract.py

@@ -405,7 +405,7 @@ class TestSubsetWarmingLevel:
         np.tile(np.arange(1, 2), 50),
         variable="tas",
         start="2000-01-01",
-        freq="AS-JAN",
+        freq="YS-JAN",
         as_dataset=True,
     ).assign_attrs(
         {
@@ -440,11 +440,20 @@ def test_kwargs(self):
         assert ds_sub.warminglevel.attrs["baseline"] == "1981-2010"
         assert ds_sub.attrs["cat:processing_level"] == "tests"
 
-    def test_outofrange(self):
-        assert xs.subset_warming_level(self.ds, wl=5) is None
-
-    def test_none(self):
-        assert xs.subset_warming_level(self.ds, wl=20) is None
+    @pytest.mark.parametrize("wl", [3.5, 5, 20, [2, 3.5, 5, 20], [3.5, 5, 20]])
+    def test_outofrange(self, wl):
+        # 3.5 is only partially covered by ds, 5 is out of range but within the csv, 20 is fully out of range
+        if not isinstance(wl, list):
+            assert xs.subset_warming_level(self.ds, wl=wl) is None
+        else:
+            ds = xs.subset_warming_level(self.ds, wl=wl)
+            assert ds.warminglevel.size == len(wl)
+            if len(wl) == 3:
+                np.testing.assert_array_equal(ds.tas.isnull().all(), [True])
+            else:
+                np.testing.assert_array_equal(
+                    ds.tas.isnull().all(dim="time"), [False, True, True, True]
+                )
 
     def test_multireals(self):
         ds = self.ds.expand_dims(
@@ -456,12 +465,12 @@ def test_multireals(self):
         )
         ds_sub = xs.subset_warming_level(
             ds,
-            wl=1,
+            wl=1.5,
             to_level="tests",
         )
         np.testing.assert_array_equal(ds_sub.time.dt.year, np.arange(1000, 1020))
         np.testing.assert_array_equal(
-            ds_sub.warminglevel_bounds[:2].dt.year, [[[1990, 2009]], [[1990, 2009]]]
+            ds_sub.warminglevel_bounds[:2].dt.year, [[[2004, 2023]], [[2004, 2023]]]
         )
         assert ds_sub.warminglevel_bounds[2].isnull().all()
 
@@ -475,7 +484,7 @@ def test_multilevels(self):
             ds_sub.warminglevel,
             ["+1Cvs1850-1900", "+2Cvs1850-1900", "+3Cvs1850-1900", "+20Cvs1850-1900"],
         )
-        np.testing.assert_array_equal(ds_sub.tas.isnull().sum("time"), [10, 0, 1, 20])
+        np.testing.assert_array_equal(ds_sub.tas.isnull().sum("time"), [20, 0, 20, 20])
 
 
 class TestResample:

xscen/extract.py

@@ -983,7 +983,7 @@ def get_warming_level(  # noqa: C901
     -------
     dict, list or str
         If `realization` is not a sequence, the output will follow the format indicated by `return_horizon`.
-        If `realization` is a sequence, the output will be a list or dictionary depending on `output`,
+        If `realization` is a sequence, the output will be of the same type,
         with values following the format indicated by `return_horizon`.
     """
     tas_src = tas_src or Path(__file__).parent / "data" / "IPCC_annual_global_tas.nc"
@@ -1107,7 +1107,7 @@ def _get_warming_level(model):
     if isinstance(realization, xr.DataArray):
         if return_horizon:
             return xr.DataArray(
-                out, dims=(realization.dims[0], "bounds"), coords=realization.coords
+                out, dims=(realization.dims[0], "wl_bounds"), coords=realization.coords
             )
         return xr.DataArray(out, dims=(realization.dims[0],), coords=realization.coords)
 
@@ -1127,6 +1127,7 @@ def subset_warming_level(
     r"""
     Subsets the input dataset with only the window of time over which the requested level of global warming
     is first reached, using the IPCC Atlas method.
+    A warming level is considered reached only if the full `window` years are available in the dataset.
 
     Parameters
     ----------
@@ -1235,23 +1236,30 @@ def subset_warming_level(
         reals = []
         for real in bounds.realization.values:
             start, end = bounds.sel(realization=real).values
-            if start is not None:
-                data = ds.sel(realization=[real], time=slice(start, end))
+            data = ds.sel(realization=[real], time=slice(start, end))
+            wl_not_reached = (
+                (start is None)
+                or (data.time.size == 0)
+                or ((data.time.dt.year[-1] - data.time.dt.year[0] + 1) != window)
+            )
+            if not wl_not_reached:
                 bnds_crd = [
                     date_cls(int(start), 1, 1),
                     date_cls(int(end) + 1, 1, 1) - datetime.timedelta(seconds=1),
                 ]
             else:
+                # In the case of not reaching the WL, data might be too short
+                # We create it again with the proper length
                 data = (
                     ds.sel(realization=[real]).isel(time=slice(0, fake_time.size))
-                    * np.NaN
+                    * np.nan
                 )
-                bnds_crd = [np.NaN, np.NaN]
+                bnds_crd = [np.nan, np.nan]
             reals.append(
                 data.expand_dims(warminglevel=wl_crd).assign_coords(
                     time=fake_time[: data.time.size],
                     warminglevel_bounds=(
-                        ("realization", "warminglevel", "bounds"),
+                        ("realization", "warminglevel", "wl_bounds"),
                         [[bnds_crd]],
                     ),
                 )
@@ -1262,19 +1270,24 @@ def subset_warming_level(
         start_yr, end_yr = get_warming_level(ds, wl=wl, return_horizon=True, **kwargs)
         # cut the window selected above and expand dims with wl_crd
         ds_wl = ds.sel(time=slice(start_yr, end_yr))
+        wl_not_reached = (
+            (start_yr is None)
+            or (ds_wl.time.size == 0)
+            or ((ds_wl.time.dt.year[-1] - ds_wl.time.dt.year[0] + 1) != window)
+        )
         if fake_time is None:
-            # WL not reached or completely outside ds time
-            if start_yr is None or ds_wl.time.size == 0:
+            # WL not reached, not in ds, or not fully contained in ds.time
+            if wl_not_reached:
                 return None
             ds_wl = ds_wl.expand_dims(warminglevel=wl_crd)
         else:
-            # WL not reached or not completely inside ds time
-            if start_yr is None or ds_wl.time.size == 0:
+            # WL not reached, not in ds, or not fully contained in ds.time
+            if wl_not_reached:
                 ds_wl = ds.isel(time=slice(0, fake_time.size)) * np.NaN
-                wlbnds = (("warminglevel", "bounds"), [[np.NaN, np.NaN]])
+                wlbnds = (("warminglevel", "wl_bounds"), [[np.NaN, np.NaN]])
             else:
                 wlbnds = (
-                    ("warminglevel", "bounds"),
+                    ("warminglevel", "wl_bounds"),
                     [
                         [
                             date_cls(int(start_yr), 1, 1),

xscen/utils.py

@@ -1082,7 +1082,7 @@ def unstack_dates(
         )
 
     # Fast track for annual
-    if base == "A":
+    if base in "YA":
         if seasons:
             seaname = seasons[first.month]
         elif anchor == "JAN":