Heat wave magnitude indicator

.zenodo.json

@@ -148,6 +148,11 @@
       "name": "Wang, Hui-Min",
       "affiliation": "National University of Singapore, Singapore, Singapore",
       "orcid": "0000-0002-5878-7542"
+    },
+    {
+      "name": "Lehner, Sebastian",
+      "affiliation": "GeoSphere Austria, Vienna, Austria",
+      "orcid": "0000-0002-7562-8172"
     }
   ],
   "keywords": [

AUTHORS.rst

@@ -46,3 +46,4 @@ Contributors
 * Javier Diez-Sierra <[email protected]> `@JavierDiezSierra <https://github.com/JavierDiezSierra>`_
 * Hui-Min Wang `@Hem-W <https://github.com/Hem-W>`_
 * Adrien Lamarche `@LamAdr <https://github.com/LamAdr>`_
+* Sebastian Lehner <[email protected]> `@seblehner <https://github.com/seblehner>`_

CHANGELOG.rst

@@ -14,11 +14,13 @@ Announcements
 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 ``heat_wave_magnitude`` : yields the magnitude of the most intensive heat wave. (:pull:`1926`).
 
 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`).
+* ``xclim.indices.run_length.windowed_max_run_sum`` accumulates positive values across runs and yields the the maximum valued run. (:pull:`1926`).
 
 Bug fixes
 ^^^^^^^^^

tests/test_indices.py

@@ -1303,6 +1303,18 @@ def test_simple(self, tasmax_series):
         np.testing.assert_array_equal(out, [10, 0, 12, 8, 0, 0, 0, 0, 0, 0, 0, 0])
 
 
+class TestHeatWaveMagnitude:
+    def test_simple(self, tasmax_series):
+        a = np.zeros(365)
+        a[15:20] += 30  # 5 days
+        a[40:42] += 50  # too short -> 0
+        a[86:96] += 30  # at the end and beginning
+        da = tasmax_series(a + K2C)
+
+        out = xci.heat_wave_magnitude(da, thresh="25 C", freq="ME")
+        np.testing.assert_array_equal(out, [25, 0, 30, 20, 0, 0, 0, 0, 0, 0, 0, 0])
+
+
 class TestHeatWaveFrequency:
     @pytest.mark.parametrize(
         "thresh_tasmin,thresh_tasmax,window,expected",

tests/test_run_length.py

@@ -124,6 +124,35 @@ def test_rle(ufunc, use_dask, index):
         np.testing.assert_array_equal(out, expected)
 
 
+@pytest.mark.parametrize("use_dask", [True, False])
+@pytest.mark.parametrize("index", ["first", "last"])
+def test_rse(ufunc, use_dask, index):
+    if use_dask and ufunc:
+        pytest.xfail("rse_1d is not implemented for dask arrays.")
+
+    values = np.zeros((10, 365, 4, 4))
+    time = pd.date_range("2000-01-01", periods=365, freq="D")
+    values[:, 1:11, ...] = 30
+    da = xr.DataArray(values, coords={"time": time}, dims=("a", "time", "b", "c"))
+
+    if ufunc:
+        pytest.xfail("rse_1d is not implemented.")
+    else:
+        if use_dask:
+            da = da.chunk({"a": 1, "b": 2})
+
+        out = rl.rse(da, index=index).mean(["a", "b", "c"])
+        if index == "last":
+            expected = np.zeros(365)
+            expected[1:10] = np.nan
+            expected[10] = 300
+        else:
+            expected = np.zeros(365)
+            expected[1] = 300
+            expected[2:11] = np.nan
+        np.testing.assert_array_equal(out, expected)
+
+
 @pytest.mark.parametrize("use_dask", [True, False])
 @pytest.mark.parametrize("index", ["first", "last"])
 def test_extract_events_identity(use_dask, index):
@@ -370,6 +399,16 @@ def test_simple(self, index):
         ) + len(a[34:45])
 
 
+class TestWindowedMaxRunSum:
+    @pytest.mark.parametrize("index", ["first", "last"])
+    def test_simple(self, index):
+        a = xr.DataArray(np.zeros(50, float), dims=("time",))
+        a[4:6] = 5  # too short
+        a[25:30] = 5  # long enough, but not max
+        a[35:45] = 5  # max sum => yields 10*5
+        assert rl.windowed_max_run_sum(a, 3, dim="time", index=index) == 50
+
+
 class TestLastRun:
     @pytest.mark.parametrize(
         "coord,expected",

tests/test_temperature.py

@@ -889,6 +889,46 @@ def test_nan_presence(self, tasmax_series):
         np.testing.assert_array_equal(hwi, [np.nan])
 
 
+class TestHeatWaveMagnitude:
+    def test_simple(self, tasmax_series):
+        tx = np.zeros(366)
+        tx[:5] = np.array([30, 30, 30, 30, 30])
+        tx = tasmax_series(tx + K2C, start="1/1/2000")
+        hwm = atmos.heat_wave_magnitude(tx, freq="YS")
+        np.testing.assert_array_equal(hwm, [25])
+
+    def test_small_window_single_day(self, tasmax_series):
+        tx = np.zeros(366)
+        tx[5:8] = np.array([30, 0, 30])
+        tx = tasmax_series(tx + K2C, start="1/1/2000")
+        hwm = atmos.heat_wave_magnitude(tx, window=1, freq="YS")
+        np.testing.assert_array_equal(hwm, [5])
+
+    def test_small_window_double_day(self, tasmax_series):
+        tx = np.zeros(366)
+        tx[5:7] = np.array([30, 30])
+        tx = tasmax_series(tx + K2C, start="1/1/2000")
+        hwm = atmos.heat_wave_magnitude(tx, window=1, freq="YS")
+        np.testing.assert_array_equal(hwm, [10])
+
+    def test_convert_units(self, tasmax_series):
+        tx = np.zeros(366)
+        tx[:5] = np.array([30, 30, 30, 30, 30])
+        tx = tasmax_series(tx, start="1/1/2000")
+        tx.attrs["units"] = "C"
+        hwm = atmos.heat_wave_magnitude(tx, freq="YS")
+        np.testing.assert_array_equal(hwm, [25])
+
+    def test_nan_presence(self, tasmax_series):
+        tx = np.zeros(366)
+        tx[:5] = np.array([30, 30, 30, 30, 30])
+        tx[-1] = np.nan
+        tx = tasmax_series(tx + K2C, start="1/1/2000")
+
+        hwm = atmos.heat_wave_magnitude(tx, freq="YS")
+        np.testing.assert_array_equal(hwm, [np.nan])
+
+
 class TestDailyFreezeThaw:
     nc_tasmax = os.path.join("NRCANdaily", "nrcan_canada_daily_tasmax_1990.nc")
     nc_tasmin = os.path.join("NRCANdaily", "nrcan_canada_daily_tasmin_1990.nc")

xclim/data/fr.json

@@ -377,6 +377,18 @@
     "title": "Indice de vague de chaleur",
     "abstract": "Nombre de jours de vagues de chaleur. Une vague de chaleur se produit lorsque la température maximale quotidienne excède un seuil donné durant un certain nombre de jours."
   },
+  "HEAT_WAVE_MAGNITUDE": {
+    "long_name": "Différence cumulative maximale entre la température maximale quotidienne et {thresh} "
+    "pour les jours faisant partie d'une vague de chaleur. Une vague de chaleur est définie comme une série d'au moins "
+    "{window} jours consécutifs où la température maximale quotidien est au-dessus de {thresh}",
+    "description": "Cumul {freq:m} maximal de la différence entre la température et {thresh} pour les jours faisant partie "
+    "d'une vague de chaleur. Une vague de chaleur est définie comme une série d'au moins {window} jours consécutifs où la "
+    "température maximale quotidien est au-dessus de {thresh}",
+    "title": "Magnitude de la vague de chaleur {freq:m} la plus intense. Une vague de chaleur est définie comme une série d'au moins "
+    "{window} jours consécutifs où la température maximale quotidien est au-dessus de {thresh}.",
+    "abstract": "Magnitude de la vague de chaleur {freq:m} la plus intense. Une vague de chaleur est définie comme une série d'au moins "
+    "{window} jours consécutifs où la température maximale quotidien est au-dessus de {thresh}.",
+  },
   "TG": {
     "long_name": "Moyenne des températures maximale et minimale quotidennes",
     "description": "Température moyenne quotidienne moyenne estimée par la moyenne des températures maximale et minimale quotidiennnes.",

xclim/indicators/atmos/_temperature.py

@@ -50,6 +50,7 @@
     "heat_spell_total_length",
     "heat_wave_frequency",
     "heat_wave_index",
+    "heat_wave_magnitude",
     "heat_wave_max_length",
     "heat_wave_total_length",
     "heating_degree_days",
@@ -215,6 +216,20 @@ class TempWithIndexing(ResamplingIndicatorWithIndexing):
     compute=indices.heat_wave_frequency,
 )
 
+heat_wave_magnitude = Temp(
+    title="Heat wave magnitude",
+    identifier="heat_wave_magnitude",
+    units="K d",
+    long_name="Maximum cumulative difference between daily maximum temperature and {thresh} for days within a heat wave. "
+    "A heat wave is defined as a series of at least {window} consecutive days with daily maximum temperature above "
+    "{thresh}.",
+ description="Magnitude of the most intensive heat wave per {freq}. The magnitude is the cumulative exceedance of daily maximum temperature over {thresh}. A heat wave is defined as a series of at least {window} consecutive days with daily maximum temperature above"
+    abstract="Magnitude of the most intensive heat wave per {freq}. A heat wave occurs when daily maximum "
+    "temperatures exceed given thresholds for a number of days.",
+    cell_methods="",
+    compute=indices.heat_wave_magnitude,
+)
+
 heat_wave_max_length = Temp(
     title="Heat wave maximum length",
     identifier="heat_wave_max_length",

xclim/indices/_threshold.py

@@ -64,6 +64,7 @@
     "growing_season_length",
     "growing_season_start",
     "heat_wave_index",
+    "heat_wave_magnitude",
     "heating_degree_days",
     "hot_spell_frequency",
     "hot_spell_max_length",
@@ -1928,6 +1929,55 @@ def heat_wave_index(
     return to_agg_units(out, tasmax, "count")
 
 
+@declare_units(tasmax="[temperature]", thresh="[temperature]")
+def heat_wave_magnitude(
+    tasmax: xarray.DataArray,
+    thresh: Quantified = "25.0 degC",
+    window: int = 3,
+    freq: str = "YS",
+    op: str = ">",
+    resample_before_rl: bool = True,
+) -> xarray.DataArray:
+    """Heat wave magnitude index.
+
+    Magnitude of the most intensive heat wave event as sum of differences between tasmax
+    and the given threshold for Heat Wave days, defined as three or more consecutive days
+    over the threshold.
+
+    Parameters
+    ----------
+    tasmax : xarray.DataArray
+        Maximum daily temperature.
+    thresh : xarray.DataArray
+        Threshold temperature on which to designate a heatwave.
+    window : int
+        Minimum number of days with temperature above threshold to qualify as a heatwave.
+    freq : str
+        Resampling frequency.
+    op : {">", ">=", "gt", "ge"}
+        Comparison operation. Default: ">".
+    resample_before_rl : bool
+        Determines if the resampling should take place before or after the run
+        length encoding (or a similar algorithm) is applied to runs.
+
+    Returns
+    -------
+    DataArray, [time]
+        Heat wave magnitude.
+    """
+    thresh = convert_units_to(thresh, tasmax)
+    over_values = (tasmax - thresh).clip(0)
+
+    out = rl.resample_and_rl(
+        over_values,
+        resample_before_rl,
+        rl.windowed_max_run_sum,
+        window=window,
+        freq=freq,
+    )
+    return to_agg_units(out, tasmax, op="integral")
+
+
 @declare_units(tas="[temperature]", thresh="[temperature]")
 def heating_degree_days(
     tas: xarray.DataArray,

xclim/indices/run_length.py

@@ -180,7 +180,8 @@ def rle(
     xr.DataArray
         Values are 0 where da is False (out of runs).
     """
-    da = da.astype(int)
+    if da.dtype == bool:
+        da = da.astype(int)
 
     # "first" case: Algorithm is applied on inverted array and output is inverted back
     if index == "first":
@@ -192,7 +193,7 @@ def rle(
     # Keep total length of each series (and also keep 0's), e.g. 100120123 -> 100N20NN3
     # Keep numbers with a 0 to the right and also the last number
     cs_s = cs_s.where(da.shift({dim: -1}, fill_value=0) == 0)
-    out = cs_s.where(da == 1, 0)  # Reinsert 0's at their original place
+    out = cs_s.where(da > 0, 0)  # Reinsert 0's at their original place
 
     # Inverting back if needed e.g. 100N20NN3 -> 3NN02N001. This is the output of
     # `rle` for 111011001 with index == "first"
@@ -409,6 +410,56 @@ def windowed_run_count(
     return out
 
 
+def windowed_max_run_sum(
+    da: xr.DataArray,
+    window: int,
+    dim: str = "time",
+    freq: str | None = None,
+    ufunc_1dim: str | bool = "from_context",
+    index: str = "first",
+) -> xr.DataArray:
+    """Return the maximum sum of consecutive float values for runs at least as long as the given window length.
+
+    Parameters
+    ----------
+    da : xr.DataArray
+        Input N-dimensional DataArray (boolean).
+    window : int
+        Minimum run length.
+        When equal to 1, an optimized version of the algorithm is used.
+    dim : str
+      Dimension along which to calculate consecutive run (default: 'time').
+    freq : str
+      Resampling frequency.
+    ufunc_1dim : Union[str, bool]
+        Use the 1d 'ufunc' version of this function : default (auto) will attempt to select optimal
+        usage based on number of data points. Using 1D_ufunc=True is typically more efficient
+        for DataArray with a small number of grid points.
+        Ignored when `window=1`. It can be modified globally through the "run_length_ufunc" global option.
+    index : {'first', 'last'}
+        If 'first', the run length is indexed with the first element in the run.
+        If 'last', with the last element in the run.
+
+    Returns
+    -------
+    xr.DataArray, [int]
+        Total number of `True` values part of a consecutive runs of at least `window` long.
+    """
+    if window == 1 and freq is None:
+        out = rle(da, dim=dim, index=index).max(dim=dim)
+
+    else:
+        d_rse = rle(da, dim=dim, index=index)
+        d_rle = rle((da > 0).astype(bool), dim=dim, index=index)
+
+        d = d_rse.where(d_rle >= window, 0)
+        if freq is not None:
+            d = d.resample({dim: freq})
+        out = d.max(dim=dim)
+
+    return out
+
+
 def _boundary_run(
     da: xr.DataArray,
     window: int,