How turn this monthly xarray dataset into an annual mean without resampling?

user_4561892 0
Viewed 1613

I have an xarray of monthly average surface temperatures read in from a server using open_dataset with decode_times=False because the calendar type is not understood by xarray.

After some manipulation, I am left with a dataset my_dataset of surface temperatures ('ts') and times ('T'):

<xarray.Dataset>
Dimensions:  (T: 1800)
Coordinates:
  * T        (T) float32 0.5 1.5 2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5 ...
Data variables:
    ts       (T) float64 246.6 247.9 250.7 260.1 271.9 281.1 283.3 280.5 ...

'T' has the following attributes:

Attributes:
    pointwidth:  1.0
    calendar:    360
    gridtype:    0
    units:       months since 0300-01-01

I would like to take this monthly data and calculate annual averages, but because the T coordinate aren't datetimes, I'm unable to use xarray.Dataset.resample. Right now, I am simply converting to a numpy array, but I would like a way to do this preserving the xarray dataset.

My current, rudimentary way:

temps = np.mean(np.array(my_dataset['ts']).reshape(-1,12),axis=1)
years = np.array(my_dataset['T'])/12

I appreciate any help, even if the best way is redefining the time coordinate to use resampling.

Edit: Requested how xarray was created, it was done via the following:

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

filename = 'http://strega.ldeo.columbia.edu:81/CMIP5/.byScenario/.abrupt4xCO2/.atmos/.mon/.ts/ACCESS1-0/r1i1p1/.ts/dods'
ds = xr.open_dataset(filename,decode_times=False)

zonal_mean = ds.mean(dim='lon')
arctic_only = zonal.where(zonal['lat'] >= 60).dropna('lat')
weights = np.cos(np.deg2rad(arctic['lat']))/np.sum(np.cos(np.deg2rad(arctic['lat'])))
my_dataset = (arctic_only * weights).sum(dim='lat')
python-xarray python-datetime numpy python
1 Answers

This is a very common problem, especially with datasets from INGRID. The reason xarray can't decode the date whose units are "months since..." is due to the underlying netcdf4-python library's refusal to parse such dates. This is discussed in a netcdf4-python github issue

The problem with time units such as "months" is that they are not well defined. In contrast to days, hours, etc. the length of a month depends on the calendar used and even varies between different months.

INGRID unfortunately refuses to accept this fact and continues to use "months" as its default unit, despite the ambiguity. So right now there is this frustrating incompatibility between INGRID and xarray / python-netcdf4.

Anyway, here is a hack to accomplish what you want without leaving xarray

# create new coordinates for month and year
ds.coords['month'] = np.ceil(ds['T'] % 12).astype('int') 
ds.coords['year'] = (ds['T'] // 12).astype('int')
# calculate monthly climatology
ds_clim = ds.groupby('month').mean(dim='T')
# calculate annual mean
ds_am = ds.groupby('year').mean(dim='T')
user_32662350
Related