Understand differing output of tx_days_above and hot_spell_frequency with window=1

[ollie-bell]

Setup Information

• Xclim version: 0.41.0
• Python version: 3.10.9
• Operating System: Windows 10 Version 21H2

Context

Am I misunderstanding the purpose of xclim.indices.hot_spell_frequency?

Number of hot spells over a given period. A hot spell is defined as an event where the maximum daily temperature exceeds a specific threshold (default: 30℃) over a minimum number of days (default: 3).

If window=1 should it not be equivalent to tx_days_above?

I'm looking for the appropriate xclim function that will compute "The number of times in each period (according to freq resampling) that tasmax exceeds X degrees for N (or more) consecutive days". This seemed to be appropriate for hot_spell_frequency but I was surprised that when N=1 it gives different results to tx_days_above.

import xarray as xr
import numpy as np
import xclim.indices
import cftime

# generate 10 years of synthetic temperature data on a 360-day calendar.
tasmax = xr.DataArray(
    np.random.randn(360*10, 5, 5)*5 + 15,
    coords={"time": [cftime.num2date(d, "days since 1999-12-01", calendar="360_day") for d in range(360*10)]},
    dims=["time", "X", "Y"],
    attrs={"units": "degC"}
)

thresh = "20 degC"
tx_days_above = xclim.indices.tx_days_above(tasmax, thresh, freq="AS-DEC", op=">")
hot_spell_frequency = xclim.indices.hot_spell_frequency(tasmax, thresh, window=1, freq="AS-DEC", op=">")

In [120]: tx_days_above.mean(dim="time")
Out[120]:
<xarray.DataArray (X: 5, Y: 5)>
array([[58.9, 65.6, 56.3, 55. , 57.9],
       [59.5, 56.7, 60.1, 52.5, 54.7],
       [54.6, 57.4, 54. , 53.6, 54.6],
       [60.4, 54.4, 54.6, 59. , 55. ],
       [52.7, 57. , 59.2, 59.7, 57.5]])
Dimensions without coordinates: X, Y

In [121]: hot_spell_frequency.mean(dim="time")
Out[121]:
<xarray.DataArray (X: 5, Y: 5)>
array([[49.4, 54. , 47.3, 46.5, 48.6],
       [50.6, 48.8, 49.9, 45.3, 46.6],
       [46.9, 47.4, 45.1, 45.5, 46.3],
       [49.9, 46.4, 47.9, 48.8, 47.5],
       [45.2, 47.4, 49.3, 48.2, 48.2]])
Dimensions without coordinates: X, Y

Code of Conduct

• I agree to follow this project's Code of Conduct

[ollie-bell]

I think I've realised that my question boils down to:

Why does windowed_run_count not equal windowed_run_events for a window size of 1?

[cuell]

Hot spell frequency counts the number of spells, which in your case, is the number of times you get a run of 1 or more days in row. If there are 2 days in a row, then that gets counted as 1. tx_days_above is the number of days in total that the threshold is exceeded.

You'll notice in your example that _hot_spell_frequency is always less than _tx_days_above.

[aulemahal]

Hi!

It as simple as windowed_run_events counts the events of minimum window days. A spell of 5 days is counted as 1 event by this function, whereas it would be an accumulation of 5 with windowed_run_count.

That's why hot_spell_frequency would be systematically smaller or equal to tx_days_above.

In xclim, we tried to use consistent vocabulary and so wherever you see "spell" this is what it means. It's a sequence of consecutive timesteps which respects a condition and has atleast a length of window. Aggregation operation are then applied on these spells and not the timesteps they are made of. "spell_frequency" means the number of spells per period.

On the other hand days means that the aggregation operation is applied independently on timesteps and there is no minimum consecutive length of values.

[ollie-bell]

Understood! Thanks @aulemahal @cuell ...

A follow up question: Is there functionality in xclim for something like windowed_run_events but instead of counting spells greater than a certain length, it counts spells equal to a certain length?

e.g. if the window is 3 days and there is a 6-day spell, this would be computed as 2 spells back-to-back rather than 1 single spell?

[Zeitsperre]

e.g. if the window is 3 days and there is a 6-day spell, this would be computed as 2 spells back-to-back rather than 1 single spell?

If the days are continuous, this would be counted as one "spell"/event.

[coxipi]

@ollie-bell There's no indicator or function that does that exactly for you, but you could play with rl._cumsum_reset_on_zero. This function gives an intermediary step in windowed_run_events. For instance:

da = xr.DataArray([1,1,1,1,1,1,0,1,1,1,1,0])
out = rl._cumsum_reset_on_zero(da, dim="dim_0")

will give [1,2,3,4,5,6,0,1,2,3,4,0]. Normally, rl.rle and rl.windowed_run_events will clean this up.

But working with this intermediary result, you can isolate all multiples of 3 (or any desired window):

out2 = (out%3==0) & (out>0)

and count the number of times out2 is True, out2.sum(dim="dim_0"). The result is 3, i.e. 2 spells in the first sequence, and 1 spell in the last one.

With a more realistic output, you could then resample as desired out2.resample({"time":"AS-DEC"}).sum(dim="time")

[bzah]

You could also play with the resampling freq to divide your time series into 3 days long periods "freq=3D" and compute hot_spell_frequency on window=3. That only works because, by default, resampling is done before counting the spells (that's configurable).

xclim.indices.hot_spell_frequency(tasmax, thresh, window=3, freq="3D", op=">").resample(time="AS-DEC").sum()