fix freq='D' behaviour of standardized indices

CHANGES.rst

@@ -29,6 +29,7 @@ Bug fixes
 ^^^^^^^^^
 * Fixed passing ``missing=0`` to ``xclim.core.calendar.convert_calendar``. (:issue:`1562`, :pull:`1563`).
 * Fix wrong `window` attributes in ``xclim.indices.standardized_precipitation_index``, ``xclim.indices.standardized_precipitation_evapotranspiration_index``. (:issue:`1552`  :pull:`1554`).
+* Fix the daily case `freq='D'` of ``xclim.stats.preprocess_standardized_index`` (:issue:`1602`  :pull:`1607`).
 * Several spelling mistakes have been corrected within the documentation and codebase. (:pull:`1576`).
 
 Internal changes

tests/test_indices.py

@@ -21,7 +21,7 @@
 import xarray as xr
 
 from xclim import indices as xci
-from xclim.core.calendar import date_range, percentile_doy
+from xclim.core.calendar import convert_calendar, date_range, percentile_doy
 from xclim.core.options import set_options
 from xclim.core.units import ValidationError, convert_units_to, units
 
@@ -456,13 +456,13 @@ def test_effective_growing_degree_days(
 
         np.testing.assert_array_equal(out, np.array([np.NaN, expected]))
 
-    # gamma reference results: Obtained with `monocongo/climate_indices` library
-    # fisk reference results: Obtained with R package `SPEI`
+    # gamma/APP reference results: Obtained with `monocongo/climate_indices` library
+    # MS/fisk/ML reference results: Obtained with R package `SPEI`
     # Using the method `APP` in XClim matches the method from monocongo, hence the very low
     # tolerance possible.
     # Repeated tests with lower tolerance means we want a more precise comparison, so we compare
     # the current version of XClim with the version where the test was implemented
-    # TODO : Add tests for SPI_daily.
+    @pytest.mark.slow
     @pytest.mark.parametrize(
         "freq, window, dist, method,  values, diff_tol",
         [
@@ -518,12 +518,96 @@ def test_effective_growing_degree_days(
                 [0.683273, 1.51189, 1.61597, 1.03875, 0.72531],
                 2e-2,
             ),
+            (
+                "D",
+                1,
+                "gamma",
+                "APP",
+                [-0.18618353, 1.44582971, 0.95985043, 0.15779587, -0.37801587],
+                2e-2,
+            ),
+            (
+                "D",
+                12,
+                "gamma",
+                "APP",
+                [-0.24417774, -0.11404418, 0.64997039, 1.07670517, 0.6462852],
+                2e-2,
+            ),
+            (
+                "D",
+                1,
+                "gamma",
+                "ML",
+                [-0.03577971, 1.30589409, 0.8863447, 0.23906544, -0.05185997],
+                2e-2,
+            ),
+            (
+                "D",
+                12,
+                "gamma",
+                "ML",
+                [-0.15846245, -0.04924534, 0.66299367, 1.09938471, 0.66095752],
+                2e-2,
+            ),
+            (
+                "D",
+                1,
+                "fisk",
+                "APP",
+                [-1.26216389, 1.03096183, 0.62985354, -0.50335153, -1.32788296],
+                2e-2,
+            ),
+            (
+                "D",
+                12,
+                "fisk",
+                "APP",
+                [-0.57109258, -0.40657737, 0.55163493, 0.97381067, 0.55580649],
+                2e-2,
+            ),
+            (
+                "D",
+                1,
+                "fisk",
+                "ML",
+                [-0.05562691, 1.30809152, 0.6954986, 0.33018744, -0.50258979],
+                2e-2,
+            ),
+            (
+                "D",
+                12,
+                "fisk",
+                "ML",
+                [-0.14151269, -0.01914608, 0.7080277, 1.01510279, 0.6954002],
+                2e-2,
+            ),
+            (
+                None,
+                1,
+                "gamma",
+                "APP",
+                [-0.18618353, 1.44582971, 0.95985043, 0.15779587, -0.37801587],
+                2e-2,
+            ),
+            (
+                None,
+                12,
+                "gamma",
+                "APP",
+                [-0.24417774, -0.11404418, 0.64997039, 1.07670517, 0.6462852],
+                2e-2,
+            ),
         ],
     )
     def test_standardized_precipitation_index(
         self, open_dataset, freq, window, dist, method, values, diff_tol
     ):
         ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1)
+        if freq == "D":
+            ds = convert_calendar(
+                ds, "366_day", missing=np.NaN
+            )  # to compare with ``climate_indices``
         pr = ds.pr.sel(time=slice("1998", "2000"))
         pr_cal = ds.pr.sel(time=slice("1950", "1980"))
         params = xci.stats.standardized_index_fit_params(
@@ -545,6 +629,7 @@ def test_standardized_precipitation_index(
         np.testing.assert_allclose(spi.values, values, rtol=0, atol=diff_tol)
 
     # See SPI version
+    @pytest.mark.slow
     @pytest.mark.parametrize(
         "freq, window, dist, method,  values, diff_tol",
         [

xclim/indices/_agro.py

@@ -1238,7 +1238,7 @@ def standardized_precipitation_index(
 
     spi = standardized_index(pr, params)
     spi.attrs = params.attrs
-    spi.attrs["freq"] = freq or xarray.infer_freq(spi.time)
+    spi.attrs["freq"] = (freq or xarray.infer_freq(spi.time)) or "undefined"
     spi.attrs["window"] = window
     spi.attrs["units"] = ""
     return spi

xclim/indices/stats.py

@@ -8,7 +8,7 @@
 import numpy as np
 import xarray as xr
 
-from xclim.core.calendar import resample_doy, select_time
+from xclim.core.calendar import compare_offsets, resample_doy, select_time
 from xclim.core.formatting import prefix_attrs, unprefix_attrs, update_history
 from xclim.core.units import convert_units_to
 from xclim.core.utils import Quantified, uses_dask
@@ -625,13 +625,21 @@ def preprocess_standardized_index(
     # We could allow a more general frequency in this function and move
     # the constraint {"D", "MS"} in specific indices such as SPI / SPEI.
     final_freq = freq or xr.infer_freq(da.time)
-    try:
-        group = {"D": "time.dayofyear", "MS": "time.month"}[final_freq]
-    except KeyError():
-        raise ValueError(
-            f"The input (following resampling if applicable) has a frequency `{final_freq}`"
-            "which is not supported for standardized indices."
+    if final_freq:
+        if final_freq == "D":
+            group = "time.dayofyear"
+        elif compare_offsets(final_freq, "==", "MS"):
+            group = "time.month"
+        else:
+            raise ValueError(
+                f"The input (following resampling if applicable) has a frequency `{final_freq}` "
+                "which is not supported for standardized indices."
+            )
+    else:
+        warnings.warn(
+            "No resampling frequency was specified and a frequency for the dataset could not be identified with ``xr.infer_freq``"
         )
+        group = "time.dayofyear"
 
     if freq is not None:
         da = da.resample(time=freq).mean(keep_attrs=True)
@@ -732,10 +740,7 @@ def standardized_index_fit_params(
         "units": "",
         "offset": offset or "",
     }
-    if indexer != {}:
-        method, args = indexer.popitem()
-    else:
-        method, args = "", []
+    method, args = ("", []) if indexer == {} else indexer.popitem()
     params.attrs["time_indexer"] = (method, *args)
 
     return params
@@ -762,8 +767,6 @@ def standardized_index(da: xr.DataArray, params: xr.DataArray):
 
     def reindex_time(da, da_ref):
         if group == "time.dayofyear":
-            da = da.rename(day="time").reindex(time=da_ref.time.dt.dayofyear)
-            da["time"] = da_ref.time
             da = resample_doy(da, da_ref)
         elif group == "time.month":
             da = da.rename(month="time").reindex(time=da_ref.time.dt.month)