Calculating anomalies from a specific climatology

[meganf11]

I am trying to calculate seasonal anomalies from monthly gridded data spanning from 1870-2100. I was excited to use xCDAT because of how easy it is to specify CORRECT seasons, such as DJF, but I am interested in taking seasonal anomalies from a 30-year baseline climatology (1971-2000). Any help is much appreciated, thank you!

[pochedls]

@meganf11 – thanks for checking out xCDAT. You're right that xCDAT currently takes departures relative to the full dataset period. There is a Feature Request open to allow the user to select the baseline period (or provide an xarray dataarray/dataset with the climatology). I outlined a possible workaround in that feature request, but I didn't try it (I imagine it doesn't "just work" and you might need to use pandas to do some sort of group subtraction).

@chengzhuzhang / @lee1043 – is there an easy workaround to take the seasonal departures with respect to a user-specified base period?

@meganf11 – if you do find a work around, can you update this discussion? Or if you're willing to try a PR for that feature request, there is a pretty good xCDAT contributing guide here (and @tomvothecoder could likely weigh in with pointers since he put together the original departures code).

[meganf11]

Thanks for all the input! I am trying the code that you suggested a few comments down, but its telling me that no bounds were found for the "T" axis and I am not sure how to interpret this. I tried going to the API for the function but had little luck.

[pochedls]

Did you open the dataset with xcdat? Usually xcdat will generate time bounds when the datset is opened (if they don't exist). Maybe you could excerpt your code and the error message in this thread?

[lee1043]

Hi @meganf11, We are excited to hear your interest in using xCDAT! Getting seasonal average correct is one of many things the xCDAT core team has spent many time on it! As @pochedls discussed above (thanks @pochedls for your comment and for pinging me!), we are yet to have the capability of getting anomaly of full time series from a specific baseline time window that is different than that of the full time series, but we definitely have a plan to add that capability.

During the meantime, I have crafted a cookbook for a workaround, which needs some extra steps as follows.

1. Get seasonal average from the original full time series
2. Get climatological seasonal average from subset of the full time series for the baseline period
3. Match time dimension: Repeat climatological seasonal average of 2 along the time axis of 1, to enable subtraction
4. Get anomaly by subtracting 3 from 1

Please feel free to make any suggestions if the cookbook is not clear enough. It is in the beginning stage and any feedback would be appreciated. Let me know if you have any question.

[pochedls]

Thanks @lee1043 – I wanted to see if we could avoid step 3 above using xr.groupby. I think I figured out a way to do this:

# imports
import xcdat as xc
# I/O
fn = '/path/to/data/data.nc'
ds = xc.open_dataset(fn)
# seasonal average
ds_season = ds.temporal.group_average('tas', freq="season")
# get seasonal average during preferred baseline
ds_reference = ds_season.sel(time=slice("1971-01-01", "2000-12-31"))
# subtract time mean of reference seasonal climatology from full seasonal time series
# note indexing on month – seasonal data is saved on Jan. 1, Apr. 1, Jul. 1, and Oct. 1 – works here, but always reliable?
tas_anomaly = ds_season.tas.groupby("time.month") - ds_reference.tas.groupby("time.month").mean("time")
# put time series in seasonal dataset
ds_season['tas_anomaly'] = tas_anomaly
# example spatial average computation
ts = ds_season.spatial.average("tas_anomaly").tas_anomaly

Your comment also makes me wonder if we should have a grower utility.

[lee1043]

@pochedls thank you for the inspiration for simplifying the code. Following on that, I updated my cookbook, which I think now it should serve well the purpose of @meganf11, using the advantage of xcdat that takes correct (i.e., continuous DJF) seasons. Using groupby only with time.month could be useful for annual cycle departure, but in season case, we would prefer using temporal.group_average with freq="season".

Below is part of the updated cookbook, which I merged above 3 & 4 steps by using combination of xcdat's temporal.group_average and groupby of xarray. I verified the anomaly result is consistent to my previous cookbook that has grower concept.

# seasonal average time series
ds_season = ds.temporal.group_average(data_var, freq="season")   

# time subsetted data
ds_reference = ds.sel(time=slice("1971-01-01", "2000-12-31"))

# climatology of subsetted data
ds_clim = ds_reference.temporal.climatology(data_var, freq="season")  

# anomaly calculation
ds_season[data_var+'_departure'] = ds_season[data_var].groupby("time.season") - ds_clim[data_var].groupby("time.season").mean("time") 
# note: "time.month" also can be used with no difference in results, but "time.season" is more conceptually precise in this case

[chengzhuzhang]

Thanks for putting together the cookbook @lee1043! And thanks for the codes for simplifying, @pochedls.

I was looking into the new feature requested by @pochedls here:. The departure function @tomvothecoder wrote (

xcdat/xcdat/temporal.py

Line 474 in 1a0afb4

def departures(

) also use the groupBy before subtracting the reference climatology. It might be worthwhile to add time slice selection function sel in temporal.climatology through a reference_period or more general sel_time_slice argument, so that both temporal.climatology and temporal.departures can have a way to select time (reference time) slice.

Regarding to grower, not sure if there is an equivalent utility in xarray

[tomvothecoder]

Yup, I use xarray's grouped arithmetic to calculate anomalies.

From Xarray example on calculating anomalies:

climatology = ds.groupby("time.month").mean("time")
anomalies = ds.groupby("time.month") - climatology

[mfwehner]

This is a not uncommon problem In the cdat days the fix was to call cdutil.setTimeBounds with an appropriate temporal argument No idea what to do in xcdat Michael Dictated to my iPhone, please excuse Siri's grammatical errorsOn Jan 24, 2023, at 10:29 AM, meganf11 ***@***.***> wrote: Thanks for all the input! I am trying the code that you suggested a few comments down, but its telling me that no bounds were found for the "T" axis and I am not sure how to interpret this. I tried going to the API for the function but had little luck. —Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>