[QUESTION] compare_zonal_mean not reading in Xarray.Dataset files?

[jmoch1214]

For fields where I've used gcpy.open_mfdataset to open a list of files and then .mean(dim='time',keep_attrs=True) to get an annual average, gcpy.compare_zonal_mean is rejecting the fields saying they are not an axarray DataArray, even though the fields are definitely Xarray Datasets and gcpy.compare_single_level successfully loads the same exact files and makes figures with no problem.

The specific error I am getting is:

---

TypeError Traceback (most recent call last)
in
1 # Zonal mean - PDF
2 pdfname = os.path.join(plotsdir,'ZonalMean_{}_{}diag.pdf'.format(desc, comparison))
----> 3 gcpy.compare_zonal_mean( refdata, refstr, devdata, devstr, varlist=varlist, pdfname=pdfname )
4 #gcpy.add_bookmarks_to_pdf( pdfname, varlist, remove_prefix='SpeciesConc')

~/Py_notebooks/gcpy/gcpy/benchmark.py in compare_zonal_mean(refdata, refstr, devdata, devstr, varlist, itime, weightsdir, pdfname, cmpres, match_cbar, pres_range, normalize_by_area, enforce_units, flip_ref, flip_dev, use_cmap_RdBu)
1050 for i in range(6):
1051 regridder = refregridder_list[i]
-> 1052 ds_ref_cmp += regridder(ds_ref_reshaped[i])
1053 else:
1054 ds_ref_cmp = ds_ref

/n/jacob_lab/Lab/python/miniconda/envs/gcpy/lib/python3.6/site-packages/xesmf/frontend.py in call(self, a)
260 regrid_func = self.regrid_dataarray
261 else:
--> 262 raise TypeError("input must be numpy array or xarray DataArray!")
263
264 return regrid_func(a)

TypeError: input must be numpy array or xarray DataArray!

[lizziel]

gcpy.open_mfdataset is a wrapper around xarray function open_mfdataset which reads in netcdf variables as dask arrays. See http://xarray.pydata.org/en/stable/dask.html for more information on using dask arrays in xarray.

For example, I used open_mfdataset to read in several files and my resulting data looks like this:

<xarray.Dataset>
Dimensions: (lat: 288, lev: 72, lon: 48)
Coordinates:

• lon (lon) float64 1.0 2.0 3.0 4.0 5.0 ... 45.0 46.0 47.0 48.0
• lat (lat) float64 1.0 2.0 3.0 4.0 5.0 ... 285.0 286.0 287.0 288.0
• lev (lev) float64 1.0 2.0 3.0 4.0 5.0 ... 69.0 70.0 71.0 72.0
  Data variables:
  Jval_O3O3P (lev, lat, lon) float32 dask.array<shape=(72, 288, 48), chunksize=(72, 288, 48)>
  Jval_O3O1D (lev, lat, lon) float32 dask.array<shape=(72, 288, 48), chunksize=(72, 288, 48)>
  Jval_O2 (lev, lat, lon) float32 dask.array<shape=(72, 288, 48), chunksize=(72, 288, 48)>
  Jval_INPN (lev, lat, lon) float32 dask.array<shape=(72, 288, 48), chunksize=(72, 288, 48)>
  Jval_MRP (lev, lat, lon) float32 dask.array<shape=(72, 288, 48), chunksize=(72, 288, 48)>

This will cause problems when regridding since xesmf expects numpy array or xarray DataArray. When you ran your code for compare_single_level you likely did not regrid and so did not run into this problem. However, when you pass the same cubed-sphere data to compare_zonal_mean it will automatically be regridded to calculate zonal mean.

To avoid this problem you can convert your dask arrays to numpy arrays or xarray DataArrays prior to passing the data to the plotting functions. For example, np.array(ds_ref['varname'].data). You can also create an issue for xesmf to expand the handling to include dask arrays.