Heat Wave Magnitude Index

[ghost]

I would like to see the following indice considered as potential contributions to the xclim library:
The Heat Wave Magnitude Index is defined as the maximum magnitude of the heat waves in a year, where heat wave is the period ≥ 3 consecutive days with maximum temperature above the daily threshold for the reference period (https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2014JD022098/)

[huard]

Hi @ninglk,

You're very welcome to submit a PR for this. Please check the documentation on how to create new indices and indicators here: https://xclim.readthedocs.io/en/stable/contributing.html#implement-features-indices-or-indicators

I suspect this is not going to be straightforward, so don't hesitate to reach out for in-depth discussion on the implementation. Don't worry about the french translation, we'll complete that part.

[Zeitsperre]

Hi @ninglk, we're just wondering if you were still interested in contributing a PR to add this indice into the library. The Xclim devs are available if you need help with formatting standards and compliance or getting things working. Please let us know!

[seblehner]

Hey @Zeitsperre @huard,

I came across this issue, because I was also in need of a Heatwave Magnitude Index. I have implemented it for myself using some functions from available xclim code that I modified.

To briefly elaborate what I did:

• I essentially started with xit.atmos.heat_wave_index and wrote a heat_wave_magnitude function that is structured similarly. The main difference is that for the magnitude, we need a summation over float values instead of a binary count.
• I wrote a function windowed_run_sum (for the lack of a better name) based on windowed_run_count from the xit.indices.run_length module that sums values instead of counting binary values. (I didn't spend time writing a ufunc version yet, since I didn't need that)
• Lastly, I wrote a rse (run summation encoding; again for the lack of my imagination) based on rle (run length encoding) and similarly modified it such that float values are summed over instead of a length count for runs.

Hence, I wanted to ask if this would still be in question for a PR.

Kind regards

[Zeitsperre]

Hi @seblehner, thanks for chiming in here. We're always motivated to hear of researchers using our libraries!

It sounds like your approach would be compatible with the library, perhaps even more so given some of the changes we've made recently to our generic indices. I'm not the expert on the vectorization code, but I imagine we (@aulemahal) would be more than capable to help you with that step.

We would more than welcome a PR for this. Feel free to re-open this issue and send along a PR that we can take a look at.

All the best,

[coxipi]

Hi! Thanks for this discussion. If you want to share your code I would be interested to see how you implemented this.

Am I understanding this is similar to a function like degree_days_above? For days with tasmax > thresh, you sum the temperatures (with the additional requirement that you must be in a heat wave)? Or something like this?

We have a PR #1778 currently implementing something that could be of use for you I believe. The context in this case is that we compute freeze rain events, and we want various statistics, including: in each spell, what is the total quantity of precipitation.

We simply use intermediate results of rle where the runs are identified, then for data points in a run, we sum the value of the dataset:

ds["event_accumulation"] = rl._cumsum_reset(
        dataam.where(runs == 1), index="first", reset_on_zero=False
    )

This gives the accumulation (of freeze rain) in each spell. You can resample this result at the end if you would want the total temperature in the heat waves

There are some specificities for the freezing rain, but I think all the logic is already there to compute what you seek.

Best

[seblehner]

Hello @coxipi,

I am happy to share it already. It sounds like what we want is similar. What I did was

1. modifying windowed_run_count (the else block part, which was relevant for me) like so:

old:

d = rle(da, dim=dim, index=index)
d = d.where(d >= window, 0)
if freq is not None:
    d = d.resample({dim: freq})
out = d.sum(dim=dim)

new (inside a new windowed_run_sum function:

d_rse = rse(da, dim=dim, index=index)
d_rle = rle(xr.where(da > 0, 1, 0), 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)

The rle and where calls are effectively the same, the additional xr.where in the argument is just to generate the binary array, since the input arg is not binary, because we need the absolute values of course. rse is a modified rle function that I'll explain in a second. The max call at the end is hardcoded and not yet generally applicable. I needed only the accumulated values from the most intensive event per resampled period. This would need adaption to cover other cases as well.

The rse function:

old (from rle):

da = da.astype(int)

...

# 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

da = da.astype(float) # changed from int to float because we need absolute values

...

# 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 > 0, 0)  # Reinsert 0's at their original place

For this I just changed the type and the one where call at the end. instead of keeping the ones, which works on a binary array, I needed to keep all positive values. This where call would perhaps also need a slight adaption to make it general, e.g. if you would want to do the same for values below a threshold.

And that's essentially it. As mentioned above, I made a new indicator heat_wave_magnitude that is structured similar to the existing heat_wave_index indicator and just adapted the functions that are called.

@Zeitsperre Just fyi, I don't have permissions to reopen this issue.

Kind regards

[seblehner]

Hey @coxipi,

I checked the current PR #1778 you mentioned, but since the implemented helper function does something different (as far as I understood it), I went ahead and opened PR #1926 with the changes I mentioned above.