Heat wave magnitude indicator

[seblehner]

Closes #994

Pull Request Checklist:

• This PR addresses an already opened issue (for bug fixes / features)
  • This PR fixes Heat Wave Magnitude Index #994
• Tests for the changes have been added (for bug fixes / features)
  • (If applicable) Documentation has been added / updated (for bug fixes / features)
• CHANGELOG.rst has been updated (with summary of main changes)
  • Link to issue (:issue:number) and pull request (:pull:number) has been added

What kind of change does this PR introduce?

• Adds a heat_wave_magnitude indicator based on a variable definition on heat waves (consecutive days of maximum temperature above a given threshold), which yields the highest magnitude in the form of accumulated maximum temperature differences above a threshold along a consecutive heat wave
• Adds a windowed_max_run_sum function to indices.run_length which yields the maximum accumulated values of runs per freq period
• Adds a rse (termed run sum encoding) function which sums positive values for runs

Does this PR introduce a breaking change?

No

TODO:

• one of the test doesn't convert units properly
• french translation (need help with that one)

[github-actions]

Welcome, new contributor!

It appears that this is your first Pull Request. To give credit where it's due, we ask that you add your information to the AUTHORS.rst and .zenodo.json:

• The relevant author information has been added to AUTHORS.rst and .zenodo.json

Please make sure you've read our contributing guide. We look forward to reviewing your Pull Request shortly ✨

[coxipi]

@aulemahal I'm now wondering if having the possibility of float input in rle might cause problems with rle_1d, since it also accepts floats, but it's a true run length encoding of floats, e.g.
[1.1, 1.1, 0.9, 0.9, 1]
will give couples of (value, length) as : (1.1,2) , (0.9, 2), (1,1).

I think this is really what is meant by a run length of floats, the proper way to generalize. I would not mind keeping the current generalization we have in rle, but because of the 1d case, it might cause confusion. Thoughts?

[coveralls]

coverage: 89.445% (+0.009%) from 89.436%
when pulling bfa2576 on seblehner:heat-wave-magnitude-indicator
into 0644f3a on Ouranosinc:main.

[aulemahal]

I think we could simply explain the nuance in the docstring ? Like explain what happens when a bool / int / float is passed in the input of either functions. In most cases, rle is not used by itself, but called from within one of the functions like rle_statistics or windowed_max_run_sum, and I would expect that function to choose the right function.

I can push such a commit.

[coxipi]

I think we could simply explain the nuance in the docstring ? Like explain what happens when a bool / int / float is passed in the input of either functions. In most cases, rle is not used by itself, but called from within one of the functions like rle_statistics or windowed_max_run_sum, and I would expect that function to choose the right function.

Ok, sounds good. After all it's not like the switch between 1d and nd is much more straightforward elsewhere, e.g. for rle_statistics. And although the aims of the old rle were more similar generally to rle_1d, it was still not a true rle as you point out in the doc

    ufunc_1dim = use_ufunc(ufunc_1dim, da, dim=dim, index=index, freq=freq)
    if ufunc_1dim:
        rl_stat = statistics_run_ufunc(da, reducer, window, dim)
    else:
        d = rle(da, dim=dim, index=index)

        def get_rl_stat(d):
            rl_stat = getattr(d.where(d >= window), reducer)(dim=dim)
            rl_stat = xr.where((d.isnull() | (d < window)).all(dim=dim), 0, rl_stat)
            return rl_stat

        if freq is None:
            rl_stat = get_rl_stat(d)
        else:
            rl_stat = d.resample({dim: freq}).map(get_rl_stat)

    return rl_stat

[coxipi]

That looks good to me. Thanks for this nice work!

On a final note, do you have a specific paper in mind on which you base this indicator? I had a discussion on various papers I can't find it now. The indicator in the paper cited in the original issue is a bit different, it has a standardization aspect with a calibration period, so I would not cite this paper, or close the issue yet, the implementation is a bit different to what is done here.

[seblehner]

That looks good to me. Thanks for this nice work!

Thank you (and everyone else) also for helping! Happy to contribute 🙌

On a final note, do you have a specific paper in mind on which you base this indicator? I had a discussion on various papers I can't find it now. The indicator in the paper cited in the original issue is a bit different, it has a standardization aspect with a calibration period, so I would not cite this paper, or close the issue yet, the implementation is a bit different to what is done here.

A colleague cited this paper: https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022EF002833 for the indicator. The paper is about a compound indicator for drought and heat, but the heatwave part follows the definition here (and is in turn referencing the Russo et al., 2014 paper, among others).

If I didn't miss anything, there are 2 points in particular about the Russo et al., 2014 definition:

• Reference period for threshold calculation: 1981--2010.
• Threshold calculated on a 31 day centered window.

The first point is already captured in the current implementation, because you can calculate the threshold for any reference period and supply the xr.DataArray as input for thresh (instead of specifying a static threshold like 25 degC).

The second point however can not as easily be rebuilt I think. I don't know how much functionality under the xclim hood already exists for such operations and haven't looked into that yet to be honest. I think it boils down to something comparable to the following two indicators:

• atmos.fraction_over_precip_doy_thresh
• atmos.fraction_over_precip_thresh

Those two do in principle something similar, but the former calculates the threshold over doy, while the latter can have a static one. That would essentially be the same missing part (together with the added window around doy) to exactly capture the Russo et al., 2014 indicator.

[coxipi]

Thanks! I just wanted to know which references should be cited, I've added both. I'll merge your PR as soon as the tests have re-passed

[Zeitsperre]

@seblehner Thanks so much for your work on this! xclim v0.53 will be released with your changes either this or next week. Have a good day!