Including dataarrays with our current dataset API model

[pochedls]

xcdat has generally operated on dataset objects, which would be comparable to opening all of the contents of a netCDF4 file. This would be similar to the following objects:

• ds = Dataset(filename) for the netCDF4 library
• ds = cdms2.open(filename) for cdms
• ds = xr.open_dataset(filename) for xarray / xcdat

We do this because the dataset (ds) object has everything you need for most calculations (i.e., a dataset contains the bounds, which can be used to calculate spatial weights for spatial averages). So, for a variable tas, we can do things like:

• calculate a spatial average (ds = ds.spatial.average("tas", ...))
• regrid (ds = ds.regridder.horizontal("tas", ...))
• calculate departures (ds = ds.temporal.departures("tas", ...))

But we'd ideally just operate directly on the variable/dataarray tas. Imagine if we could do the following:

ds = xc.open_dataset(filename)
tas = ds['tas']
tas_global_mean = tas.spatial.average()
tas_global_mean_anomaly = tas_global_mean.temporal.departures(freq="month")

(Or you could imagine functions like tas_global = spatial_average(tas))...

This is how a lot of people are used to operating. CDAT actually operated this way. The problem has been getting bounds into the dataarray object. This has been a longstanding issue/request for xarray.

Could we help address that longstanding issue (or would some of the simpler approaches work for us)?

[pochedls]

One approach outlined in xarray's issue could potentially be implemented (note you need to download a file listed at the bottom to run this example

Setup

import xcdat as xc
import numpy as np
import xarray as xr

def getBoundedDataArray(ds, v):
    # loop over coordinates to get coordinates and coordinate bounds
    coords = {}
    for c_key in ds[v].cf.axes.keys():
        try:
            # get dimension key (e.g., "time", "lat", "lon")
            dim_key = xc.get_dim_keys(ds, c_key)
            # add axis to coordinate dict
            coords[dim_key] = ds[dim_key]
            # get coordinate dtype
            dim_value_dtype = ds[dim_key].dtype
            # create a bounds dtype (based on coordinate dtype)
            bounds_dtype = np.dtype([('start', dim_value_dtype), ('stop', dim_value_dtype)])
            # get the bounds for axis
            bnds = ds.bounds.get_bounds(axis=c_key)
            # convert to expected form for bounds_dtype
            newbnds = [tuple(row) for row in bnds.values]
            # create new bounds object
            newbnds = np.array(newbnds, dtype=bounds_dtype)
            # add bounds to coordinate dict
            coords[bnds.name] = (dim_key, newbnds)
        except:
            continue
    # return dataarray with bounds
    return xr.DataArray(ds[v], coords=coords, dims=ds[v].dims, attrs=ds[v].attrs)

Result

ds = xc.open_dataset('tas_Amon_MIROC-ES2L_historical_r2i1p1f2_gn_185001-201412.nc')
da = getBoundedDataArray(ds, 'tas')
da

<xarray.DataArray 'tas' (time: 1980, lat: 64, lon: 128)>
array([[[263.05438, 262.90805, 262.76172, ..., 263.49677, 263.3493 ,
263.20184],
...
[241.66898, 241.59831, 241.52194, ..., 241.88553, 241.80688,
241.73735]]], dtype=float32)
Coordinates:

• lon (lon) float64 0.0 2.812 5.625 8.438 ... 348.8 351.6 354.4 357.2
  lon_bnds (lon) [('start', '<f8'), ('stop', '<f8')] (-1.40625, 1.40625) ...
• lat (lat) float64 -87.86 -85.1 -82.31 -79.53 ... 82.31 85.1 87.86
  lat_bnds (lat) [('start', '<f8'), ('stop', '<f8')] (-90., -86.57774751)...
• time (time) object 1850-01-16 12:00:00 ... 2014-12-16 12:00:00
  time_bnds (time) [('start', 'O'), ('stop', 'O')] (cftime.DatetimeGregori...
  Attributes:
  standard_name: air_temperature
  long_name: Near-Surface Air Temperature
  ...

Selection works...

da.isel(time=0).time_bnds.values

array((cftime.DatetimeGregorian(1850, 1, 1, 0, 0, 0, 0, has_year_zero=False), cftime.DatetimeGregorian(1850, 2, 1, 0, 0, 0, 0, has_year_zero=False)), dtype=[('start', 'O'), ('stop', 'O')])

File download

https://aims3.llnl.gov/thredds/fileServer/css03_data/CMIP6/CMIP/MIROC/MIROC-ES2L/historical/r2i1p1f2/Amon/tas/gn/v20190823/tas_Amon_MIROC-ES2L_historical_r2i1p1f2_gn_185001-201412.nc.

[pochedls]

There's other potential solutions discussed at xarray, but I think this basic example might work for our purposes and could help us understand the issue better (to potentially contribute to xarray's work on this issue).

[pochedls]

One other note on this example. You might be able to register the __call__ function to enable CDAT-like behavior (the arguments in getBoundedDataArray would need to change to (self, v)):

xr.Dataset.__call__ = getBoundedDataArray

ds = xc.open_dataset('tas_Amon_MIROC-ES2L_historical_r2i1p1f2_gn_185001-201412.nc')
tas = ds('tas')