[Bug]: ValueError: Cannot generate bounds for multidimensional coordinates.

[pochedls]

What happened?

I'm running into an issue working with a number of relatively basic (seemingly okay) netCDF files. I get the following error: ValueError: Cannot generate bounds for multidimensional coordinates.

What did you expect to happen?

I don't think that an error should be generated in these cases.

Minimal Complete Verifiable Example

import xcdat

fn = '/p/css03/esgf_publish/CMIP6/CMIP/CNRM-CERFACS/CNRM-CM6-1/historical/r3i1p1f2/Amon/tas/gr/v20190125/tas_Amon_CNRM-CM6-1_historical_r3i1p1f2_gr_185001-201412.nc'

ds = xcdat.open_dataset(fn)

Relevant log output

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/code/xcdat/xcdat/bounds.py:178, in BoundsAccessor.get_bounds(self, axis)
    177 try:
--> 178     bounds_key = coord_var.attrs["bounds"]
    179 except KeyError:

KeyError: 'bounds'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/code/xcdat/xcdat/bounds.py:145, in BoundsAccessor.add_missing_bounds(self, width)
    144 try:
--> 145     self.get_bounds(axis)
    146 except KeyError:

File ~/code/xcdat/xcdat/bounds.py:180, in BoundsAccessor.get_bounds(self, axis)
    179 except KeyError:
--> 180     raise KeyError(
    181         f"The coordinate variable '{coord_var.name}' has no 'bounds' attr. "
    182         "Set the 'bounds' attr to the name of the bounds data variable."
    183     )
    185 try:

KeyError: "The coordinate variable 'height' has no 'bounds' attr. Set the 'bounds' attr to the name of the bounds data variable."

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/code/xcdat/xcdat/bounds.py:178, in BoundsAccessor.get_bounds(self, axis)
    177 try:
--> 178     bounds_key = coord_var.attrs["bounds"]
    179 except KeyError:

KeyError: 'bounds'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
File ~/code/xcdat/xcdat/bounds.py:220, in BoundsAccessor.add_bounds(self, axis, width)
    219 try:
--> 220     self.get_bounds(axis)
    221     raise ValueError(
    222         f"{axis} bounds already exist. Drop them first to add new bounds."
    223     )

File ~/code/xcdat/xcdat/bounds.py:180, in BoundsAccessor.get_bounds(self, axis)
    179 except KeyError:
--> 180     raise KeyError(
    181         f"The coordinate variable '{coord_var.name}' has no 'bounds' attr. "
    182         "Set the 'bounds' attr to the name of the bounds data variable."
    183     )
    185 try:

KeyError: "The coordinate variable 'height' has no 'bounds' attr. Set the 'bounds' attr to the name of the bounds data variable."

During handling of the above exception, another exception occurred:

ValueError                                Traceback (most recent call last)
Input In [1], in <cell line: 5>()
      1 import xcdat
      3 fn = '/p/css03/esgf_publish/CMIP6/CMIP/CNRM-CERFACS/CNRM-CM6-1/historical/r3i1p1f2/Amon/tas/gr/v20190125/tas_Amon_CNRM-CM6-1_historical_r3i1p1f2_gr_185001-201412.nc'
----> 5 ds = xcdat.open_dataset(fn)

File ~/code/xcdat/xcdat/dataset.py:95, in open_dataset(path, data_var, add_bounds, decode_times, center_times, lon_orient, **kwargs)
     92 else:
     93     ds = xr.open_dataset(path, decode_times=False, **kwargs)
---> 95 ds = _postprocess_dataset(ds, data_var, center_times, add_bounds, lon_orient)
     97 return ds

File ~/code/xcdat/xcdat/dataset.py:482, in _postprocess_dataset(dataset, data_var, center_times, add_bounds, lon_orient)
    479         raise ValueError("This dataset does not have a time coordinates to center.")
    481 if add_bounds:
--> 482     dataset = dataset.bounds.add_missing_bounds()
    484 if lon_orient is not None:
    485     if dataset.cf.dims.get("X") is not None:

File ~/code/xcdat/xcdat/bounds.py:147, in BoundsAccessor.add_missing_bounds(self, width)
    145             self.get_bounds(axis)
    146         except KeyError:
--> 147             self._dataset = self.add_bounds(axis, width)
    148 return self._dataset

File ~/code/xcdat/xcdat/bounds.py:225, in BoundsAccessor.add_bounds(self, axis, width)
    221     raise ValueError(
    222         f"{axis} bounds already exist. Drop them first to add new bounds."
    223     )
    224 except KeyError:
--> 225     dataset = self._add_bounds(axis, width)
    227 return dataset

File ~/code/xcdat/xcdat/bounds.py:269, in BoundsAccessor._add_bounds(self, axis, width)
    267 # Validate coordinate shape and dimensions
    268 if coord_var.ndim != 1:
--> 269     raise ValueError("Cannot generate bounds for multidimensional coordinates.")
    270 if coord_var.shape[0] <= 1:
    271     raise ValueError("Cannot generate bounds for a coordinate of length <= 1.")

ValueError: Cannot generate bounds for multidimensional coordinates.

Anything else we need to know?

No response

Environment

/home/pochedley1/bin/anaconda3/envs/xcdat/lib/python3.10/site-packages/_distutils_hack/init.py:33: UserWarning: Setuptools is replacing distutils.
warnings.warn("Setuptools is replacing distutils.")

INSTALLED VERSIONS

commit: None
python: 3.10.5 | packaged by conda-forge | (main, Jun 14 2022, 07:04:59) [GCC 10.3.0]
python-bits: 64
OS: Linux
OS-release: 3.10.0-1160.62.1.el7.x86_64
machine: x86_64
processor: x86_64
byteorder: little
LC_ALL: None
LANG: en_US.UTF-8
LOCALE: ('en_US', 'UTF-8')
libhdf5: 1.12.1
libnetcdf: 4.8.1

xarray: 2022.6.0
pandas: 1.4.3
numpy: 1.22.4
scipy: 1.8.1
netCDF4: 1.6.0
pydap: None
h5netcdf: None
h5py: None
Nio: None
zarr: None
cftime: 1.6.1
nc_time_axis: None
PseudoNetCDF: None
rasterio: None
cfgrib: None
iris: None
bottleneck: None
dask: 2022.7.1
distributed: 2022.7.1
matplotlib: 3.5.2
cartopy: None
seaborn: None
numbagg: None
fsspec: 2022.5.0
cupy: None
pint: None
sparse: 0.13.0
flox: None
numpy_groupies: None
setuptools: 63.2.0
pip: 22.2
conda: None
pytest: None
IPython: 8.4.0
sphinx: None

[pochedls]

The issue seems to be that add_missing_bounds loops over axes, which are statically defined in axis.py (here) as X, Y, T, and Z.

When add_missing_bounds gets to Z it finds a singleton value (2 for 2m above ground) with no bounds and throws the error above (here).

Perhaps we could simply not return bounds in these cases.

[pochedls]

Or perhaps the coordinate shape validation could be skipped if the axis in question is not a dimension (i.e., is not in ds.dims).

[tomvothecoder]

Hi @pochedls, thanks for reporting this issue. It sounds like you're on the right track with possible solutions.
Would you be willing to attempt a PR and I'll review?