Computing 'seasonal means' spanning 4 months (with resample or groupy)

[oliviermarti]

Climatologists often use 'seasonal means', i.e. means over 3 months. Useful periods are DJF for December-January-February, MAM, JJA and SON. groupby or resample are nice functions to compute these seasonal means.

sea for example :
https://xarray.pydata.org/en/stable/examples/monthly-means.html
https://stackoverflow.com/questions/59234745/is-there-any-easy-way-to-compute-seasonal-mean-with-xarray

But some studies need means over 4 months : DJFM, MAMJ, JJAS and SOND. Would it be feasible that these 4-month periods are recognized by groupby and resample ?

For resample, we define 3-month means with a syntax like : resample(time='QS-DEC')

A resampling over 4 months is more tricky : it is not a real resampling, as some months are repeated. We still need 4 seasonal values ...

Thanks,

Olivier

[mathause]

You can use xr.dot to 'emulate' groupby (xr.dot only works for DataArrays but it should get you there). xr.dot does the same as (air.air * seasons).sum("time") but more efficient.

import xarray as xr

# load example dataset and create monthly means
air = xr.tutorial.open_dataset("air_temperature")
air = air.resample(time="M").mean()

# find months belonging to season
DJFM = air.time.dt.month.isin([12, 1, 2, 3])
MAMJ = air.time.dt.month.isin([3, 4, 5, 6])

concat_dim = xr.DataArray(["DJFM", "MAMJ"], dims="season")

seasons = xr.concat([DJFM, MAMJ], dim=concat_dim)

# calc mean over seasons
xr.dot(air.air, seasons)

[oliviermarti]

This will compute wrong results for DJFM : il will average the JFM values at the beginning of the year with December of the same year. It should use Dec of previous year.

Is there a way to have a correct computation ? it would really help !

[dcherian]

Totally correct!

@oliviermarti you might find some useful code here: https://github.com/pydata/xarray/pull/8509/files#diff-055aed795aa245af3db06b1b2026d7034fcec0d8604effdd501eb1ac54db0e4eR141

[mathause]

Also maybe https://xcdat.readthedocs.io/en/latest/generated/xarray.Dataset.temporal.group_average.html (but I have no experience with xCDAT)

[trexfeathers]

This should work, which I tested on this file:

from ncdata.iris_xarray import cubes_to_xarray, cubes_from_xarray
import iris
from iris.analysis import MEAN
from iris.coord_categorisation import add_season, add_season_year
from iris.cube import CubeList
import xarray as xr


dataset = xr.open_dataset("my_file.nc")
(cube,) = cubes_from_xarray(dataset)


def season_cube(seasons: list[str]):
    cube_new = cube.copy()
    add_season(cube_new, "time", name="season", seasons=seasons)
    add_season_year(cube_new, "time", name=f"season_year", seasons=seasons)
    cube_mean = cube_new.aggregated_by(["season", "season_year"], MEAN)

    seasons_four = [s for s in seasons if len(s) == 4]
    cube_four = cube_mean.extract(iris.Constraint(season=seasons_four))

    return cube_four


all_fragments = CubeList()
for seasons in (
    ["DJFM", "AM", "JJAS", "ON"],
    ["JF", "MAMJ", "JA", "SOND"],
):
    cube_generated = season_cube(seasons)
    all_fragments.extend(cube_generated.slices_over(["time"]))

cube_seasons = all_fragments.merge_cube()

dataset_seasons = cubes_to_xarray(cube_seasons)
print(dataset_seasons)
print(dataset_seasons.season.values)
print(dataset_seasons.season_year.values)

(The ncdata package has been created to allow Xarray and Iris to augment each other by combining their strengths).

[trexfeathers]

But it don’t take into account the correct length of each month :-(

Yes, Iris is not making any assumptions about this, just meaning all the data points provided. The Cube.aggregated_by() function can perform a weighted mean by including the weights= argument. What you input for weights= will depend on how much time each point represents.

[dcherian]

for seasons in (
    ["DJFM", "AM", "JJAS", "ON"],
    ["JF", "MAMJ", "JA", "SOND"],
):

This is quite clever! I've been wondering about how to do the sliding window resampling.

[oliviermarti]

Iris correctly jump over year change :-) But it don’t take into account the correct length of each month :-( dataset.Var.values = array([ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30., 31., 32., 33., 34., 35., 36.]) dataset_seasons.Var.values = array([ 2. , 4.5, 7.5, 10.5, 13.5, 16.5, 19.5, 22.5, 25.5, 28.5, 31.5, 34.5, 36. ]) Olivier

Le 17 avr. 2024 à 13:36, Martin Yeo ***@***.***> a écrit : This should work:

— Olivier Marti LSCE Bât 714 p. 1049 MERMAID Team Office : +33 1 69 08 77 27 Mobile : +33 6 45 36 43 74 ***@***.***