Support freq="W" in standardized indices (#1952)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [ ] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #xyz
- [ ] Tests for the changes have been added (for bug fixes / features)
- [ ] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [ ] CHANGELOG.rst has been updated (with summary of main changes)
- [ ] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

### What kind of change does this PR introduce?

* Add support for weekly standardized indices

### Does this PR introduce a breaking change?

No

### Other information:

I realize that xclim and climate_indices treat zero-inflated
distributions differently. I will investigate which library does it
correctly and include changes if needed

EDIT: I believe `xclim` has the correct implementation. The probability
of zeroes should be determined in the fitting procedure, not by using
the full dataset. The idea is when you compute a CDF, a zero value
should be mapped to the `prob_of_zero` in the calibration period. That
is the same logic as using the fitting params of the distribution in the
calibration period to compute the CDF in the full dataset.

I confirmed this is also how it's done in the case of the R package
(SPEI), although I find some steps a bit weird, but I'm not too familiar
with R.

CHANGELOG.rst

@@ -29,6 +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" calendars using numpy's `datetime64` dtype are supported for this mode. (:issue:`1892`, :pull:`1952`)
 
 Bug fixes
 ^^^^^^^^^

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,13 @@ 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", use_cftime=False
+ )
  pr = ds.pr.sel(time=slice("1998", "2000"))
  pr_cal = ds.pr.sel(time=slice("1950", "1980"))
  fitkwargs = {}

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" (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
@@ -1204,7 +1206,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(
@@ -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
@@ -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(

xclim/indices/stats.py

@@ -540,6 +540,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
@@ -675,6 +677,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}` "
@@ -775,8 +779,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.")
@@ -929,6 +932,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})