Repository: Ouranosinc/xclim

Issue:
Implement a 'doy_to_days_after' and reverse functionality

### Description
Certain indices outputting 'day of year' or `doy` values are not directly usable for ensemble percentiles and statistics due to its circular nature.   

Averaging `doy` results can potentially be handled by users via https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.stats.circmean.html.  However I have not easily found a percentiles option.

The other solution is to  create offset functions allowing units to become a 'days after %mm%dd' or similar which will allow valid ensemble statistics to be calculated..   Final results (mean, percentiles etc) could potentially be restored to original `doy` via a reverse function.

To streamline wps calls in `finch` it would be potentially interesting to have an `units` parameter for the calculation of indicators that currently return to `doy` to output `days_after` instead.  Non exhaustive list of indices include  `degree_day_exceedance`, `first_day_*`, `last_day_*`  

Note I think xarray may interpret 'days_since' as a cftime flag and try to automatically 'help' in certain cases. This is why I suggest a `days_after` unit for now but I am open to debate / validation

 

 





```


Discussion:

Patch:
diff --git a/xclim/core/calendar.py b/xclim/core/calendar.py
index e2936660..86a3787f 100644
--- a/xclim/core/calendar.py
+++ b/xclim/core/calendar.py
@@ -872,7 +872,7 @@ def doy_to_days_since(
       Passing `start` only makes sense if `da` has a yearly sampling frequency.
     calendar: str, optional
       The calendar to use when computing the new interval.
-      If None (default), the calendar of the `time` axis is used.
+      If None (default), the calendar attribute of the data or of its `time` axis is used.
       All time coordinates of `da` must exist in this calendar.
       No check is done to ensure doy values exist in this calendar.
 
@@ -897,7 +897,7 @@ def doy_to_days_since(
     array([-86, 92])
     """
     base_calendar = get_calendar(da)
-    calendar = calendar or base_calendar
+    calendar = calendar or da.attrs.get("calendar", base_calendar)
     dac = convert_calendar(da, calendar)
 
     base_doy, start_doy, doy_max = _doy_days_since_doys(dac.time, start)
@@ -916,6 +916,7 @@ def doy_to_days_since(
         else:
             out.attrs.update(units="days after time coordinate")
 
+    out.attrs.pop("is_dayofyear", None)
     out.attrs.update(calendar=calendar)
     return convert_calendar(out, base_calendar).rename(da.name)
 
@@ -975,4 +976,5 @@ def days_since_to_doy(
     out.attrs.update(
         {k: v for k, v in da.attrs.items() if k not in ["units", "calendar"]}
     )
+    out.attrs.update(calendar=calendar, is_dayofyear=True)
     return convert_calendar(out, base_calendar).rename(da.name)
diff --git a/xclim/indices/_hydrology.py b/xclim/indices/_hydrology.py
index c15952f7..f5088642 100644
--- a/xclim/indices/_hydrology.py
+++ b/xclim/indices/_hydrology.py
@@ -1,6 +1,7 @@
 import numpy as np
 import xarray
 
+from xclim.core.calendar import get_calendar
 from xclim.core.units import declare_units, rate2amount, to_agg_units
 
 from . import generic
@@ -122,7 +123,7 @@ def snd_max_doy(snd: xarray.DataArray, freq: str = "AS-JUL") -> xarray.DataArray
 
     # Compute doymax. Will return first time step if all snow depths are 0.
     out = generic.select_resample_op(snd, op=generic.doymax, freq=freq)
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(snd))
 
     # Mask arrays that miss at least one non-null snd.
     return out.where(~valid)
diff --git a/xclim/indices/_threshold.py b/xclim/indices/_threshold.py
index 6afe6c39..030fabb0 100644
--- a/xclim/indices/_threshold.py
+++ b/xclim/indices/_threshold.py
@@ -4,6 +4,7 @@ from typing import Optional
 import numpy as np
 import xarray
 
+from xclim.core.calendar import get_calendar
 from xclim.core.units import (
     convert_units_to,
     declare_units,
@@ -192,7 +193,7 @@ def continuous_snow_cover_end(
         .map(rl.season, window=window, dim="time", coord="dayofyear")
         .end
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(snd))
     return out
 
 
@@ -239,7 +240,7 @@ def continuous_snow_cover_start(
         )
         .start
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(snd))
     return out
 
 
@@ -466,7 +467,7 @@ def freshet_start(
     thresh = convert_units_to(thresh, tas)
     over = tas > thresh
     out = over.resample(time=freq).map(rl.first_run, window=window, coord="dayofyear")
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tas))
     return out
 
 
@@ -549,7 +550,7 @@ def growing_season_end(
         dim="time",
         coord="dayofyear",
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tas))
     return out
 
 
@@ -742,7 +743,7 @@ def last_spring_frost(
         dim="time",
         coord="dayofyear",
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tas))
     return out
 
 
@@ -790,7 +791,7 @@ def first_day_below(
         dim="time",
         coord="dayofyear",
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tasmin))
     return out
 
 
@@ -838,7 +839,7 @@ def first_day_above(
         dim="time",
         coord="dayofyear",
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tasmin))
     return out
 
 
@@ -882,7 +883,7 @@ def first_snowfall(
         dim="time",
         coord="dayofyear",
     )
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(prsn))
     return out
 
 
@@ -1720,7 +1721,7 @@ def degree_days_exceedance_date(
         )
 
     out = c.clip(0).resample(time=freq).map(_exceedance_date)
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(tas))
     return out
 
 
diff --git a/xclim/indices/generic.py b/xclim/indices/generic.py
index 6f4d5473..78aafb82 100644
--- a/xclim/indices/generic.py
+++ b/xclim/indices/generic.py
@@ -11,6 +11,8 @@ from typing import Union
 import numpy as np
 import xarray as xr
 
+from xclim.core.calendar import get_calendar
+
 __all__ = [
     "select_time",
     "select_resample_op",
@@ -84,7 +86,7 @@ def doymax(da: xr.DataArray):
     """Return the day of year of the maximum value."""
     i = da.argmax(dim="time")
     out = da.time.dt.dayofyear[i]
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(da))
     return out
 
 
@@ -92,7 +94,7 @@ def doymin(da: xr.DataArray):
     """Return the day of year of the minimum value."""
     i = da.argmin(dim="time")
     out = da.time.dt.dayofyear[i]
-    out.attrs["units"] = ""
+    out.attrs.update(units="", is_dayofyear=True, calendar=get_calendar(da))
     return out
 
 


Test patch:
diff --git a/xclim/testing/tests/test_calendar.py b/xclim/testing/tests/test_calendar.py
index eed8ccce..e263f1ef 100644
--- a/xclim/testing/tests/test_calendar.py
+++ b/xclim/testing/tests/test_calendar.py
@@ -384,12 +384,18 @@ def test_datetime_to_decimal_year(source_cal, exp180):
 def test_doy_to_days_since():
     # simple test
     time = date_range("2020-07-01", "2022-07-01", freq="AS-JUL")
-    da = xr.DataArray([190, 360, 3], dims=("time",), coords={"time": time})
+    da = xr.DataArray(
+        [190, 360, 3],
+        dims=("time",),
+        coords={"time": time},
+        attrs={"is_dayofyear": True, "calendar": "default"},
+    )
 
     out = doy_to_days_since(da)
     np.testing.assert_array_equal(out, [7, 178, 186])
 
     assert out.attrs["units"] == "days after 07-01"
+    assert "is_dayofyear" not in out.attrs
 
     da2 = days_since_to_doy(out)
     xr.testing.assert_identical(da, da2)
@@ -403,11 +409,19 @@ def test_doy_to_days_since():
     np.testing.assert_array_equal(out, [8, 178, 186])
 
     da2 = days_since_to_doy(out)  # calendar read from attribute
+    da2.attrs.pop("calendar")  # drop for identicality
+    da.attrs.pop("calendar")  # drop for identicality
     xr.testing.assert_identical(da, da2)
 
     # with start
     time = date_range("2020-12-31", "2022-12-31", freq="Y")
-    da = xr.DataArray([190, 360, 3], dims=("time",), coords={"time": time}, name="da")
+    da = xr.DataArray(
+        [190, 360, 3],
+        dims=("time",),
+        coords={"time": time},
+        name="da",
+        attrs={"is_dayofyear": True, "calendar": "default"},
+    )
 
     out = doy_to_days_since(da, start="01-02")
     np.testing.assert_array_equal(out, [188, 358, 1])
@@ -418,7 +432,13 @@ def test_doy_to_days_since():
 
     # finer freq
     time = date_range("2020-01-01", "2020-03-01", freq="MS")
-    da = xr.DataArray([15, 33, 66], dims=("time",), coords={"time": time})
+    da = xr.DataArray(
+        [15, 33, 66],
+        dims=("time",),
+        coords={"time": time},
+        name="da",
+        attrs={"is_dayofyear": True, "calendar": "default"},
+    )
 
     out = doy_to_days_since(da)
     assert out.attrs["units"] == "days after time coordinate"
diff --git a/xclim/testing/tests/test_fwi.py b/xclim/testing/tests/test_fwi.py
index ae9679b2..737db50b 100644
--- a/xclim/testing/tests/test_fwi.py
+++ b/xclim/testing/tests/test_fwi.py
@@ -170,7 +170,7 @@ def test_fire_weather_ufunc_overwintering():
         pr=convert_units_to(ds.pr, "mm/d"),
     )
     season_mask_all = fire_season(ds.tas, method="WF93", temp_end_thresh="4 degC")
-    season_mask_all_LA08 = fire_season(ds.tas, snd=ds.snd, method="LA08")
+    season_mask_all_LA08 = fire_season(ds.tas, snd=ds.swe, method="LA08")
     season_mask_yr = fire_season(ds.tas, method="WF93", freq="YS")
 
     # Mask is computed correctly and parameters are passed
@@ -189,7 +189,7 @@ def test_fire_weather_ufunc_overwintering():
     out2 = fire_weather_ufunc(
         tas=ds.tas,
         pr=ds.pr,
-        snd=ds.snd,
+        snd=ds.swe,
         lat=ds.lat,
         season_method="LA08",
         overwintering=True,
diff --git a/xclim/testing/tests/test_generic.py b/xclim/testing/tests/test_generic.py
index 8fb03e13..413f692f 100644
--- a/xclim/testing/tests/test_generic.py
+++ b/xclim/testing/tests/test_generic.py
@@ -137,4 +137,8 @@ def test_doyminmax(q_series):
     dmn = generic.doymin(q)
     assert dmx.values == [40]
     assert dmn.values == [50]
-    assert dmx.units == ""
+    for da in [dmx, dmn]:
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in da.attrs.keys()
+        assert da.attrs["units"] == ""
+        assert da.attrs["is_dayofyear"]
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index 36cad653..87ea550a 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -282,6 +282,10 @@ class TestLastSpringFrost:
 
         lsf = xci.last_spring_frost(tas)
         assert lsf == 180
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in lsf.attrs.keys()
+        assert lsf.attrs["units"] == ""
+        assert lsf.attrs["is_dayofyear"]
 
 
 class TestFirstDayBelow:
@@ -298,6 +302,10 @@ class TestFirstDayBelow:
 
         fdb = xci.first_day_below(tas)
         assert np.isnan(fdb)
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in fdb.attrs.keys()
+        assert fdb.attrs["units"] == ""
+        assert fdb.attrs["is_dayofyear"]
 
 
 class TestFirstDayAbove:
@@ -317,6 +325,10 @@ class TestFirstDayAbove:
 
         fdb = xci.first_day_above(tas)
         assert np.isnan(fdb)
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in fdb.attrs.keys()
+        assert fdb.attrs["units"] == ""
+        assert fdb.attrs["is_dayofyear"]
 
 
 class TestDaysOverPrecipThresh:
@@ -378,6 +390,10 @@ class TestFreshetStart:
         tg = tas_series(tg + K2C, start="1/1/2000")
         out = xci.freshet_start(tg, window=w)
         assert out[0] == tg.indexes["time"][20].dayofyear
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in out.attrs.keys()
+        assert out.attrs["units"] == ""
+        assert out.attrs["is_dayofyear"]
 
     def test_no_start(self, tas_series):
         tg = np.zeros(365) - 1
@@ -414,6 +430,10 @@ class TestGrowingSeasonEnd:
         tas = tas.where(~tas.time.isin(warm_period.time), 280)
         gs_end = xci.growing_season_end(tas, mid_date=mid_date)
         np.testing.assert_array_equal(gs_end, expected)
+        for attr in ["units", "is_dayofyear", "calendar"]:
+            assert attr in gs_end.attrs.keys()
+        assert gs_end.attrs["units"] == ""
+        assert gs_end.attrs["is_dayofyear"]
 
 
 class TestGrowingSeasonLength:
@@ -1680,6 +1700,11 @@ def test_degree_days_exceedance_date(tas_series):
     )
     assert out[0] == 256
 
+    for attr in ["units", "is_dayofyear", "calendar"]:
+        assert attr in out.attrs.keys()
+    assert out.attrs["units"] == ""
+    assert out.attrs["is_dayofyear"]
+
 
 @pytest.mark.parametrize("method,exp", [("binary", [1, 1, 1, 1, 1, 0, 0, 0, 0, 0])])
 def test_snowfall_approximation(pr_series, tasmax_series, method, exp):
@@ -1705,6 +1730,10 @@ def test_first_snowfall(prsn_series):
     prsn = prsn_series(30 - abs(np.arange(366) - 180), start="2000-01-01")
     out = xci.first_snowfall(prsn, thresh="15 kg m-2 s-1", freq="YS")
     assert out[0] == 166
+    for attr in ["units", "is_dayofyear", "calendar"]:
+        assert attr in out.attrs.keys()
+    assert out.attrs["units"] == ""
+    assert out.attrs["is_dayofyear"]
 
 
 def test_last_snowfall(prsn_series):
@@ -1739,6 +1768,10 @@ def test_continous_snow_cover_start(snd_series):
     out = xci.continuous_snow_cover_start(snd)
     assert len(out) == 2
     np.testing.assert_array_equal(out, [snd.time.dt.dayofyear[0].data + 2, np.nan])
+    for attr in ["units", "is_dayofyear", "calendar"]:
+        assert attr in out.attrs.keys()
+    assert out.attrs["units"] == ""
+    assert out.attrs["is_dayofyear"]
 
 
 def test_continuous_snow_cover_end(snd_series):
@@ -1756,6 +1789,10 @@ def test_continuous_snow_cover_end(snd_series):
     assert len(out) == 2
     doy = snd.time.dt.dayofyear[0].data
     np.testing.assert_array_equal(out, [(doy + 219) % 366, np.nan])
+    for attr in ["units", "is_dayofyear", "calendar"]:
+        assert attr in out.attrs.keys()
+    assert out.attrs["units"] == ""
+    assert out.attrs["is_dayofyear"]
 
 
 def test_high_precip_low_temp(pr_series, tasmin_series):
diff --git a/xclim/testing/tests/test_run_length.py b/xclim/testing/tests/test_run_length.py
index eac4767b..e7ed252f 100644
--- a/xclim/testing/tests/test_run_length.py
+++ b/xclim/testing/tests/test_run_length.py
@@ -283,7 +283,7 @@ def test_keep_longest_run_synthetic():
 
 def test_keep_longest_run_data():
     era5 = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc")
-    cond = era5.snd > 0.001
+    cond = era5.swe > 0.001
     lrun = rl.keep_longest_run(cond, "time")
     np.testing.assert_array_equal(
         lrun.isel(time=slice(651, 658), location=2),