Improved algorithm for index 'dry_spell_total_length' with new option…

AUTHORS.rst

@@ -31,3 +31,4 @@ Contributors
 * Jamie Quinn <jamiejquinn@jamiejquinn.com> `@JamieJQuinn <https://github.com/JamieJQuinn>`_
 * Tom Keel <thomas.keel.18@ucl.ac.uk> `@Thomasjkeel <https://github.com/Thomasjkeel>`_
 * Maliko Tanguy <malngu@ceh.ac.uk> `@malngu <https://github.com/malngu>`_
+* Yannick Rousseau `@yrouranos <https://github.com/yrouranos>`_

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
 ^^^^^^^^^^^^^
@@ -27,6 +27,7 @@ New features and enhancements
     - ``xclim.sdba.adjustment.PrincipalComponent`` was modified to have a simpler signature. The "full" method for finding the best PC orientation was added. (:issue:`697`).
 * 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`).
+* 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.
 * Added ``properties.py`` and ``measures.py`` in order to perform diagnostic tests of sdba (:issue:`424`, :pull:`967`).
 * Update how ``percentile_doy`` rechunk the input data to preserve the initial chunk size. This should make the computation memory footprint more predictable (:issue:`932`, :pull:`987`).
 

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": {

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,

xclim/indices/_agro.py

@@ -1,12 +1,12 @@
 # noqa: D100
-
 from typing import Optional
 
 import numpy as np
 import xarray
 
 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
@@ -647,11 +647,18 @@ 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",
+    **indexer,
 ) -> 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.
+    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
     ----------
@@ -660,21 +667,40 @@ 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.
+    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
     -------
     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)
 
-    mask = pram.rolling(time=window, center=True).sum() < thresh
-    out = (mask.rolling(time=window, center=True).sum() >= 1).resample(time=freq).sum()
+    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")
 
 

xclim/testing/tests/test_indices.py

@@ -2469,91 +2469,116 @@ def test_water_budget(pr_series, tasmin_series, tasmax_series):
     np.testing.assert_allclose(out[1, 0], [8.5746025 / 86400], rtol=1e-1)
 
 
-class TestDrySpell:
-    def test_dry_spell(self, 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"
+@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):
 
-        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")
+    pr = pr_series(np.array(pr), start="1981-01-01", units="mm/day")
 
-        np.testing.assert_allclose(events[0], [2], rtol=1e-1)
-        np.testing.assert_allclose(total_d[0], [12], rtol=1e-1)
+    out_events, out_total_d_sum, out_total_d_max = outs
 
-    def test_dry_spell_frequency_op(self, pr_series):
-        pr = pr_series(
-            np.array(
-                [
-                    29.012,
-                    0.1288,
-                    0.0253,
-                    0.0035,
-                    4.9147,
-                    1.4186,
-                    1.014,
-                    0.5622,
-                    0.8001,
-                    10.5823,
-                    2.8879,
-                    8.2635,
-                    0.292,
-                    0.5242,
-                    0.2426,
-                    1.3934,
-                    0.0,
-                    0.4633,
-                    0.1862,
-                    0.0034,
-                    2.4591,
-                    3.8547,
-                    3.1983,
-                    3.0442,
-                    7.422,
-                    14.8854,
-                    13.4334,
-                    0.0012,
-                    0.0782,
-                    31.2916,
-                    0.0379,
-                ]
-            )
-        )
-        pr.attrs["units"] = "mm/day"
+    events = xci.dry_spell_frequency(
+        pr, thresh=f"{thresh1} mm", window=window, freq="YS"
+    )
+    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"
+    )
 
-        test_sum = xci.dry_spell_frequency(
-            pr, thresh="1 mm", window=3, freq="MS", op="sum"
-        )
-        test_max = xci.dry_spell_frequency(
-            pr, thresh="1 mm", window=3, freq="MS", op="max"
+    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):
+    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(
+            [
+                29.012,
+                0.1288,
+                0.0253,
+                0.0035,
+                4.9147,
+                1.4186,
+                1.014,
+                0.5622,
+                0.8001,
+                10.5823,
+                2.8879,
+                8.2635,
+                0.292,
+                0.5242,
+                0.2426,
+                1.3934,
+                0.0,
+                0.4633,
+                0.1862,
+                0.0034,
+                2.4591,
+                3.8547,
+                3.1983,
+                3.0442,
+                7.422,
+                14.8854,
+                13.4334,
+                0.0012,
+                0.0782,
+                31.2916,
+                0.0379,
+            ]
         )
+    )
+    pr.attrs["units"] = "mm/day"
+
+    test_sum = xci.dry_spell_frequency(pr, thresh="1 mm", window=3, freq="MS", op="sum")
+    test_max = xci.dry_spell_frequency(pr, thresh="1 mm", window=3, freq="MS", op="max")
 
-        np.testing.assert_allclose(test_sum[0], [2], rtol=1e-1)
-        np.testing.assert_allclose(test_max[0], [3], rtol=1e-1)
+    np.testing.assert_allclose(test_sum[0], [2], rtol=1e-1)
+    np.testing.assert_allclose(test_max[0], [3], rtol=1e-1)
 
 
 class TestRPRCTot:

xclim/testing/tests/test_precip.py

@@ -464,23 +464,55 @@ 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 = 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 = 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[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 total "
         "precipitation on a window of 7 days is under 3 mm."
     ) 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
+    )
+
+
+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():