[Bug]: Metadata regridding issues related to multidimensional coords

[pochedls]

What happened?

xcdat isn't recognizing dimensions when regridding with xesmf, but regridding does work when using xesmf directly.

Minimal Complete Verifiable Example (MVCE)

xcdat_example.tar.gz

# %% imports
import xcdat as xc
import xesmf as xe

# %% parameters
fn_source = 'source.nc'
fn_target = 'target.nc'

# %% get source / target grids
ds = xc.open_dataset('source.nc')
ngrid = xc.open_dataset('target.nc')

# %% try regridding with xcdat
dsr = ds.regridder.horizontal('tmp2m', ngrid, tool="xesmf", method="bilinear")

KeyError: "No 'X' axis dimension coordinate variables were found in the xarray object. Make sure dimension coordinate variables exist, they are one dimensional, and their CF 'axis' or 'standard_name' attrs are correctly set."

# %% inspect dataset
print(ds.tmp2m)

<xarray.DataArray 'tmp2m' (y: 1059, x: 1799)> Size: 8MB
array([[296.45013, 296.45013, 296.38763, ..., 299.76263, 299.70013, 299.70013],
...,
[282.76263, 282.76263, 282.76263, ..., 270.07513, 270.01263, 270.01263]],
shape=(1059, 1799), dtype=float32)
Coordinates:
latitude (y, x) float64 15MB ...
longitude (y, x) float64 15MB ...

print(ds.longitude)

<xarray.DataArray 'longitude' (y: 1059, x: 1799)> Size: 15MB
[1905141 values with dtype=float64]
Coordinates:
heightAboveGround float64 8B ...
latitude (y, x) float64 15MB ...
longitude (y, x) float64 15MB ...
time object 8B ...
Dimensions without coordinates: y, x
Attributes:
units: degrees_east
standard_name: longitude
long_name: longitude
axis: X

# %% try regridding with xesmf
regridder = xe.Regridder(ds, ngrid, "bilinear")
dsr = regridder(ds)

Log Output

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
Cell In[4], [line 2](vscode-notebook-cell:?execution_count=4&line=2)
      1 # %% try regridding with xcdat
----> [2](vscode-notebook-cell:?execution_count=4&line=2) dsr = ds.regridder.horizontal("tmp2m", ngrid, tool="xesmf", method="bilinear")

File /global/u2/v/vo13/xCDAT/xcdat/xcdat/regridder/accessor.py:169, in RegridderAccessor.horizontal(self, data_var, output_grid, tool, **options)
    164 except KeyError as e:
    165     raise ValueError(
    166         f"Tool {e!s} does not exist, valid choices {list(HORIZONTAL_REGRID_TOOLS)}"
    167     ) from e
--> [169](https://vscode-remote+ssh-002dremote-002bperlmutter.vscode-resource.vscode-cdn.net/global/u2/v/vo13/xCDAT/xcdat/xcdat/regridder/accessor.py:169) input_grid = _get_input_grid(self._ds, data_var, ["X", "Y"])
    170 regridder = regrid_tool(input_grid, output_grid, **options)
    171 output_ds = regridder.horizontal(data_var, self._ds)

File /global/u2/v/vo13/xCDAT/xcdat/xcdat/regridder/accessor.py:377, in _get_input_grid(ds, data_var, dup_check_dims)
    374 all_coords = set(ds.coords.keys())
    376 for dimension in dup_check_dims:
--> [377](https://vscode-remote+ssh-002dremote-002bperlmutter.vscode-resource.vscode-cdn.net/global/u2/v/vo13/xCDAT/xcdat/xcdat/regridder/accessor.py:377)     coords = get_dim_coords(ds, dimension)
    379     if isinstance(coords, xr.Dataset):
    380         coord = set([get_dim_coords(ds[data_var], dimension).name])

File /global/u2/v/vo13/xCDAT/xcdat/xcdat/axis.py:141, in get_dim_coords(obj, axis)
    134     raise ValueError(
    135         f"This DataArray has more than one dimension {dim_coord_keys} mapped to the "
    136         f"'{axis}' axis, which is an unexpected behavior. Try dropping extraneous "
    137         "dimensions from the DataArray first (might affect data shape)."
    138     )
    140 if len(dim_coord_keys) == 0:
--> [141](https://vscode-remote+ssh-002dremote-002bperlmutter.vscode-resource.vscode-cdn.net/global/u2/v/vo13/xCDAT/xcdat/xcdat/axis.py:141)     raise KeyError(
    142         f"No '{axis}' axis dimension coordinate variables were found in the "
    143         f"xarray object. Make sure dimension coordinate variables exist, they are "
    144         "one dimensional, and their CF 'axis' or 'standard_name' attrs are "
    145         "correctly set."
    146     )
    148 dim_coords = obj[
    149     dim_coord_keys if len(dim_coord_keys) > 1 else dim_coord_keys[0]
    150 ].copy()
    152 return dim_coords

KeyError: "No 'X' axis dimension coordinate variables were found in the xarray object. Make sure dimension coordinate variables exist, they are one dimensional, and their CF 'axis' or 'standard_name' attrs are correctly set."

Anything else we need to know?

This may be related to #718 and #747.

Environment

xCDAT : 0.10.1
xesmf: 0.8.10

[tomvothecoder]

Jason is starting to look at this issue. I am also investigating as it might be good to include with v0.11.0.

[tomvothecoder]

@pochedls @jasonb5

Root Cause

The issue stems from the dataset having multidimensional coordinates. Notice how latitude and longitude have (y, x) dimensions:

<xarray.Dataset> Size: 23MB
Dimensions:            (y: 1059, x: 1799)
Coordinates:
    heightAboveGround  float64 8B ...
    latitude           (y, x) float32 8MB ...
    longitude          (y, x) float32 8MB ...
    time               object 8B ...
Dimensions without coordinates: y, x
Data variables:
    tmp2m              (y, x) float32 8MB ...
    forecast_times     int64 8B ...

In _get_input_grid(), we use xc.get_dim_coords() to map to axes. Xarray does not recognize multidimensional coordinates as index keys, as noted in the code comment below.

xcdat/xcdat/axis.py

Lines 120 to 123 in e9937b6

	# Get the object's index keys, with each being a dimension.
	# NOTE: xarray does not include multidimensional coordinates as index keys.
	# Example: ["lat", "lon", "time"]
	index_keys = obj.indexes.keys()

As a result, xCDAT cannot retrieve the axes, leading to the error: KeyError: "No 'X' axis dimension coordinate variables were found in the xarray object. Make sure dimension coordinate variables exist, they are one dimensional, and their CF 'axis' or 'standard_name' attrs are correctly set."

Workaround/Solution

I'm pretty sure we've ran into this in the past before.

We need to extend support for mapping to multidimensional coords (possibly use cf_xarray directly). I'm looking at two approaches:

1. Extend get_dim_coords() -- more invasive
2. Update get_input_grid() -- less invasive, but less maintainable overall

Example:

ds.cf

Coordinates:
             CF Axes:   X: ['longitude']
                        Y: ['latitude']
                        T: ['time']
                        Z: n/a

      CF Coordinates:   longitude: ['longitude']
                        latitude: ['latitude']
                        vertical: ['heightAboveGround']
                        time: ['time']

       Cell Measures:   area, volume: n/a

      Standard Names:   height: ['heightAboveGround']
                        latitude: ['latitude']
                        longitude: ['longitude']
                        time: ['time']

              Bounds:   n/a

       Grid Mappings:   n/a

Data Variables:
       Cell Measures:   area, volume: n/a

      Standard Names:   air_temperature: ['tmp2m']

              Bounds:   n/a

       Grid Mappings:   n/a

[tomvothecoder]

Another possible issue: xCDAT doesn't support adding bounds for multidimensional coords too.

ds = xc.open_dataset("source.nc", add_bounds=("X","Y")) # Default

print(ds)

<xarray.Dataset> Size: 23MB
Dimensions:            (y: 1059, x: 1799)
Coordinates:
    heightAboveGround  float64 8B ...
    latitude           (y, x) float32 8MB ...
    longitude          (y, x) float32 8MB ...
    time               object 8B ...
Dimensions without coordinates: y, x
Data variables:
    tmp2m              (y, x) float32 8MB ...
    forecast_times     int64 8B ...

[pochedls]

Another possible issue: xCDAT doesn't support adding bounds for multidimensional coords too.

I think this is an issue for methods like conservative_normed, but should not affect bilinear regridding. Note that I could get regridding to work using xesmf directly.

[jasonb5]

This issue does sound familiar.

I'm thinking we could add a flag (default false) to get_dim_coords to allow returning multi-dimensional coordinates and have the xesmf regridder utilize that.

I believe regrid2 cannot handle multi-dimensional coordinates, not sure if it does so gracefully.

If there's a use case for xCDAT to generate bounds for these grids we could work on adding that.