Merge pull request #258 from Ouranosinc/fix-79

Accept float in days_over_precip_thresh

tests/conftest.py

@@ -82,6 +82,28 @@ def _pr_series(values, start="7/1/2000"):
     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"):

tests/test_indices.py

@@ -277,6 +277,30 @@ def test_simple(self, pr_series, per_doy):
             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):

xclim/indices/_multivariate.py

@@ -92,19 +92,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
 
@@ -594,7 +585,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
@@ -605,30 +596,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")
@@ -646,9 +627,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
@@ -659,28 +640,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
@@ -715,20 +687,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")
@@ -763,20 +727,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")
@@ -811,18 +767,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")
@@ -858,19 +808,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")
@@ -905,20 +848,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")
@@ -953,20 +888,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")
@@ -1062,18 +989,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
 

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
@@ -473,7 +472,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.
 
@@ -492,6 +491,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