Improve SPI performance

CHANGES.rst

@@ -6,10 +6,22 @@ v0.46.0
 --------------------
 Contributors to this version: Éric Dupuis (:user:`coxipi`).
 
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* Standardized indices (``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index``)  (:issue:`1270`, :issue:`1416`, :issue:`1474`, :pull:`1311`) were changed:
+    * More optimized and noticeably faster calculation.
+    * Can be computed in two steps: First compute fit parameters with ``xclim.indices.standardized_index_fit_params``, then use the output in the standardized indices functions.
+    * The standardized index values are now clipped to ± 8.21. This reflects the ``float64`` precision of the computation when cumulative distributed function values are inverted to a normal distribution.
+    * An offset parameter is now accessible ``xclim.indices.standardized_precipitation_evapotranspiration_index``.
+
 Internal changes
 ^^^^^^^^^^^^^^^^
 * Change "degK" to "K" used to design Kelvin units. (:pull:`1475`).
 
+Breaking changes
+^^^^^^^^^^^^^^^^
+* ``xclim.indices.standardized_precipitation_index`` and ``xclim.indices.standardized_precipitation_evapotranspiration_index`` will no longer accept two datasets (data and calibration data). Instead, a single dataset along start and end of calibration dates (`cal_start`, `cal_end`) is expected. (:issue:`1270` ,:pull:`1311`)
+
 v0.45.0 (2023-09-05)
 --------------------
 Contributors to this version: David Huard (:user:`huard`), Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Juliette Lavoie (:user:`juliettelavoie`), Gabriel Rondeau-Genesse (:user:`RondeauG`), Marco Braun (:user:`vindelico`), Éric Dupuis (:user:`coxipi`).

tests/test_indices.py

@@ -526,10 +526,10 @@ def test_standardized_precipitation_index(
         ds = open_dataset("sdba/CanESM2_1950-2100.nc").isel(location=1)
         pr = ds.pr.sel(time=slice("1998", "2000"))
         pr_cal = ds.pr.sel(time=slice("1950", "1980"))
-        spi = xci.standardized_precipitation_index(
-            pr, pr_cal, freq, window, dist, method
+        params = xci.standardized_index_fit_params(
+            pr_cal, freq=freq, window=window, dist=dist, method=method
         )
-
+        spi = xci.standardized_precipitation_index(pr, params=params)
         # Only a few moments before year 2000 are tested
         spi = spi.isel(time=slice(-11, -1, 2))
 
@@ -610,11 +610,17 @@ def test_standardized_precipitation_evapotranspiration_index(
             tas = tasmax - 2.5
             tasmin = tasmax - 5
             wb = xci.water_budget(pr, None, tasmin, tasmax, tas)
-            wb_cal = wb.sel(time=slice("1950", "1980"))
-            wb = wb.sel(time=slice("1998", "2000"))
 
+        params = xci.standardized_index_fit_params(
+            wb.sel(time=slice("1950", "1980")),
+            freq=freq,
+            window=window,
+            dist=dist,
+            method=method,
+            offset="1 mm/d",
+        )
         spei = xci.standardized_precipitation_evapotranspiration_index(
-            wb, wb_cal, freq, window, dist, method
+            wb.sel(time=slice("1998", "2000")), params=params
         )
 
         # Only a few moments before year 2000 are tested
@@ -625,6 +631,51 @@ def test_standardized_precipitation_evapotranspiration_index(
 
         np.testing.assert_allclose(spei.values, values, rtol=0, atol=diff_tol)
 
+    def test_standardized_index_modularity(self, open_dataset):
+        freq, window, dist, method = "MS", 6, "gamma", "APP"
+        ds = (
+            open_dataset("sdba/CanESM2_1950-2100.nc")
+            .isel(location=1)
+            .sel(time=slice("1950", "2000"))
+        )
+        pr = ds.pr
+        # generate water budget
+        with xr.set_options(keep_attrs=True):
+            tasmax = ds.tasmax
+            tas = tasmax - 2.5
+            tasmin = tasmax - 5
+            wb = xci.water_budget(pr, None, tasmin, tasmax, tas)
+
+        params = xci.standardized_index_fit_params(
+            wb.sel(time=slice("1950", "1980")),
+            freq=freq,
+            window=window,
+            dist=dist,
+            method=method,
+            offset="1 mm/d",
+        )
+        spei1 = xci.standardized_precipitation_evapotranspiration_index(
+            wb.sel(time=slice("1998", "2000")), params=params
+        )
+
+        spei2 = xci.standardized_precipitation_evapotranspiration_index(
+            wb,
+            freq=freq,
+            window=window,
+            dist=dist,
+            method=method,
+            offset="1 mm/d",
+            cal_start="1950",
+            cal_end="1980",
+        ).sel(time=slice("1998", "2000"))
+
+        # In the previous computation, the first {window-1} values are NaN because the rolling is performed on the period [1998,2000].
+        # Here, the computation is performed on the period [1950,2000], *then* subsetted to [1998,2000], so it doesn't have NaNs
+        # for the first values
+        spei2[{"time": slice(0, window - 1)}] = np.nan
+
+        np.testing.assert_allclose(spei1.values, spei2.values, rtol=0, atol=1e-4)
+
 
 class TestDailyFreezeThawCycles:
     @pytest.mark.parametrize(