Repository: Ouranosinc/xclim

Issue:
Percentile threshold indices : Frequencies?
For indices such as 'R95p' (with the following definition taken from icclim)
-  Days with RR > 95th percentile of daily amounts (very wet days)
(days)

- Let RRwj be the daily precipitation amount at wet day w (RR ≥ 1.0 mm)
of period j and let RRwn95 be the 95th percentile of precipitation at wet
days in the 1961–1990 period. Then counted is the number of days where:

- RRwj > RRwn95

Does the definition of RRwn95 (the percentile threshold value calculated from the reference period) change with the frequency? i.e if the R95p resampling value is seasonal does the RRwn95 value get calculated seasonally? or is the percentile threshold always annual?  ICCLIM docs are unclear on thisI

Discussion:

Patch:
diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py
index ef34fe6c..04473704 100644
--- a/xclim/indices/_multivariate.py
+++ b/xclim/indices/_multivariate.py
@@ -91,19 +91,10 @@ def cold_spell_duration_index(tasmin, tn10, window=6, freq="YS"):
     >>> tn10 = percentile_doy(historical_tasmin, per=.1)
     >>> cold_spell_duration_index(reference_tasmin, tn10)
     """
-    if "dayofyear" not in tn10.coords.keys():
-        raise AttributeError("tn10 should have dayofyear coordinates.")
-
-    # The day of year value of the tasmin series.
-    doy = tasmin.indexes["time"].dayofyear
-
     tn10 = utils.convert_units_to(tn10, tasmin)
-    # If calendar of `tn10` is different from `tasmin`, interpolate.
-    tn10 = utils.adjust_doy_calendar(tn10, tasmin)
 
-    # Create an array with the shape and coords of tasmin, but with values set to tx90 according to the doy index.
-    thresh = xr.full_like(tasmin, np.nan)
-    thresh.data = tn10.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(tn10, tasmin)
 
     below = tasmin < thresh
 
@@ -568,7 +559,7 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"):
     pr : xarray.DataArray
       Mean daily precipitation flux [Kg m-2 s-1] or [mm/day]
     per : xarray.DataArray
-      Daily percentile of wet day precipitation flux [Kg m-2 s-1] or [mm/day]
+      Daily percentile of wet day precipitation flux [Kg m-2 s-1] or [mm/day].
     thresh : str
        Precipitation value over which a day is considered wet [Kg m-2 s-1] or [mm/day].
     freq : str, optional
@@ -579,30 +570,20 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"):
     xarray.DataArray
       Count of days with daily precipitation above the given percentile [days]
 
-    Notes
-    -----
-    The percentile should be computed for a 5 day window centered on each calendar day for a reference period.
-
     Example
     -------
-    >>> p75 = percentile_doy(historical_pr, per=0.75)
+    >>> p75 = historical_pr.quantile(.75, dim="time", keep_attrs=True)
     >>> r75p = days_over_precip_thresh(pr, p75)
     """
-    if "dayofyear" not in per.coords.keys():
-        raise AttributeError("percentile should have dayofyear coordinates.")
-
     per = utils.convert_units_to(per, pr)
     thresh = utils.convert_units_to(thresh, pr)
 
-    per = utils.adjust_doy_calendar(per, pr)
-    mper = np.maximum(per, thresh)
+    tp = np.maximum(per, thresh)
+    if "dayofyear" in per.coords:
+        # Create time series out of doy values.
+        tp = utils.resample_doy(tp, pr)
 
-    # create array of percentile with pr shape and coords
-    tp = xr.full_like(pr, np.nan)
-    doy = tp.time.dt.dayofyear.values
-    tp.data = mper.sel(dayofyear=doy)
-
-    # compute the days where precip is both over the wet day threshold and the percentile threshold.
+    # Compute the days where precip is both over the wet day threshold and the percentile threshold.
     over = pr > tp
 
     return over.resample(time=freq).sum(dim="time")
@@ -620,9 +601,9 @@ def fraction_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"):
     Parameters
     ----------
     pr : xarray.DataArray
-      Mean daily precipitation flux [Kg m-2 s-1] or [mm/day]
+      Mean daily precipitation flux [Kg m-2 s-1] or [mm/day].
     per : xarray.DataArray
-      Daily percentile of wet day precipitation flux [Kg m-2 s-1] or [mm/day]
+      Daily percentile of wet day precipitation flux [Kg m-2 s-1] or [mm/day].
     thresh : str
        Precipitation value over which a day is considered wet [Kg m-2 s-1] or [mm/day].
     freq : str, optional
@@ -633,28 +614,19 @@ def fraction_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"):
     xarray.DataArray
       Fraction of precipitation over threshold during wet days days.
 
-    Notes
-    -----
-    The percentile should be computed for a 5 day window centered on each calendar day for a reference period.
     """
-    if "dayofyear" not in per.coords.keys():
-        raise AttributeError("percentile should have dayofyear coordinates.")
-
     per = utils.convert_units_to(per, pr)
     thresh = utils.convert_units_to(thresh, pr)
 
-    per = utils.adjust_doy_calendar(per, pr)
-    mper = np.maximum(per, thresh)
-
-    # create array of percentile with pr shape and coords
-    tp = xr.full_like(pr, np.nan)
-    doy = tp.time.dt.dayofyear.values
-    tp.data = mper.sel(dayofyear=doy)
+    tp = np.maximum(per, thresh)
+    if "dayofyear" in per.coords:
+        # Create time series out of doy values.
+        tp = utils.resample_doy(tp, pr)
 
     # Total precip during wet days over period
     total = pr.where(pr > thresh).resample(time=freq).sum(dim="time")
 
-    # compute the days where precip is both over the wet day threshold and the percentile threshold.
+    # Compute the days where precip is both over the wet day threshold and the percentile threshold.
     over = pr.where(pr > tp).resample(time=freq).sum(dim="time")
 
     return over / total
@@ -689,20 +661,12 @@ def tg90p(tas, t90, freq="YS"):
     >>> t90 = percentile_doy(historical_tas, per=0.9)
     >>> hot_days = tg90p(tas, t90)
     """
-    if "dayofyear" not in t90.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
-
     t90 = utils.convert_units_to(t90, tas)
 
-    # adjustment of t90 to tas doy range
-    t90 = utils.adjust_doy_calendar(t90, tas)
-
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tas, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t90.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t90, tas)
 
-    # compute the cold days
+    # Identify the days over the 90th percentile
     over = tas > thresh
 
     return over.resample(time=freq).sum(dim="time")
@@ -737,20 +701,12 @@ def tg10p(tas, t10, freq="YS"):
     >>> t10 = percentile_doy(historical_tas, per=0.1)
     >>> cold_days = tg10p(tas, t10)
     """
-    if "dayofyear" not in t10.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
-
     t10 = utils.convert_units_to(t10, tas)
 
-    # adjustment of t10 to tas doy range
-    t10 = utils.adjust_doy_calendar(t10, tas)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t10, tas)
 
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tas, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t10.sel(dayofyear=doy)
-
-    # compute the cold days
+    # Identify the days below the 10th percentile
     below = tas < thresh
 
     return below.resample(time=freq).sum(dim="time")
@@ -785,18 +741,12 @@ def tn90p(tasmin, t90, freq="YS"):
     >>> t90 = percentile_doy(historical_tas, per=0.9)
     >>> hot_days = tg90p(tas, t90)
     """
-    if "dayofyear" not in t90.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
     t90 = utils.convert_units_to(t90, tasmin)
-    # adjustment of t90 to tas doy range
-    t90 = utils.adjust_doy_calendar(t90, tasmin)
 
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tasmin, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t90.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t90, tasmin)
 
-    # compute the cold days
+    # Identify the days with min temp above 90th percentile.
     over = tasmin > thresh
 
     return over.resample(time=freq).sum(dim="time")
@@ -832,19 +782,12 @@ def tn10p(tasmin, t10, freq="YS"):
     >>> t10 = percentile_doy(historical_tas, per=0.1)
     >>> cold_days = tg10p(tas, t10)
     """
-    if "dayofyear" not in t10.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
     t10 = utils.convert_units_to(t10, tasmin)
 
-    # adjustment of t10 to tas doy range
-    t10 = utils.adjust_doy_calendar(t10, tasmin)
-
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tasmin, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t10.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t10, tasmin)
 
-    # compute the cold days
+    # Identify the days below the 10th percentile
     below = tasmin < thresh
 
     return below.resample(time=freq).sum(dim="time")
@@ -879,20 +822,12 @@ def tx90p(tasmax, t90, freq="YS"):
     >>> t90 = percentile_doy(historical_tas, per=0.9)
     >>> hot_days = tg90p(tas, t90)
     """
-    if "dayofyear" not in t90.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
-
     t90 = utils.convert_units_to(t90, tasmax)
 
-    # adjustment of t90 to tas doy range
-    t90 = utils.adjust_doy_calendar(t90, tasmax)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t90, tasmax)
 
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tasmax, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t90.sel(dayofyear=doy)
-
-    # compute the cold days
+    # Identify the days with max temp above 90th percentile.
     over = tasmax > thresh
 
     return over.resample(time=freq).sum(dim="time")
@@ -927,20 +862,12 @@ def tx10p(tasmax, t10, freq="YS"):
     >>> t10 = percentile_doy(historical_tas, per=0.1)
     >>> cold_days = tg10p(tas, t10)
     """
-    if "dayofyear" not in t10.coords.keys():
-        raise AttributeError("t10 should have dayofyear coordinates.")
-
     t10 = utils.convert_units_to(t10, tasmax)
 
-    # adjustment of t10 to tas doy range
-    t10 = utils.adjust_doy_calendar(t10, tasmax)
-
-    # create array of percentile with tas shape and coords
-    thresh = xr.full_like(tasmax, np.nan)
-    doy = thresh.time.dt.dayofyear.values
-    thresh.data = t10.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(t10, tasmax)
 
-    # compute the cold days
+    # Identify the days below the 10th percentile
     below = tasmax < thresh
 
     return below.resample(time=freq).sum(dim="time")
@@ -1036,18 +963,8 @@ def warm_spell_duration_index(tasmax, tx90, window=6, freq="YS"):
     precipitation, J. Geophys. Res., 111, D05109, doi: 10.1029/2005JD006290.
 
     """
-    if "dayofyear" not in tx90.coords.keys():
-        raise AttributeError("tx90 should have dayofyear coordinates.")
-
-    # The day of year value of the tasmax series.
-    doy = tasmax.indexes["time"].dayofyear
-
-    # adjustment of tx90 to tasmax doy range
-    tx90 = utils.adjust_doy_calendar(tx90, tasmax)
-
-    # Create an array with the shape and coords of tasmax, but with values set to tx90 according to the doy index.
-    thresh = xr.full_like(tasmax, np.nan)
-    thresh.data = tx90.sel(dayofyear=doy)
+    # Create time series out of doy values.
+    thresh = utils.resample_doy(tx90, tasmax)
 
     above = tasmax > thresh
 
diff --git a/xclim/utils.py b/xclim/utils.py
index b7235bc7..e24d1aba 100644
--- a/xclim/utils.py
+++ b/xclim/utils.py
@@ -197,7 +197,6 @@ def convert_units_to(source, target, context=None):
       Target array of values to which units must conform.
     context : str
 
-
     Returns
     -------
     out
@@ -463,7 +462,7 @@ def adjust_doy_calendar(source, target):
     Parameters
     ----------
     source : xarray.DataArray
-      Array with `dayofyear` coordinates.
+      Array with `dayofyear` coordinate.
     target : xarray.DataArray
       Array with `time` coordinate.
 
@@ -482,6 +481,38 @@ def adjust_doy_calendar(source, target):
     return _interpolate_doy_calendar(source, doy_max)
 
 
+def resample_doy(doy, arr):
+    """Create a temporal DataArray where each day takes the value defined by the day-of-year.
+
+    Parameters
+    ----------
+    doy : xarray.DataArray
+      Array with `dayofyear` coordinate.
+    arr : xarray.DataArray
+      Array with `time` coordinate.
+
+    Returns
+    -------
+    xarray.DataArray
+      An array with the same `time` dimension as `arr` whose values are filled according to the day-of-year value in
+      `doy`.
+    """
+    if "dayofyear" not in doy.coords:
+        raise AttributeError("`doy` should have dayofyear coordinates.")
+
+    # Adjust calendar
+    adoy = adjust_doy_calendar(doy, arr)
+
+    # Create array with arr shape and coords
+    out = xr.full_like(arr, np.nan)
+
+    # Fill with values from `doy`
+    d = out.time.dt.dayofyear.values
+    out.data = adoy.sel(dayofyear=d)
+
+    return out
+
+
 def get_daily_events(da, da_value, operator):
     r"""
     function that returns a 0/1 mask when a condition is True or False


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
index 7aeeff08..846da743 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -82,6 +82,28 @@ def pr_series():
     return _pr_series
 
 
+@pytest.fixture
+def pr_ndseries():
+    def _pr_series(values, start="1/1/2000"):
+        nt, nx, ny = np.atleast_3d(values).shape
+        time = pd.date_range(start, periods=nt, freq=pd.DateOffset(days=1))
+        x = np.arange(nx)
+        y = np.arange(ny)
+        return xr.DataArray(
+            values,
+            coords=[time, x, y],
+            dims=("time", "x", "y"),
+            name="pr",
+            attrs={
+                "standard_name": "precipitation_flux",
+                "cell_methods": "time: sum over day",
+                "units": "kg m-2 s-1",
+            },
+        )
+
+    return _pr_series
+
+
 @pytest.fixture
 def q_series():
     def _q_series(values, start="1/1/2000"):
diff --git a/tests/test_indices.py b/tests/test_indices.py
index d48b249a..026f1323 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -257,6 +257,30 @@ class TestDaysOverPrecipThresh:
             out[0], (3 + 4 + 6 + 7) / (3 + 4 + 5 + 6 + 7)
         )
 
+    def test_quantile(self, pr_series):
+        a = np.zeros(365)
+        a[:8] = np.arange(8)
+        pr = pr_series(a, start="1/1/2000")
+
+        # Create synthetic percentile
+        pr0 = pr_series(np.ones(365) * 5, start="1/1/2000")
+        per = pr0.quantile(0.5, dim="time", keep_attrs=True)
+        per.attrs["units"] = "kg m-2 s-1"  # This won't be needed with xarray 0.13
+
+        out = xci.days_over_precip_thresh(pr, per, thresh="2 kg/m**2/s")
+        np.testing.assert_array_almost_equal(
+            out[0], 2
+        )  # Only days 6 and 7 meet criteria.
+
+    def test_nd(self, pr_ndseries):
+        pr = pr_ndseries(np.ones((300, 2, 3)))
+        pr0 = pr_ndseries(np.zeros((300, 2, 3)))
+        per = pr0.quantile(0.5, dim="time", keep_attrs=True)
+        per.attrs["units"] = "kg m-2 s-1"  # This won't be needed with xarray 0.13
+
+        out = xci.days_over_precip_thresh(pr, per, thresh="0.5 kg/m**2/s")
+        np.testing.assert_array_almost_equal(out, 300)
+
 
 class TestFreshetStart:
     def test_simple(self, tas_series):