filter_events with (resample_before_rl=False)

[kobebryant432]

Addressing a Problem?

Currently filter_events resamples the data (if freq is provided) before looking for events by default.
If using yearly frequency this can cut a valid event in half.
For other rle there is the resampel_before_rl option to account for this - although this only covers the cases in which the actual data is not required.

It would be nice to have this consistent so that windowed_run_count functionality matches the find_event.mean(dim="event").event_length even when windowed_run_count is using:

cond = xarray.DataArray(
    np.array(
        [
            [False, False, True, True, False, True, True, True, False, False],
            [True, True, True, True, True, False, False, False, False, False],
        ]
    ).T,
    dims=["time", "lat"],
    coords={
        "time": xarray.date_range("1999-12-30", periods=10, freq="D"), 
        "lat": [40.0, 41.0],
    },
)
data = xarray.DataArray(
    np.array(
        [
            [26, 26, 31, 31, 20, 26, 31, 31, 20, 20],
            [26, 26, 26, 26, 26, 20, 20, 20, 20, 20],
        ]
    ).T,
    dims=["time", "lat"],
    coords={
        "time": xarray.date_range("1999-12-30", periods=10, freq="D"), 
        "lat": [40.0, 41.0],
    },
    attrs={"units": "degC"},
)

window = 3
freq = "YS"
events = rl.find_events(
    cond,
    window=window, 
    freq=freq,
    data=data,
)
event_mean = events.mean(dim="event")

out = rl.resample_and_rl(
    cond,
    resample_before_rl=False,
    compute=rl.windowed_run_count,
    window=window, 
    freq=freq,
)

print(event_mean.event_length.values)
print(out.values)

These values are not equal and to my knowledge can currently not be aligned.
They are equal when resampe_before_rl=True (as expected).

Potential Solution

Adding a resample_before_rl option does not seem to easily work currently.
Any ideas on how to integrate this?

Additional context

No response

Contribution

• I would be willing/able to open a Pull Request to contribute this feature.

Code of Conduct

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

[kobebryant432]

I guess a neat way to avoid this is using freq="YS-JUL" for event types that happen mainly in the winter (assuming Northern Hemisphere - an assumption that is implicitly made also in the cold_spell_... indicators as this is the default here). But if events can happen year round this is not a solution.

[aulemahal]

Hi!

At first I thought it would be a simple resample call, but unless I am missing something, there's no way to do something like resample(event_start=FREQ).construct('time') similarly to a rolling(time=N).construct(). Moreover, you can't resample on a data variable, only on coordinates!

I tried to implement something in branch resample-find-events, here:

xclim/src/xclim/indices/run_length.py

Line 1750 in 094097a

def group_events_by_period(events: xr.Dataset, freq: str):

The function uses the event_start min and max to generate a time coordinate at the specified frequency an groups (reshapes) the events dataset to divide the events by period of their start.

I don't think there's an easy way to compute this lazily as we need the min and max of the event_start. For now, the function will raise upon receiving a dask array. To implement this lazily, I think more acrobatics are needed and the function would need the original time axis.

Moreover, as the original time axis is not used, the final periods can be different from the resample_before_rl=True case if a year has no events. Maybe that's problematic ?

[aulemahal]

Oh and as it doesn't support dask, I included dropping the extra space along the event dimension.

It is usable as:

import xclim as xc
from xclim.testing import open_dataset

ds = open_dataset("ERA5/daily_surface_cancities_1990-1993.nc")

events_before = xc.indices.run_length.find_events(ds.tas > 290, 3, data=ds.tas, freq='YS', resample_before_rl=True)
events_after = xc.indices.run_length.find_events(ds.tas > 290, 3, data=ds.tas, freq='YS', resample_before_rl=False)

print(events_before.dims)
# {'time': 4, 'location': 5, 'event': 92}

print(events_after.dims)
# {'time': 3, 'location': 5, 'event': 18}

[kobebryant432]

Thanks for looking into this, indeed it felt like it should be rather easy but turned out to be a bit messy.

For my application using freq="YS-JUL" is the simpler solution, which I will use. I am using dask as well so, although this is a neat solution, it wouldn't cover the potential need.

Moreover, as the original time axis is not used, the final periods can be different from the resample_before_rl=True case if a year has no events. Maybe that's problematic ?

I would say this is rather problematic yes, intuitively I would expect to have as many years irrespective of the resampling method as the freq is the same.

Overall, I do think that this is a rather rare case as most of these cases could be avoided by centering the year differently.