Heat spells (#1885)

<!--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
- [x] Tests for the changes have been added (for bug fixes / features)
- [x] (If applicable) Documentation has been added / updated (for bug
fixes / features)
- [x] CHANGELOG.rst has been updated (with summary of main changes)
- [x] Link to issue (:issue:`number`) and pull request (:pull:`number`)
has been added

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

* 3 new indicators `heat_spell_frequency`, `heat_spell_max_length` and
`heat_spell_total_length`. The "heat_wave" indicators test that a
bivariate condition on each day and then find consecutive runs of a
minimum length. The new functions adhere to xclim's concept of "spell" :
the condition is tested against a rolling N-day window and if it is
fulfilled *all days* of the window are considered part of the spell.
* New generic function `spell_mask` that encapsulate the spell notion.
It returns a boolean mask of the days part of spells. The code has been
added two features:
+ Multivariate : input data can be a list (and so must be the
thresholds). The way we check for the conditions across variables is
decided by new argument `var_reducer` ('all' or 'any').
+ Weights : if the spell stat is a mean, weights can be given (used by
the INSPQ, but not for the actual PC indicators).
* Indicator parameters can now be assigned new names. See examples in
how the heat spell indicators are created. We could already do something
similar with data variables. The internals of `Indicator` have changed
accordingly (sorry).


### Does this PR introduce a breaking change?
No.

### Other information:
I kinda decided that "heat" meant "tasmax and tasmin".

CHANGELOG.rst

@@ -4,7 +4,16 @@ Changelog
 
 v0.53.0 (unreleased)
 --------------------
-Contributors to this version: Adrien Lamarche (:user:`LamAdr`), Trevor James Smith (:user:`Zeitsperre`), Éric Dupuis (:user:`coxipi`).
+Contributors to this version: Adrien Lamarche (:user:`LamAdr`), Trevor James Smith (:user:`Zeitsperre`),  Éric Dupuis (:user:`coxipi`), Pascal Bourgault (:user:`aulemahal`).
+
+New indicators
+^^^^^^^^^^^^^^
+* New ``heat_spell_frequency``, ``heat_spell_max_length`` and ``heat_spell_total_length`` : spell length statistics on a bivariate condition that uses the average over a window by default. (:pull:`1885`).
+
+New features and enhancements
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+* New generic ``xclim.indices.generic.spell_mask``  that returns a mask of which days are part of a spell. Supports multivariate conditions and weights. Used in new generic index ``xclim.indices.generic.bivariate_spell_length_statistics`` that extends ``spell_length_statistics`` to two variables.  (:pull:`1885`).
+* Indicator parameters can now be assigned a new name, different from the argument name in the compute function. (:pull:`1885`).
 
 Bug fixes
 ^^^^^^^^^
@@ -32,6 +41,7 @@ Internal changes
 * Many ``DeprecationWarning`` and ``FutureWarning`` messages emitted from `xarray` and `pint` have been addressed. (:issue:`1719`, :pull:`1881`).
 * The codebase has been adjusted to address many `pylint`-related warnings and errors. In some cases, `casting` was used to redefine some `numpy` and `xarray` objects. (:issue:`1719`, :pull:`1881`).
 * ``xclim.core`` now uses absolute imports for clarity and some objects commonly used in the module have been moved to hidden submodules. (:issue:`1719`, :pull:`1881`).
+* ``xclim.core.indicator.Parameter`` has a new attribute ``compute_name`` while ``xclim.core.indicator.Indicator`` lost its ``_variable_mapping``. The translation from parameter (and variable) names in the indicator to the names on the compute function is handled by ``Indicator._get_compute_args``. (:pull:`1885`).
 
 v0.52.0 (2024-08-08)
 --------------------

tests/test_formatting.py

@@ -35,7 +35,7 @@ def test_indicator_docstring():
     assert (
         doc[41]
         == "    Total number of series of at least {window} consecutive days with daily minimum temperature above "
-        "{thresh_tasmin} and daily maximum temperature above {thresh_tasmax} (heat_wave_events), "
+        "{thresh_tasmin} and daily maximum temperature above {thresh_tasmax}, "
         "with additional attributes: **description**: {freq} number of heat wave events within a given period. "
         "A heat wave occurs when daily minimum and maximum temperatures exceed {thresh_tasmin} and {thresh_tasmax}, "
         "respectively, over at least {window} days."

tests/test_generic.py

@@ -663,3 +663,108 @@ def test_select_time_errors(self):
 
         with pytest.raises(TypeError):
             select_time(da, doy_bounds=(300, 203, 202))
+
+
+class TestSpellMask:
+    def test_single_variable(self):
+        data = xr.DataArray([0, 1, 2, 3, 2, 1, 0, 0], dims=("time",))
+
+        out = generic.spell_mask(data, 3, "min", ">=", 2)
+        np.testing.assert_array_equal(
+            out, np.array([0, 0, 1, 1, 1, 0, 0, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask(data, 3, "max", ">=", 2)
+        np.testing.assert_array_equal(
+            out, np.array([1, 1, 1, 1, 1, 1, 1, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask(data, 2, "mean", ">=", 2)
+        np.testing.assert_array_equal(
+            out, np.array([0, 0, 1, 1, 1, 0, 0, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask(data, 3, "mean", ">", 2, weights=[0.2, 0.4, 0.4])
+        np.testing.assert_array_equal(
+            out, np.array([0, 1, 1, 1, 1, 0, 0, 0]).astype(bool)
+        )
+
+    def test_multiple_variables(self):
+        data1 = xr.DataArray([0, 1, 2, 3, 2, 1, 0, 0], dims=("time",))
+        data2 = xr.DataArray([1, 2, 3, 2, 1, 0, 0, 0], dims=("time",))
+
+        out = generic.spell_mask([data1, data2], 3, "min", ">=", [2, 2])
+        np.testing.assert_array_equal(
+            out, np.array([0, 0, 0, 0, 0, 0, 0, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask(
+            [data1, data2], 3, "min", ">=", [2, 2], var_reducer="any"
+        )
+        np.testing.assert_array_equal(
+            out, np.array([0, 1, 1, 1, 1, 0, 0, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask([data1, data2], 2, "mean", ">=", [2, 2])
+        np.testing.assert_array_equal(
+            out, np.array([0, 0, 1, 1, 0, 0, 0, 0]).astype(bool)
+        )
+
+        out = generic.spell_mask(
+            [data1, data2], 3, "mean", ">", [2, 1.5], weights=[0.2, 0.4, 0.4]
+        )
+        np.testing.assert_array_equal(
+            out, np.array([0, 1, 1, 1, 1, 0, 0, 0]).astype(bool)
+        )
+
+    def test_errors(self):
+        data = xr.DataArray([0, 1, 2, 3, 2, 1, 0, 0], dims=("time",))
+
+        # Threshold must be seq
+        with pytest.raises(ValueError, match="must be a sequence of the same length"):
+            generic.spell_mask([data, data], 3, "min", "<=", 2)
+
+        # Threshold must be same length
+        with pytest.raises(ValueError, match="must be a sequence of the same length"):
+            generic.spell_mask([data, data], 3, "min", "<=", [2])
+
+        # Weights must have win_reducer = 'mean'
+        with pytest.raises(
+            ValueError, match="is only supported if 'win_reducer' is 'mean'"
+        ):
+            generic.spell_mask(data, 3, "min", "<=", 2, weights=[1, 2, 3])
+
+        # Weights must have same length as window
+        with pytest.raises(ValueError, match="Weights have a different length"):
+            generic.spell_mask(data, 3, "mean", "<=", 2, weights=[1, 2])
+
+
+def test_spell_length_statistics_multi(tasmin_series, tasmax_series):
+    tn = tasmin_series(
+        np.zeros(
+            365,
+        )
+        + 270,
+        start="2001-01-01",
+    )
+    tx = tasmax_series(
+        np.zeros(
+            365,
+        )
+        + 270,
+        start="2001-01-01",
+    )
+
+    outc, outs, outm = generic.bivariate_spell_length_statistics(
+        tn,
+        "0 °C",
+        tx,
+        "1°C",
+        window=5,
+        win_reducer="min",
+        op="<",
+        spell_reducer=["count", "sum", "max"],
+        freq="YS",
+    )
+    xr.testing.assert_equal(outs, outm)
+    np.testing.assert_allclose(outc, 1)

tests/test_indicators.py

@@ -709,8 +709,7 @@ def test_indicator_from_dict():
     assert ind.parameters["threshold"].description == "A threshold temp"
     # Injection of parameters
     assert ind.injected_parameters["op"] == "<"
-    # Default value for input variable injected and meta injected
-    assert ind._variable_mapping["data"] == "tas"
+    assert ind.parameters["tas"].compute_name == "data"
     assert signature(ind).parameters["tas"].default == "tas"
     assert ind.parameters["tas"].units == "[temperature]"
 

tests/test_temperature.py

@@ -610,6 +610,100 @@ def test_3d_data_with_nans(self, open_dataset):
         assert np.isnan(gdd.values[0, -1, -1])
 
 
+class TestHeatSpellFrequency:
+    def test_1d(self, tasmax_series, tasmin_series):
+        tn1 = np.zeros(366)
+        tx1 = np.zeros(366)
+        tn1[:10] = np.array([20, 23, 23, 23, 20, 20, 23, 23, 23, 23])
+        tx1[:10] = np.array([29, 31, 31, 31, 28, 28, 31, 31, 31, 31])
+
+        tn = tasmin_series(tn1 + K2C, start="1/1/2000")
+        tx = tasmax_series(tx1 + K2C, start="1/1/2000")
+
+        hsf = atmos.heat_spell_frequency(
+            tn,
+            tx,
+            thresh_tasmin="22.1 C",
+            thresh_tasmax="30.1 C",
+            freq="YS",
+        )
+        np.testing.assert_allclose(hsf.values[:1], 2)
+
+        hsf = atmos.heat_spell_frequency(
+            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C", window=5, freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 1)
+
+        # no hs
+        hsf = atmos.heat_spell_frequency(
+            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C", freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 0)
+
+
+class TestHeatSpellMaxLength:
+    def test_1d(self, tasmax_series, tasmin_series):
+        tn1 = np.zeros(366)
+        tx1 = np.zeros(366)
+        tn1[:10] = np.array([20, 23, 23, 23, 20, 20, 23, 23, 23, 23])
+        tx1[:10] = np.array([29, 31, 31, 31, 28, 28, 31, 31, 31, 31])
+
+        tn = tasmin_series(tn1 + K2C, start="1/1/2000")
+        tx = tasmax_series(tx1 + K2C, start="1/1/2000")
+
+        hsf = atmos.heat_spell_max_length(
+            tn,
+            tx,
+            thresh_tasmin="22.1 C",
+            thresh_tasmax="30.1 C",
+            freq="YS",
+        )
+        np.testing.assert_allclose(hsf.values[:1], 4)
+
+        hsf = atmos.heat_spell_max_length(
+            tn,
+            tx,
+            thresh_tasmin="22 C",
+            thresh_tasmax="30 C",
+            window=5,
+            freq="YS",
+        )
+        np.testing.assert_allclose(hsf.values[:1], 5)
+
+        # no hs
+        hsf = atmos.heat_spell_max_length(
+            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C", freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 0)
+
+
+class TestHeatSpellTotalLength:
+    def test_1d(self, tasmax_series, tasmin_series):
+        tn1 = np.zeros(366)
+        tx1 = np.zeros(366)
+        tn1[:10] = np.array([20, 23, 23, 23, 20, 20, 23, 23, 23, 23])
+        tx1[:10] = np.array([29, 31, 31, 31, 28, 28, 31, 31, 31, 31])
+
+        tn = tasmin_series(tn1 + K2C, start="1/1/2000")
+        tx = tasmax_series(tx1 + K2C, start="1/1/2000")
+
+        hsf = atmos.heat_spell_total_length(
+            tn, tx, thresh_tasmin="22.1 C", thresh_tasmax="30.1 C", freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 7)
+
+        hsf = atmos.heat_spell_total_length(
+            tn, tx, thresh_tasmin="22 C", thresh_tasmax="30 C", window=5, freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 5)
+
+        # no hs
+        hsf = atmos.heat_spell_total_length(
+            tn, tx, thresh_tasmin="40 C", thresh_tasmax="40 C", freq="YS"
+        )
+        np.testing.assert_allclose(hsf.values[:1], 0)
+
+
 class TestHeatWaveFrequency:
     def test_1d(self, tasmax_series, tasmin_series):
         tn1 = np.zeros(366)