From b1fbbf796f16f61c0d6c21d611ddf2c7409917cd Mon Sep 17 00:00:00 2001
From: yrouranos <rousseau.yannick@ouranos.ca>
Date: Mon, 25 Oct 2021 12:29:28 -0400
Subject: [PATCH 1/9] Improved algorithm for index 'dry_spell_total_length'
 with new options (operator 'max', possibility to specify start/end dates) and
 calculation for days near dataset boundaries.

---
 xclim/indices/_agro.py              | 59 ++++++++++++++++++++++++++---
 xclim/testing/tests/test_indices.py | 10 ++++-
 xclim/testing/tests/test_precip.py  | 23 +++++++++--
 3 files changed, 81 insertions(+), 11 deletions(-)

diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py
index 8df821607..1898d30d9 100644
--- a/xclim/indices/_agro.py
+++ b/xclim/indices/_agro.py
@@ -1,5 +1,6 @@
 # noqa: D100
 
+import datetime
 from typing import Optional
 
 import numpy as np
@@ -628,11 +629,17 @@ def dry_spell_frequency(
 
 @declare_units(pr="[precipitation]", thresh="[length]")
 def dry_spell_total_length(
-    pr: xarray.DataArray, thresh: str = "1.0 mm", window: int = 3, freq: str = "YS"
+    pr: xarray.DataArray,
+    thresh: str = "1.0 mm",
+    window: int = 3,
+    op: str = "sum",
+    freq: str = "YS",
+    start_date: str = "",
+    end_date: str = "",
 ) -> xarray.DataArray:
     """
-    Return the total number of days in dry periods of n days and more, during which the accumulated precipitation
-    on a window of n days is under the threshold.
+    Return the total number of days in dry periods of n days and more, during which the maximum or accumulated
+    precipitation on a window of n days is under the threshold.
 
     Parameters
     ----------
@@ -641,9 +648,15 @@ def dry_spell_total_length(
     thresh : str
       Accumulated precipitation value under which a period is considered dry.
     window : int
-      Number of days where the accumulated precipitation is under threshold.
+      Number of days where the maximum or accumulated precipitation is under threshold.
+    op : {"max", "sum"}
+        Reduce operation.
     freq : str
       Resampling frequency.
+    start_date : str
+        First day of year to consider ("mm-dd").
+    end_date : str
+        Last day of year to consider ("mm-dd").
 
     Returns
     -------
@@ -653,8 +666,42 @@ def dry_spell_total_length(
     pram = rate2amount(pr, out_units="mm")
     thresh = convert_units_to(thresh, pram)
 
-    mask = pram.rolling(time=window, center=True).sum() < thresh
-    out = (mask.rolling(time=window, center=True).sum() >= 1).resample(time=freq).sum()
+    pram = xarray.where(pram < (thresh if op == "max" else 0), 0, pram)
+    pram.attrs["units"] = "mm"
+
+    dry = None
+    for i in range(2):
+        pram_i = pram if dry is None else pram.sortby("time", ascending=False)
+        if op == "max":
+            mask_i = xarray.DataArray(pram_i.rolling(time=window).max() < thresh)
+        else:
+            mask_i = xarray.DataArray(pram_i.rolling(time=window).sum() < thresh)
+        dry_i = xarray.DataArray(mask_i.rolling(time=window).sum() >= 1).shift(
+            time=-(window - 1)
+        )
+        if dry is None:
+            dry = dry_i
+        else:
+            dry_i = dry_i.sortby("time", ascending=True)
+            dt = (dry.time - dry.time[0]).dt.days
+            dry = xarray.where(dt > len(dry.time) - window - 1, dry_i, dry)
+
+    doy_start = 1
+    if start_date != "":
+        doy_start = datetime.datetime.strptime(start_date, "%m-%d").timetuple().tm_yday
+    doy_end = 365
+    if end_date != "":
+        doy_end = datetime.datetime.strptime(end_date, "%m-%d").timetuple().tm_yday
+    if doy_end >= doy_start:
+        doy = (pram.time.dt.dayofyear >= doy_start) & (
+            pram.time.dt.dayofyear <= doy_end
+        )
+    else:
+        doy = (pram.time.dt.dayofyear <= doy_end) | (
+            pram.time.dt.dayofyear >= doy_start
+        )
+
+    out = (dry & doy).astype(float).resample(time=freq).sum("time")
 
     return to_agg_units(out, pram, "count")
 
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index 17fb43b93..cd2bd9c77 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -2290,7 +2290,13 @@ def test_dry_spell(pr_series):
     pr.attrs["units"] = "mm/day"
 
     events = xci.dry_spell_frequency(pr, thresh="3 mm", window=7, freq="YS")
-    total_d = xci.dry_spell_total_length(pr, thresh="3 mm", window=7, freq="YS")
+    total_d_sum = xci.dry_spell_total_length(
+        pr, thresh="3 mm", window=7, op="sum", freq="YS"
+    )
+    total_d_max = xci.dry_spell_total_length(
+        pr, thresh="3 mm", window=7, op="max", freq="YS"
+    )
 
     np.testing.assert_allclose(events[0], [2], rtol=1e-1)
-    np.testing.assert_allclose(total_d[0], [12], rtol=1e-1)
+    np.testing.assert_allclose(total_d_sum[0], [12], rtol=1e-1)
+    np.testing.assert_allclose(total_d_max[0], [20], rtol=1e-1)
diff --git a/xclim/testing/tests/test_precip.py b/xclim/testing/tests/test_precip.py
index 71d15bd8b..9fa023c82 100644
--- a/xclim/testing/tests/test_precip.py
+++ b/xclim/testing/tests/test_precip.py
@@ -446,10 +446,23 @@ def test_dry_spell():
     pr = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc").pr
 
     events = atmos.dry_spell_frequency(pr, thresh="3 mm", window=7, freq="YS")
-    total_d = atmos.dry_spell_total_length(pr, thresh="3 mm", window=7, freq="YS")
+    total_d_sum = atmos.dry_spell_total_length(
+        pr, thresh="3 mm", window=7, op="sum", freq="YS"
+    )
+    total_d_max = atmos.dry_spell_total_length(
+        pr, thresh="3 mm", window=7, op="max", freq="YS"
+    )
+
+    total_d_sum = list(
+        total_d_sum.sel(location=slice("Halifax")).squeeze()[0:2].astype(int).values
+    )
+    total_d_max = list(
+        total_d_max.sel(location=slice("Halifax")).squeeze()[0:2].astype(int).values
+    )
 
     np.testing.assert_allclose(events[0:2, 0], [5, 8], rtol=1e-1)
-    np.testing.assert_allclose(total_d[0:2, 0], [50, 67], rtol=1e-1)
+    np.testing.assert_allclose(total_d_sum, [50, 60], rtol=1e-1)
+    np.testing.assert_allclose(total_d_max, [76, 97], rtol=1e-1)
 
     assert (
         "The annual number of dry periods of 7 days and more, during which the accumulated "
@@ -457,5 +470,9 @@ def test_dry_spell():
     ) in events.description
     assert (
         "The annual number of days in dry periods of 7 days and more"
-        in total_d.description
+        in total_d_sum.description
+    )
+    assert (
+        "The annual number of days in dry periods of 7 days and more"
+        in total_d_max.description
     )

From 74d3a3b2fed4356e288e4767fda26c9f4e14f83d Mon Sep 17 00:00:00 2001
From: yrouranos <rousseau.yannick@ouranos.ca>
Date: Mon, 25 Oct 2021 13:47:10 -0400
Subject: [PATCH 2/9] Improved algorithm for index 'dry_spell_total_length'
 with new options (operator 'max', possibility to specify start/end dates) and
 calculation for days near dataset boundaries.

---
 AUTHORS.rst | 1 +
 1 file changed, 1 insertion(+)

diff --git a/AUTHORS.rst b/AUTHORS.rst
index 755e26856..ad902b9f2 100644
--- a/AUTHORS.rst
+++ b/AUTHORS.rst
@@ -27,3 +27,4 @@ Contributors
 * Clair Barnes <clair.barnes.16@ucl.ac.uk> `@clairbarnes <https://github.com/clairbarnes>`_
 * Raquel Alegre <raquel.alegre@gmail.com> `@raquel-ucl <https://github.com/raquel-ucl>`_
 * Abel Aoun <aoun.abel@gmail.com> `@bzah <https://github.com/bzah>`_
+* Yannick Rousseau `@yrouranos <https://github.com/yrouranos>`_

From 21733d0616aa3b950bf22966cbf0752a6a1bf7ab Mon Sep 17 00:00:00 2001
From: yrouranos <rousseau.yannick@ouranos.ca>
Date: Thu, 4 Nov 2021 20:22:23 -0400
Subject: [PATCH 3/9] Improved algorithm for index 'dry_spell_total_length'
 with new options (operator 'max', possibility to specify start/end dates) and
 calculation for days near dataset boundaries.

---
 xclim/indices/_agro.py              | 30 +++++----
 xclim/testing/tests/test_indices.py | 95 ++++++++++++++++++-----------
 xclim/testing/tests/test_precip.py  |  8 +--
 3 files changed, 79 insertions(+), 54 deletions(-)

diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py
index 1898d30d9..2e0f1aafb 100644
--- a/xclim/indices/_agro.py
+++ b/xclim/indices/_agro.py
@@ -673,12 +673,10 @@ def dry_spell_total_length(
     for i in range(2):
         pram_i = pram if dry is None else pram.sortby("time", ascending=False)
         if op == "max":
-            mask_i = xarray.DataArray(pram_i.rolling(time=window).max() < thresh)
+            mask_i = pram_i.rolling(time=window).max() < thresh
         else:
-            mask_i = xarray.DataArray(pram_i.rolling(time=window).sum() < thresh)
-        dry_i = xarray.DataArray(mask_i.rolling(time=window).sum() >= 1).shift(
-            time=-(window - 1)
-        )
+            mask_i = pram_i.rolling(time=window).sum() < thresh
+        dry_i = (mask_i.rolling(time=window).sum() >= 1).shift(time=-(window - 1))
         if dry is None:
             dry = dry_i
         else:
@@ -692,14 +690,20 @@ def dry_spell_total_length(
     doy_end = 365
     if end_date != "":
         doy_end = datetime.datetime.strptime(end_date, "%m-%d").timetuple().tm_yday
-    if doy_end >= doy_start:
-        doy = (pram.time.dt.dayofyear >= doy_start) & (
-            pram.time.dt.dayofyear <= doy_end
-        )
-    else:
-        doy = (pram.time.dt.dayofyear <= doy_end) | (
-            pram.time.dt.dayofyear >= doy_start
-        )
+    bisextile_year_fix = xarray.where(pram.time.dt.year % 4 == 0, 1, 0)
+    doy_start = (
+        ([doy_start] * len(pram.time))
+        if doy_start < 60
+        else (doy_start + bisextile_year_fix)
+    )
+    doy_end = (
+        ([doy_end] * len(pram.time)) if doy_end < 60 else (doy_end + bisextile_year_fix)
+    )
+    doy = xarray.where(
+        doy_end >= doy_start,
+        (pram.time.dt.dayofyear >= doy_start) & (pram.time.dt.dayofyear <= doy_end),
+        (pram.time.dt.dayofyear <= doy_end) | (pram.time.dt.dayofyear >= doy_start),
+    )
 
     out = (dry & doy).astype(float).resample(time=freq).sum("time")
 
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index cd2bd9c77..ae2c3f58b 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -2261,42 +2261,67 @@ def test_water_budget(pr_series, tasmin_series, tasmax_series):
 
 
 def test_dry_spell(pr_series):
-    pr = pr_series(
-        np.array(
-            [
-                1.01,
-                1.01,
-                1.01,
-                1.01,
-                1.01,
-                1.01,
-                0.01,
-                0.01,
-                0.01,
-                0.51,
-                0.51,
-                0.75,
-                0.75,
-                0.51,
-                0.01,
-                0.01,
-                0.01,
-                1.01,
-                1.01,
-                1.01,
-            ]
-        )
-    )
-    pr.attrs["units"] = "mm/day"
 
-    events = xci.dry_spell_frequency(pr, thresh="3 mm", window=7, freq="YS")
-    total_d_sum = xci.dry_spell_total_length(
-        pr, thresh="3 mm", window=7, op="sum", freq="YS"
+    pr_arr_1 = (
+        [1.01] * 6
+        + [0.01] * 3
+        + [0.51] * 2
+        + [0.75] * 2
+        + [0.51]
+        + [0.01] * 3
+        + [1.01] * 3
     )
-    total_d_max = xci.dry_spell_total_length(
-        pr, thresh="3 mm", window=7, op="max", freq="YS"
+    pr_arr_2 = (
+        [0.01] * 6
+        + [1.01] * 3
+        + [0.51] * 2
+        + [0.75] * 2
+        + [0.51]
+        + [0.01] * 3
+        + [0.01] * 3
     )
+    pr_arr_3 = [3.01] * 358 + [0.99] * 14 + [3.01] * 358
+    pr_1 = pr_series(np.array(pr_arr_1))
+    pr_2 = pr_series(np.array(pr_arr_2))
+    pr_3 = pr_series(np.array(pr_arr_3), start="1981-01-01")
+    pr_1.attrs["units"] = "mm/day"
+    pr_2.attrs["units"] = "mm/day"
+    pr_3.attrs["units"] = "mm/day"
+
+    for i in range(1, 4):
+
+        if i == 1:
+            thresh, window = 3, 7
+            out_events = [2]
+            out_total_d_sum = [12]
+            out_total_d_max = [20]
+            pr = pr_1
+        elif i == 2:
+            thresh, window = 3, 7
+            out_events = [2]
+            out_total_d_sum = [18]
+            out_total_d_max = [20]
+            pr = pr_2
+        else:
+            thresh, window = 1, 14
+            out_events = [0]
+            out_total_d_sum = out_total_d_max = [7, 7]
+            pr = pr_3
+
+        events = xci.dry_spell_frequency(
+            pr, thresh=str(thresh) + " mm", window=window, freq="YS"
+        )
+        total_d_sum = xci.dry_spell_total_length(
+            pr,
+            thresh=str(thresh * (window if i == 3 else 1)) + " mm",
+            window=window,
+            op="sum",
+            freq="YS",
+        )
+        total_d_max = xci.dry_spell_total_length(
+            pr, thresh=str(thresh) + " mm", window=window, op="max", freq="YS"
+        )
 
-    np.testing.assert_allclose(events[0], [2], rtol=1e-1)
-    np.testing.assert_allclose(total_d_sum[0], [12], rtol=1e-1)
-    np.testing.assert_allclose(total_d_max[0], [20], rtol=1e-1)
+        np.testing.assert_allclose(events[0], out_events, rtol=1e-1)
+        np.testing.assert_allclose(total_d_sum[0], out_total_d_sum, rtol=1e-1)
+        np.testing.assert_allclose(total_d_max[0], out_total_d_max, rtol=1e-1)
diff --git a/xclim/testing/tests/test_precip.py b/xclim/testing/tests/test_precip.py
index 9fa023c82..f1302fb3e 100644
--- a/xclim/testing/tests/test_precip.py
+++ b/xclim/testing/tests/test_precip.py
@@ -453,12 +453,8 @@ def test_dry_spell():
         pr, thresh="3 mm", window=7, op="max", freq="YS"
     )
 
-    total_d_sum = list(
-        total_d_sum.sel(location=slice("Halifax")).squeeze()[0:2].astype(int).values
-    )
-    total_d_max = list(
-        total_d_max.sel(location=slice("Halifax")).squeeze()[0:2].astype(int).values
-    )
+    total_d_sum = total_d_sum.sel(location=slice("Halifax"), drop=True).squeeze()[0:2]
+    total_d_max = total_d_max.sel(location=slice("Halifax"), drop=True).squeeze()[0:2]
 
     np.testing.assert_allclose(events[0:2, 0], [5, 8], rtol=1e-1)
     np.testing.assert_allclose(total_d_sum, [50, 60], rtol=1e-1)

From a540a79e6b30f643b6217bd09815d6a9d8d00c03 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Wed, 1 Dec 2021 11:30:52 -0500
Subject: [PATCH 4/9] Adapt to new indexer - modifications to simplify code

---
 xclim/indices/_agro.py              | 59 ++++++-----------------------
 xclim/testing/tests/test_indices.py | 25 +++++++-----
 xclim/testing/tests/test_precip.py  |  4 +-
 3 files changed, 29 insertions(+), 59 deletions(-)

diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py
index 598e2dca9..5f80e35fd 100644
--- a/xclim/indices/_agro.py
+++ b/xclim/indices/_agro.py
@@ -11,7 +11,7 @@
 from xclim.core.units import convert_units_to, declare_units, rate2amount, to_agg_units
 from xclim.core.utils import DayOfYearStr
 from xclim.indices._threshold import first_day_above, first_day_below, freshet_start
-from xclim.indices.generic import aggregate_between_dates, day_lengths
+from xclim.indices.generic import aggregate_between_dates, day_lengths, select_time
 
 # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start
 # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html
@@ -653,8 +653,7 @@ def dry_spell_total_length(
     window: int = 3,
     op: str = "sum",
     freq: str = "YS",
-    start_date: str = "",
-    end_date: str = "",
+    **indexer,
 ) -> xarray.DataArray:
     """
     Return the total number of days in dry periods of n days and more, during which the maximum or accumulated
@@ -672,10 +671,10 @@ def dry_spell_total_length(
         Reduce operation.
     freq : str
       Resampling frequency.
-    start_date : str
-        First day of year to consider ("mm-dd").
-    end_date : str
-        Last day of year to consider ("mm-dd").
+    indexer :
+        Indexing parameters to compute the indicator on a temporal subset of the data.
+        It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`.
+        Indexing is done after finding the dry days, but before finding the spells.
 
     Returns
     -------
@@ -685,47 +684,13 @@ def dry_spell_total_length(
     pram = rate2amount(pr, out_units="mm")
     thresh = convert_units_to(thresh, pram)
 
-    pram = xarray.where(pram < (thresh if op == "max" else 0), 0, pram)
-    pram.attrs["units"] = "mm"
-
-    dry = None
-    for i in range(2):
-        pram_i = pram if dry is None else pram.sortby("time", ascending=False)
-        if op == "max":
-            mask_i = pram_i.rolling(time=window).max() < thresh
-        else:
-            mask_i = pram_i.rolling(time=window).sum() < thresh
-        dry_i = (mask_i.rolling(time=window).sum() >= 1).shift(time=-(window - 1))
-        if dry is None:
-            dry = dry_i
-        else:
-            dry_i = dry_i.sortby("time", ascending=True)
-            dt = (dry.time - dry.time[0]).dt.days
-            dry = xarray.where(dt > len(dry.time) - window - 1, dry_i, dry)
-
-    doy_start = 1
-    if start_date != "":
-        doy_start = datetime.datetime.strptime(start_date, "%m-%d").timetuple().tm_yday
-    doy_end = 365
-    if end_date != "":
-        doy_end = datetime.datetime.strptime(end_date, "%m-%d").timetuple().tm_yday
-    bisextile_year_fix = xarray.where(pram.time.dt.year % 4 == 0, 1, 0)
-    doy_start = (
-        ([doy_start] * len(pram.time))
-        if doy_start < 60
-        else (doy_start + bisextile_year_fix)
-    )
-    doy_end = (
-        ([doy_end] * len(pram.time)) if doy_end < 60 else (doy_end + bisextile_year_fix)
-    )
-    doy = xarray.where(
-        doy_end >= doy_start,
-        (pram.time.dt.dayofyear >= doy_start) & (pram.time.dt.dayofyear <= doy_end),
-        (pram.time.dt.dayofyear <= doy_end) | (pram.time.dt.dayofyear >= doy_start),
-    )
-
-    out = (dry & doy).astype(float).resample(time=freq).sum("time")
+    # Pad with NaNs i.e. assuming dry days at the end.
+    pram_pad = pram.pad(time=(0, window))
+    mask = getattr(pram_pad.rolling(time=window), op)() < thresh
+    dry = (mask.rolling(time=window).sum() >= 1).shift(time=-(window - 1))
+    dry = dry.isel(time=slice(0, pram.time.size)).astype(float)
 
+    out = select_time(dry, **indexer).resample(time=freq).sum("time")
     return to_agg_units(out, pram, "count")
 
 
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index 354443b0a..4cd090bed 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -2403,12 +2403,9 @@ def test_dry_spell(pr_series):
         + [0.01] * 3
     )
     pr_arr_3 = [3.01] * 358 + [0.99] * 14 + [3.01] * 358
-    pr_1 = pr_series(np.array(pr_arr_1))
-    pr_2 = pr_series(np.array(pr_arr_2))
-    pr_3 = pr_series(np.array(pr_arr_3), start="1981-01-01")
-    pr_1.attrs["units"] = "mm/day"
-    pr_2.attrs["units"] = "mm/day"
-    pr_3.attrs["units"] = "mm/day"
+    pr_1 = pr_series(np.array(pr_arr_1), units="mm/day")
+    pr_2 = pr_series(np.array(pr_arr_2), units="mm/day")
+    pr_3 = pr_series(np.array(pr_arr_3), start="1981-01-01", units="mm/day")
 
     for i in range(1, 4):
 
@@ -2431,17 +2428,17 @@ def test_dry_spell(pr_series):
             pr = pr_3
 
         events = xci.dry_spell_frequency(
-            pr, thresh=str(thresh) + " mm", window=window, freq="YS"
+            pr, thresh=f"{thresh} mm", window=window, freq="YS"
         )
         total_d_sum = xci.dry_spell_total_length(
             pr,
-            thresh=str(thresh * (window if i == 3 else 1)) + " mm",
+            thresh=f"{thresh * (window if i == 3 else 1)} mm",
             window=window,
             op="sum",
             freq="YS",
         )
         total_d_max = xci.dry_spell_total_length(
-            pr, thresh=str(thresh) + " mm", window=window, op="max", freq="YS"
+            pr, thresh=f"{thresh} mm", window=window, op="max", freq="YS"
         )
 
         np.testing.assert_allclose(events[0], out_events, rtol=1e-1)
@@ -2449,7 +2446,15 @@ def test_dry_spell(pr_series):
         np.testing.assert_allclose(total_d_max[0], out_total_d_max, rtol=1e-1)
 
 
-def test_dry_spell_frequency_op(self, pr_series):
+def test_dry_spell_total_length_indexer(pr_series):
+    pr = pr_series([1] * 5 + [0] * 10 + [1] * 350, start="1900-01-01", units="mm/d")
+    out = xci.dry_spell_total_length(
+        pr, window=7, op="sum", thresh="3 mm", freq="MS", date_bounds=("01-10", "12-31")
+    )
+    np.testing.assert_allclose(out, [9] + [0] * 11)
+
+
+def test_dry_spell_frequency_op(pr_series):
     pr = pr_series(
         np.array(
             [
diff --git a/xclim/testing/tests/test_precip.py b/xclim/testing/tests/test_precip.py
index 87a82edfb..8f6ea14fc 100644
--- a/xclim/testing/tests/test_precip.py
+++ b/xclim/testing/tests/test_precip.py
@@ -475,8 +475,8 @@ def test_dry_spell():
         pr, thresh="3 mm", window=7, op="max", freq="YS"
     )
 
-    total_d_sum = total_d_sum.sel(location=slice("Halifax"), drop=True).squeeze()[0:2]
-    total_d_max = total_d_max.sel(location=slice("Halifax"), drop=True).squeeze()[0:2]
+    total_d_sum = total_d_sum.sel(location="Halifax", drop=True).isel(time=slice(0, 2))
+    total_d_max = total_d_max.sel(location="Halifax", drop=True).isel(time=slice(0, 2))
 
     np.testing.assert_allclose(events[0:2, 0], [5, 8], rtol=1e-1)
     np.testing.assert_allclose(total_d_sum, [50, 60], rtol=1e-1)

From 0348d50ec5f4145498022d745739e062898bd939 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Wed, 1 Dec 2021 11:43:29 -0500
Subject: [PATCH 5/9] Add indexing tests

---
 xclim/testing/tests/test_precip.py | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/xclim/testing/tests/test_precip.py b/xclim/testing/tests/test_precip.py
index 8f6ea14fc..360b8c70d 100644
--- a/xclim/testing/tests/test_precip.py
+++ b/xclim/testing/tests/test_precip.py
@@ -464,8 +464,8 @@ def test_liquid_precip_ratio():
         assert "where temperature is above 33 degf." in out.description
 
 
-def test_dry_spell():
-    pr = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc").pr
+def test_dry_spell(atmosds):
+    pr = atmosds.pr
 
     events = atmos.dry_spell_frequency(pr, thresh="3 mm", window=7, freq="YS")
     total_d_sum = atmos.dry_spell_total_length(
@@ -496,6 +496,25 @@ def test_dry_spell():
     )
 
 
+def test_dry_spell_total_length_indexer(pr_series):
+    pr = pr_series(
+        [np.NaN] + [1] * 4 + [0] * 10 + [1] * 350, start="1900-01-01", units="mm/d"
+    )
+    out = atmos.dry_spell_total_length(
+        pr,
+        window=7,
+        op="sum",
+        thresh="3 mm",
+        freq="MS",
+    )
+    np.testing.assert_allclose(out, [np.NaN] + [0] * 11)
+
+    out = atmos.dry_spell_total_length(
+        pr, window=7, op="sum", thresh="3 mm", freq="MS", date_bounds=("01-10", "12-31")
+    )
+    np.testing.assert_allclose(out, [9] + [0] * 11)
+
+
 def test_dry_spell_frequency_op():
     pr = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc").pr
     test_sum = atmos.dry_spell_frequency(

From 9a11bf592b18667d4385e126fe2c9278c971af58 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Thu, 20 Jan 2022 17:40:48 -0500
Subject: [PATCH 6/9] Docstring ajustment

---
 xclim/indices/_agro.py | 27 ++++++++++++++++++---------
 1 file changed, 18 insertions(+), 9 deletions(-)

diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py
index 5f80e35fd..26c701e1c 100644
--- a/xclim/indices/_agro.py
+++ b/xclim/indices/_agro.py
@@ -1,6 +1,4 @@
 # noqa: D100
-
-import datetime
 from typing import Optional
 
 import numpy as np
@@ -656,8 +654,10 @@ def dry_spell_total_length(
     **indexer,
 ) -> xarray.DataArray:
     """
-    Return the total number of days in dry periods of n days and more, during which the maximum or accumulated
-    precipitation on a window of n days is under the threshold.
+    Total length of dry spells
+
+    Total number of days in dry periods of a minimum length, during which the maximum or
+    accumulated precipitation within a window of the same length is under a threshold.
 
     Parameters
     ----------
@@ -668,23 +668,32 @@ def dry_spell_total_length(
     window : int
       Number of days where the maximum or accumulated precipitation is under threshold.
     op : {"max", "sum"}
-        Reduce operation.
+      Reduce operation.
     freq : str
       Resampling frequency.
     indexer :
-        Indexing parameters to compute the indicator on a temporal subset of the data.
-        It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`.
-        Indexing is done after finding the dry days, but before finding the spells.
+      Indexing parameters to compute the indicator on a temporal subset of the data.
+      It accepts the same arguments as :py:func:`xclim.indices.generic.select_time`.
+      Indexing is done after finding the dry days, but before finding the spells.
 
     Returns
     -------
     xarray.DataArray
       The {freq} total number of days in dry periods of minimum {window} days.
+
+    Notes
+    -----
+    The algorithm assumes days before and after the timeseries are "wet", meaning that
+    the condition for being considered part of a dry spell is stricter on the edges. For
+    example, with `window=3` and `op='sum'`, the first day of the series is considered
+    part of a dry spell only if the accumulated precipitation within the first 3 days is
+    under the threshold. In comparison, a day in the middle of the series is considered
+    part of a dry spell if any of the three 3-day periods of which it is part are
+    considered dry (so a total of five days are included in the computation, compared to only 3.)
     """
     pram = rate2amount(pr, out_units="mm")
     thresh = convert_units_to(thresh, pram)
 
-    # Pad with NaNs i.e. assuming dry days at the end.
     pram_pad = pram.pad(time=(0, window))
     mask = getattr(pram_pad.rolling(time=window), op)() < thresh
     dry = (mask.rolling(time=window).sum() >= 1).shift(time=-(window - 1))

From 0d58d994ba77e6628bc7358d24b59b4baa121f6a Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Thu, 20 Jan 2022 17:52:58 -0500
Subject: [PATCH 7/9] upd hist

---
 HISTORY.rst | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/HISTORY.rst b/HISTORY.rst
index 832fab886..569289602 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -4,7 +4,7 @@ History
 
 0.33.0 (unreleased)
 -------------------
-Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Tom Keel (:user:`Thomasjkeel`), Jeremy Fyke (:user:`JeremyFyke`), David Huard (:user:`huard`), Abel Aoun (:user:`bzah`), Juliette Lavoie (:user:`juliettelavoie`).
+Contributors to this version: Trevor James Smith (:user:`Zeitsperre`), Pascal Bourgault (:user:`aulemahal`), Tom Keel (:user:`Thomasjkeel`), Jeremy Fyke (:user:`JeremyFyke`), David Huard (:user:`huard`), Abel Aoun (:user:`bzah`), Juliette Lavoie (:user:`juliettelavoie`), Yannick Rousseau (:user:`yrouranos`).
 
 Announcements
 ^^^^^^^^^^^^^
@@ -23,6 +23,7 @@ New features and enhancements
 * New ``xclim.indices.stats.parametric_cdf`` function to facilitate the computation of return periods over DataArrays of statistical distribution parameters (:issue:`876`, :pull:`984`).
 * Add ``copy`` parameter to ``percentile_doy`` to control if the array input can be dumped after computing percentiles (:issue:`932`, :pull:`985`).
 * Added ``properties.py`` and ``measures.py`` in order to perform diagnostic tests of sdba (:pull:`967`, :issue:`424`).
+* New improved algorithm for ``dry_spell_total_length``, performing the temporal indexing at the right moment and with control on the aggregation operator (``op``) for determining the dry spells.
 
 Breaking changes
 ^^^^^^^^^^^^^^^^

From 3ed1239588c6ecda909b7d557aca020a6cfef220 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Fri, 28 Jan 2022 15:08:18 -0500
Subject: [PATCH 8/9] edit indicator metadata

---
 xclim/data/fr.json                | 4 ++--
 xclim/indicators/atmos/_precip.py | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/xclim/data/fr.json b/xclim/data/fr.json
index 1c3d84675..a855612ad 100644
--- a/xclim/data/fr.json
+++ b/xclim/data/fr.json
@@ -947,8 +947,8 @@
   },
   "DRY_SPELL_TOTAL_LENGTH": {
     "title": "Durée totale en jours des périodes séches",
-    "abstract": "Durée totale en jours des périodes séches de n jours et plus, pendant lesquelles les précipitations accumulées sur une fenêtre de n jours sont inférieures à un seuil donné.",
-    "description": "Durée totale {freq:f} en jours des périodes séches de {window} jours et plus, pendant lesquelles les précipitations accumulées sur une fenêtre de {window} jours sont inférieures à {thresh}.",
+    "abstract": "Durée totale en jours des périodes séches de n jours et plus, pendant lesquelles les précipitations accumulées ou maximales sur une fenêtre de n jours sont inférieures à un seuil donné.",
+    "description": "Durée totale {freq:f} en jours des périodes séches de {window} jours et plus, pendant lesquelles les précipitations {op:fpl} sur une fenêtre de {window} jours sont inférieures à {thresh}.",
     "long_name": "Durée totale {freq:f} en jours des périodes sèches de {window} jours et plus"
   },
   "COLD_AND_DRY_DAYS": {
diff --git a/xclim/indicators/atmos/_precip.py b/xclim/indicators/atmos/_precip.py
index 7cc05206f..b1c18a5b3 100644
--- a/xclim/indicators/atmos/_precip.py
+++ b/xclim/indicators/atmos/_precip.py
@@ -370,8 +370,8 @@ class HrPrecip(Hourly):
 
 dry_spell_total_length = Precip(
     identifier="dry_spell_total_length",
-    description="The {freq} number of days in dry periods of {window} days and more, during which the accumulated "
-    "precipitation on a window of {window} days is under {thresh}.",
+    description="The {freq} number of days in dry periods of {window} days and more, during which the {op}"
+    "precipitation within windows of {window} days is under {thresh}.",
     units="days",
     cell_methods="",
     compute=indices.dry_spell_total_length,

From e19a0d685c5f47d2cdc436cf4b3b0fd1aa769e55 Mon Sep 17 00:00:00 2001
From: Pascal Bourgault <bourgault.pascal@ouranos.ca>
Date: Fri, 28 Jan 2022 15:39:49 -0500
Subject: [PATCH 9/9] parametrize tests

---
 xclim/indices/_agro.py              |   3 +-
 xclim/testing/tests/test_indices.py | 111 +++++++++++++---------------
 2 files changed, 54 insertions(+), 60 deletions(-)

diff --git a/xclim/indices/_agro.py b/xclim/indices/_agro.py
index 26c701e1c..9b670c25e 100644
--- a/xclim/indices/_agro.py
+++ b/xclim/indices/_agro.py
@@ -6,10 +6,11 @@
 
 import xclim.indices as xci
 import xclim.indices.run_length as rl
+from xclim.core.calendar import select_time
 from xclim.core.units import convert_units_to, declare_units, rate2amount, to_agg_units
 from xclim.core.utils import DayOfYearStr
 from xclim.indices._threshold import first_day_above, first_day_below, freshet_start
-from xclim.indices.generic import aggregate_between_dates, day_lengths, select_time
+from xclim.indices.generic import aggregate_between_dates, day_lengths
 
 # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start
 # See https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html
diff --git a/xclim/testing/tests/test_indices.py b/xclim/testing/tests/test_indices.py
index 985d96f20..2271bb56c 100644
--- a/xclim/testing/tests/test_indices.py
+++ b/xclim/testing/tests/test_indices.py
@@ -2471,68 +2471,61 @@ def test_water_budget(pr_series, tasmin_series, tasmax_series):
     np.testing.assert_allclose(out[1, 0], [8.5746025 / 86400], rtol=1e-1)
 
 
-def test_dry_spell(pr_series):
-
-    pr_arr_1 = (
-        [1.01] * 6
-        + [0.01] * 3
-        + [0.51] * 2
-        + [0.75] * 2
-        + [0.51]
-        + [0.01] * 3
-        + [1.01] * 3
+@pytest.mark.parametrize(
+    "pr,thresh1,thresh2,window,outs",
+    [
+        (
+            [1.01] * 6
+            + [0.01] * 3
+            + [0.51] * 2
+            + [0.75] * 2
+            + [0.51]
+            + [0.01] * 3
+            + [1.01] * 3,
+            3,
+            3,
+            7,
+            (2, 12, 20),
+        ),
+        (
+            [0.01] * 6
+            + [1.01] * 3
+            + [0.51] * 2
+            + [0.75] * 2
+            + [0.51]
+            + [0.01] * 3
+            + [0.01] * 3,
+            3,
+            3,
+            7,
+            (2, 18, 20),
+        ),
+        ([3.01] * 358 + [0.99] * 14 + [3.01] * 358, 1, 14, 14, (0, 7, 7)),
+    ],
+)
+def test_dry_spell(pr_series, pr, thresh1, thresh2, window, outs):
+
+    pr = pr_series(np.array(pr), start="1981-01-01", units="mm/day")
+
+    out_events, out_total_d_sum, out_total_d_max = outs
+
+    events = xci.dry_spell_frequency(
+        pr, thresh=f"{thresh1} mm", window=window, freq="YS"
     )
-    pr_arr_2 = (
-        [0.01] * 6
-        + [1.01] * 3
-        + [0.51] * 2
-        + [0.75] * 2
-        + [0.51]
-        + [0.01] * 3
-        + [0.01] * 3
+    total_d_sum = xci.dry_spell_total_length(
+        pr,
+        thresh=f"{thresh2} mm",
+        window=window,
+        op="sum",
+        freq="YS",
+    )
+    total_d_max = xci.dry_spell_total_length(
+        pr, thresh=f"{thresh1} mm", window=window, op="max", freq="YS"
     )
-    pr_arr_3 = [3.01] * 358 + [0.99] * 14 + [3.01] * 358
-    pr_1 = pr_series(np.array(pr_arr_1), units="mm/day")
-    pr_2 = pr_series(np.array(pr_arr_2), units="mm/day")
-    pr_3 = pr_series(np.array(pr_arr_3), start="1981-01-01", units="mm/day")
-
-    for i in range(1, 4):
-
-        if i == 1:
-            thresh, window = 3, 7
-            out_events = [2]
-            out_total_d_sum = [12]
-            out_total_d_max = [20]
-            pr = pr_1
-        elif i == 2:
-            thresh, window = 3, 7
-            out_events = [2]
-            out_total_d_sum = [18]
-            out_total_d_max = [20]
-            pr = pr_2
-        else:
-            thresh, window = 1, 14
-            out_events = [0]
-            out_total_d_sum = out_total_d_max = [7, 7]
-            pr = pr_3
-
-        events = xci.dry_spell_frequency(
-            pr, thresh=f"{thresh} mm", window=window, freq="YS"
-        )
-        total_d_sum = xci.dry_spell_total_length(
-            pr,
-            thresh=f"{thresh * (window if i == 3 else 1)} mm",
-            window=window,
-            op="sum",
-            freq="YS",
-        )
-        total_d_max = xci.dry_spell_total_length(
-            pr, thresh=f"{thresh} mm", window=window, op="max", freq="YS"
-        )
 
-        np.testing.assert_allclose(events[0], out_events, rtol=1e-1)
-        np.testing.assert_allclose(total_d_sum[0], out_total_d_sum, rtol=1e-1)
-        np.testing.assert_allclose(total_d_max[0], out_total_d_max, rtol=1e-1)
+    np.testing.assert_allclose(events[0], [out_events], rtol=1e-1)
+    np.testing.assert_allclose(total_d_sum[0], [out_total_d_sum], rtol=1e-1)
+    np.testing.assert_allclose(total_d_max[0], [out_total_d_max], rtol=1e-1)
 
 
 def test_dry_spell_total_length_indexer(pr_series):