diff --git a/tests/conftest.py b/tests/conftest.py index 7aeeff089..846da7439 100644 --- a/tests/conftest.py +++ b/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"): diff --git a/tests/test_indices.py b/tests/test_indices.py index d48b249a7..026f1323d 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -257,6 +257,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): diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py index ef34fe6c3..04473704b 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 b7235bc75..e24d1aba4 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