from intake import open_catalog

cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/atmosphere.yaml")

ds = cat['era5_hourly_reanalysis_single_levels_sa'].to_dask()
ds

# with flox, this mean calculation should be straightforward
import flox.xarray  # required

# ideally we would use the default method="cohorts" but that's not very optimal at the moment
mean = ds.groupby("time.dayofyear").mean(method="map-reduce")

# A user would do ds.groupby('time.dayofyear') - mean
# but this is what Xarray does under the hood, and it is clearer
# about what's happening
# The chunking to 1 *should* make this work better below, but wouldn't expect the average user to do it.
anomaly = ds - mean.sel(dayofyear=ds.time.dt.dayofyear)
anomaly