From d51ca372fdb5ae2102b2b2e5890850fa8161cca7 Mon Sep 17 00:00:00 2001 From: David Huard Date: Wed, 21 Aug 2019 14:39:44 -0400 Subject: [PATCH 1/5] create resample_doy function. Accept float in days_over_precip_thresh --- tests/test_indices.py | 10 ++++++++++ xclim/indices/_multivariate.py | 19 ++++++------------ xclim/utils.py | 35 ++++++++++++++++++++++++++++++++-- 3 files changed, 49 insertions(+), 15 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index d48b249a7..5359235fd 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -257,6 +257,16 @@ def test_simple(self, pr_series, per_doy): out[0], (3 + 4 + 6 + 7) / (3 + 4 + 5 + 6 + 7) ) + def test_float(self, pr_series): + a = np.zeros(365) + a[:8] = np.arange(8) + pr = pr_series(a, start="1/1/2000") + + out = xci.days_over_precip_thresh(pr, "5 kg/m**2/s", thresh="2 kg/m**2/s") + np.testing.assert_array_almost_equal( + out[0], 2 + ) # Only days 6 and 7 meet criteria. + class TestFreshetStart: def test_simple(self, tas_series): diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py index ef34fe6c3..5f678d078 100644 --- a/xclim/indices/_multivariate.py +++ b/xclim/indices/_multivariate.py @@ -567,8 +567,8 @@ 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] + per : xarray.DataArray or string (scalar with units) + 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 @@ -588,21 +588,14 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"): >>> p75 = percentile_doy(historical_pr, per=0.75) >>> 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 isinstance(per, xr.DataArray): + 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") diff --git a/xclim/utils.py b/xclim/utils.py index b7235bc75..f7c368676 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.keys(): + 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 From 37a58990d48cd5c2a6501d7e2ecc43ee2e0fb4f2 Mon Sep 17 00:00:00 2001 From: David Huard Date: Wed, 21 Aug 2019 15:50:52 -0400 Subject: [PATCH 2/5] Update xclim/indices/_multivariate.py Co-Authored-By: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> --- xclim/indices/_multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py index 5f678d078..fbe71b0a1 100644 --- a/xclim/indices/_multivariate.py +++ b/xclim/indices/_multivariate.py @@ -567,7 +567,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 or string (scalar with units) + per : xarray.DataArray or str 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]. From a8ab6f9980e3453d5a9db321ca2b17ff8ab0746b Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 22 Aug 2019 09:02:42 -0400 Subject: [PATCH 3/5] used resample_doy in multivariate indices --- xclim/indices/_multivariate.py | 127 +++++++++------------------------ 1 file changed, 32 insertions(+), 95 deletions(-) diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py index 5f678d078..3bc55cccd 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 @@ -567,7 +558,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 or string (scalar with units) + per : xarray.DataArray or string 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]. @@ -593,6 +584,7 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"): tp = np.maximum(per, thresh) if isinstance(per, xr.DataArray): + # Create time series out of doy values. tp = utils.resample_doy(tp, pr) # Compute the days where precip is both over the wet day threshold and the percentile threshold. @@ -613,9 +605,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] - per : xarray.DataArray - Daily percentile of wet day 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 or str + 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 @@ -636,18 +628,15 @@ def fraction_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"): 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 isinstance(per, xr.DataArray): + # 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 @@ -682,20 +671,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 time series out of doy values. + thresh = utils.resample_doy(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) - - # compute the cold days + # Identify the days over the 90th percentile over = tas > thresh return over.resample(time=freq).sum(dim="time") @@ -730,20 +711,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 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) + # Create time series out of doy values. + thresh = utils.resample_doy(t10, tas) - # compute the cold days + # Identify the days below the 10th percentile below = tas < thresh return below.resample(time=freq).sum(dim="time") @@ -778,18 +751,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") @@ -825,19 +792,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 time series out of doy values. + thresh = utils.resample_doy(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) - - # compute the cold days + # Identify the days below the 10th percentile below = tasmin < thresh return below.resample(time=freq).sum(dim="time") @@ -872,20 +832,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") @@ -920,20 +872,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") @@ -1029,18 +973,11 @@ 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 From 6f55eded429104c047ec1a299333b82987ea3782 Mon Sep 17 00:00:00 2001 From: David Huard Date: Thu, 12 Sep 2019 16:06:00 -0400 Subject: [PATCH 4/5] removed support for string percentile treshold. --- tests/conftest.py | 22 ++++++++++++++++++++++ tests/test_indices.py | 18 ++++++++++++++++-- xclim/indices/_multivariate.py | 21 ++++----------------- xclim/utils.py | 2 +- 4 files changed, 43 insertions(+), 20 deletions(-) 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 5359235fd..026f1323d 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -257,16 +257,30 @@ def test_simple(self, pr_series, per_doy): out[0], (3 + 4 + 6 + 7) / (3 + 4 + 5 + 6 + 7) ) - def test_float(self, pr_series): + def test_quantile(self, pr_series): a = np.zeros(365) a[:8] = np.arange(8) pr = pr_series(a, start="1/1/2000") - out = xci.days_over_precip_thresh(pr, "5 kg/m**2/s", thresh="2 kg/m**2/s") + # 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 0de14d414..c8ba79246 100644 --- a/xclim/indices/_multivariate.py +++ b/xclim/indices/_multivariate.py @@ -558,7 +558,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 or str + per : xarray.DataArray 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]. @@ -570,10 +570,6 @@ 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) @@ -583,7 +579,7 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"): thresh = utils.convert_units_to(thresh, pr) tp = np.maximum(per, thresh) - if isinstance(per, xr.DataArray): + if "dayofyear" in per.coords: # Create time series out of doy values. tp = utils.resample_doy(tp, pr) @@ -606,7 +602,7 @@ def fraction_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 or str + per : xarray.DataArray 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]. @@ -618,18 +614,12 @@ 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) tp = np.maximum(per, thresh) - if isinstance(per, xr.DataArray): + if "dayofyear" in per.coords: # Create time series out of doy values. tp = utils.resample_doy(tp, pr) @@ -973,9 +963,6 @@ def warm_spell_duration_index(tasmax, tx90, window=6, freq="YS"): precipitation, J. Geophys. Res., 111, D05109, doi: 10.1029/2005JD006290. """ - # The day of year value of the tasmax series. - doy = tasmax.indexes["time"].dayofyear - # Create time series out of doy values. thresh = utils.resample_doy(tx90, tasmax) diff --git a/xclim/utils.py b/xclim/utils.py index f7c368676..e24d1aba4 100644 --- a/xclim/utils.py +++ b/xclim/utils.py @@ -497,7 +497,7 @@ def resample_doy(doy, arr): 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.keys(): + if "dayofyear" not in doy.coords: raise AttributeError("`doy` should have dayofyear coordinates.") # Adjust calendar From e7470094a65feea405d28a9de4079e48243fd493 Mon Sep 17 00:00:00 2001 From: David Huard Date: Fri, 13 Sep 2019 15:12:26 -0400 Subject: [PATCH 5/5] modified days_over_precip_thresh example to use quantile. --- xclim/indices/_multivariate.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py index c8ba79246..04473704b 100644 --- a/xclim/indices/_multivariate.py +++ b/xclim/indices/_multivariate.py @@ -572,7 +572,7 @@ def days_over_precip_thresh(pr, per, thresh="1 mm/day", freq="YS"): 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) """ per = utils.convert_units_to(per, pr)