From d5f3436b8bcb03835927e342e2bed8cd282a5164 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 11 Oct 2024 13:14:38 -0400 Subject: [PATCH 1/7] add support for freq="W" --- CHANGELOG.rst | 1 + xclim/indices/stats.py | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 4368d61a5..f9cc4583a 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,6 +23,7 @@ New features and enhancements * Indicator parameters can now be assigned a new name, different from the argument name in the compute function. (:pull:`1885`). * ``xclim.indices.run_length.windowed_max_run_sum`` accumulates positive values across runs and yields the the maximum valued run. (:pull:`1926`). * Helper function ``xclim.indices.helpers.make_hourly_temperature`` to estimate hourly temperatures from daily min and max temperatures. (:pull:`1909`). +* `xclim.indices.stats.standardized_index` now supports a weekly resampling frequency. Only standard calendar is supported for this mode. (:pull:``) Bug fixes ^^^^^^^^^ diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index a69a5f88c..93f26ca0b 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -657,6 +657,8 @@ def preprocess_standardized_index( group = "time.dayofyear" elif compare_offsets(final_freq, "==", "MS"): group = "time.month" + elif compare_offsets(final_freq, "==", "W"): + group = "time.week" else: raise ValueError( f"The input (following resampling if applicable) has a frequency `{final_freq}` " @@ -910,6 +912,9 @@ def reindex_time(da, da_ref, group): elif group == "time.month": da = da.rename(month="time").reindex(time=da_ref.time.dt.month) da["time"] = da_ref.time + elif group == "time.week": + da = da.rename(week="time").reindex(time=da_ref.time.dt.week) + da["time"] = da_ref.time # I don't think rechunking is necessary here, need to check return da if not uses_dask(da) else da.chunk({"time": -1}) From 284323d5de85cf2dd7ef30ed6f2efc2f095f1a0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Fri, 11 Oct 2024 13:20:59 -0400 Subject: [PATCH 2/7] updage changelog --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index f9cc4583a..7185c9fa5 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -23,7 +23,7 @@ New features and enhancements * Indicator parameters can now be assigned a new name, different from the argument name in the compute function. (:pull:`1885`). * ``xclim.indices.run_length.windowed_max_run_sum`` accumulates positive values across runs and yields the the maximum valued run. (:pull:`1926`). * Helper function ``xclim.indices.helpers.make_hourly_temperature`` to estimate hourly temperatures from daily min and max temperatures. (:pull:`1909`). -* `xclim.indices.stats.standardized_index` now supports a weekly resampling frequency. Only standard calendar is supported for this mode. (:pull:``) +* `xclim.indices.stats.standardized_index` now supports a weekly resampling frequency. Only standard calendar is supported for this mode. (:issue:`1892`, :pull:`1952`) Bug fixes ^^^^^^^^^ From 693529f073f6f7c4312baa8c4fa74e9a134cab07 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 17 Oct 2024 10:29:52 -0400 Subject: [PATCH 3/7] weekly tests, remove PWM (not working w/out lmoments3) --- tests/test_indices.py | 37 ++++++++++++++++++++++++++++++++++++- xclim/indices/_agro.py | 4 ++-- xclim/indices/stats.py | 5 +++-- 3 files changed, 41 insertions(+), 5 deletions(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index 0d48e060c..c58b2394b 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -623,6 +623,38 @@ class TestStandardizedIndices: [-0.24417774, -0.11404418, 0.64997039, 1.07670517, 0.6462852], 2e-2, ), + ( + "W", + 1, + "gamma", + "APP", + [0.64820146, 0.04991201, -1.62956493, 1.08898709, -0.01741762], + 2e-2, + ), + ( + "W", + 12, + "gamma", + "APP", + [-1.08683311, -0.47230036, -0.7884111, 0.3341876, 0.06282969], + 2e-2, + ), + ( + "W", + 1, + "gamma", + "ML", + [0.64676962, -0.06904886, -1.60493289, 1.07864037, -0.01415902], + 2e-2, + ), + ( + "W", + 12, + "gamma", + "ML", + [-1.08627775, -0.46491398, -0.77806462, 0.31759127, 0.03794528], + 2e-2, + ), ], ) def test_standardized_precipitation_index( @@ -636,8 +668,11 @@ def test_standardized_precipitation_index( pytest.skip("Skipping SPI/ML/D on older numpy") ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1) if freq == "D": - ds = ds.convert_calendar("366_day", missing=np.nan) # to compare with ``climate_indices`` + ds = ds.convert_calendar("366_day", missing=np.nan) + elif freq == "W": + # only standard calendar supported with freq="W" + ds = ds.convert_calendar("standard", missing=np.nan, align_on="year") pr = ds.pr.sel(time=slice("1998", "2000")) pr_cal = ds.pr.sel(time=slice("1950", "1980")) fitkwargs = {} diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index 02a1a32ad..d337a4f07 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1204,7 +1204,7 @@ def standardized_precipitation_index( :cite:cts:`mckee_relationship_1993` """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( @@ -1298,7 +1298,7 @@ def standardized_precipitation_evapotranspiration_index( """ fitkwargs = fitkwargs or {} - dist_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} if dist in dist_methods: if method not in dist_methods[dist]: raise NotImplementedError( diff --git a/xclim/indices/stats.py b/xclim/indices/stats.py index 93f26ca0b..fef30292e 100644 --- a/xclim/indices/stats.py +++ b/xclim/indices/stats.py @@ -524,6 +524,8 @@ def _fit_start(x, dist: str, **fitkwargs: Any) -> tuple[tuple, dict]: loc0 = (x1 * xn - x2**2) / (x1 + xn - 2 * x2) loc0 = loc0 if loc0 < x1 else (0.9999 * x1 if x1 > 0 else 1.0001 * x1) x_pos = x - loc0 + # TODO: change this? + # not necessary for log-logistic, according to SPEI package x_pos = x_pos[x_pos > 0] # method of moments: # LHS is computed analytically with the two-parameters log-logistic distribution @@ -758,8 +760,7 @@ def standardized_index_fit_params( "Pass a value for `floc` in `fitkwargs`." ) - # "WPM" method doesn't seem to work for gamma or pearson3 - dist_and_methods = {"gamma": ["ML", "APP", "PWM"], "fisk": ["ML", "APP"]} + dist_and_methods = {"gamma": ["ML", "APP"], "fisk": ["ML", "APP"]} dist = get_dist(dist) if dist.name not in dist_and_methods: raise NotImplementedError(f"The distribution `{dist.name}` is not supported.") From 3af6756731949f39ea1ec1c51bbcb9a0378afd2f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Thu, 17 Oct 2024 10:31:19 -0400 Subject: [PATCH 4/7] Update CHANGELOG.rst Co-authored-by: Trevor James Smith <10819524+Zeitsperre@users.noreply.github.com> --- CHANGELOG.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index fc8be7c53..72c04f10f 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -28,7 +28,7 @@ New features and enhancements * ``xclim.indices.run_length.runs_with_holes`` allows to input a condition that must be met for a run to start and a second condition that must be met for the run to stop. (:pull:`1778`). * New generic compute function ``xclim.indices.generic.thresholded_events`` that finds events based on a threshold condition and returns basic stats for each. See also ``xclim.indices.run_length.find_events``. (:pull:`1778`). * ``xclim.core.units.rate2amount`` and ``xclim.core.units.amount2rate`` can now also accept quantities (pint objects or strings), in which case the ``dim`` argument must be the ``time`` coordinate through which we can find the sampling rate. (:pull:`1778`). -* `xclim.indices.stats.standardized_index` now supports a weekly resampling frequency. Only standard calendar is supported for this mode. (:issue:`1892`, :pull:`1952`) +* ``xclim.indices.stats.standardized_index`` now supports a weekly resampling frequency. Only "standard" calendar is supported for this mode. (:issue:`1892`, :pull:`1952`) Bug fixes ^^^^^^^^^ From 8c74372e9944f2f40c1dfa6a681bf2981e3b65f0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= Date: Thu, 17 Oct 2024 10:39:09 -0400 Subject: [PATCH 5/7] better docstrings --- xclim/indices/_agro.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index d337a4f07..f19f1a413 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1162,6 +1162,8 @@ def standardized_precipitation_index( * N-month SPI / N-day SPI is determined by choosing the `window = N` and the appropriate frequency `freq`. * Supported statistical distributions are: ["gamma", "fisk"], where "fisk" is scipy's implementation of a log-logistic distribution + * Supported frequencies are daily ("D"), weekly ("W"), and monthly ("MS"). + Weekly frequency will only work if the input array has a "standard" calendar. * If `params` is given as input, it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. * "APP" method only supports two-parameter distributions. Parameter `loc` needs to be fixed to use method `APP`. * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in @@ -1267,10 +1269,8 @@ def standardized_precipitation_evapotranspiration_index( dist : {'gamma', 'fisk'} Name of the univariate distribution. (see :py:mod:`scipy.stats`). method : {'APP', 'ML'} - Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate), or - `PWM` (probability weighted moments). - The approximate method uses a deterministic function that doesn't involve any optimization. Available methods - vary with the distribution: 'gamma':{'APP', 'ML', 'PWM'}, 'fisk':{'APP', 'ML'} + Name of the fitting method, such as `ML` (maximum likelihood), `APP` (approximate). The approximate method + uses a deterministic function that doesn't involve any optimization. fitkwargs : dict, optional Kwargs passed to ``xclim.indices.stats.fit`` used to impose values of certains parameters (`floc`, `fscale`). cal_start : DateStr, optional From d9f006531cf98d365eb2cbc5f0ca514305dfa6e3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=C3=89ric=20Dupuis?= <71575674+coxipi@users.noreply.github.com> Date: Thu, 17 Oct 2024 11:30:09 -0400 Subject: [PATCH 6/7] better calendar options and doc Co-authored-by: Pascal Bourgault --- CHANGELOG.rst | 2 +- tests/test_indices.py | 2 +- xclim/indices/_agro.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 49af2b195..2a7db84cb 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -29,7 +29,7 @@ New features and enhancements * ``xclim.indices.run_length.runs_with_holes`` allows to input a condition that must be met for a run to start and a second condition that must be met for the run to stop. (:pull:`1778`). * New generic compute function ``xclim.indices.generic.thresholded_events`` that finds events based on a threshold condition and returns basic stats for each. See also ``xclim.indices.run_length.find_events``. (:pull:`1778`). * ``xclim.core.units.rate2amount`` and ``xclim.core.units.amount2rate`` can now also accept quantities (pint objects or strings), in which case the ``dim`` argument must be the ``time`` coordinate through which we can find the sampling rate. (:pull:`1778`). -* ``xclim.indices.stats.standardized_index`` now supports a weekly resampling frequency. Only "standard" calendar is supported for this mode. (:issue:`1892`, :pull:`1952`) +* ``xclim.indices.stats.standardized_index`` now supports a weekly resampling frequency. Only "standard" calendars using numpy's `datetime64` dtype are supported for this mode. (:issue:`1892`, :pull:`1952`) Bug fixes ^^^^^^^^^ diff --git a/tests/test_indices.py b/tests/test_indices.py index 0a75f5aaf..b774e9083 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -672,7 +672,7 @@ def test_standardized_precipitation_index( ds = ds.convert_calendar("366_day", missing=np.nan) elif freq == "W": # only standard calendar supported with freq="W" - ds = ds.convert_calendar("standard", missing=np.nan, align_on="year") + ds = ds.convert_calendar("standard", missing=np.nan, align_on="year", use_cftime=False) pr = ds.pr.sel(time=slice("1998", "2000")) pr_cal = ds.pr.sel(time=slice("1950", "1980")) fitkwargs = {} diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py index f19f1a413..3a9cf0139 100644 --- a/xclim/indices/_agro.py +++ b/xclim/indices/_agro.py @@ -1163,7 +1163,7 @@ def standardized_precipitation_index( * Supported statistical distributions are: ["gamma", "fisk"], where "fisk" is scipy's implementation of a log-logistic distribution * Supported frequencies are daily ("D"), weekly ("W"), and monthly ("MS"). - Weekly frequency will only work if the input array has a "standard" calendar. + Weekly frequency will only work if the input array has a "standard" (non-cftime) calendar. * If `params` is given as input, it overrides the `cal_start`, `cal_end`, `freq` and `window`, `dist` and `method` options. * "APP" method only supports two-parameter distributions. Parameter `loc` needs to be fixed to use method `APP`. * The standardized index is bounded by ±8.21. 8.21 is the largest standardized index as constrained by the float64 precision in From fd80bdf366fab32954a4b21bf4ffebb806ca120c Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 17 Oct 2024 15:32:35 +0000 Subject: [PATCH 7/7] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- tests/test_indices.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/tests/test_indices.py b/tests/test_indices.py index b774e9083..2b707f8c1 100644 --- a/tests/test_indices.py +++ b/tests/test_indices.py @@ -672,7 +672,9 @@ def test_standardized_precipitation_index( ds = ds.convert_calendar("366_day", missing=np.nan) elif freq == "W": # only standard calendar supported with freq="W" - ds = ds.convert_calendar("standard", missing=np.nan, align_on="year", use_cftime=False) + ds = ds.convert_calendar( + "standard", missing=np.nan, align_on="year", use_cftime=False + ) pr = ds.pr.sel(time=slice("1998", "2000")) pr_cal = ds.pr.sel(time=slice("1950", "1980")) fitkwargs = {}