metpy cross_section function does not recognize coordinates in ERA5

Question

I use the ERA5 dataset in order to plot a cross section of vertical velocity over southern Africa. A sample of the data is found here. I try to replicate the example given here, however, the cross_section function does not work, as there is missing information in the dataset's coordinates. The code I use is the following:

from matplotlib import pyplot
from matplotlib.cm import get_cmap
from __future__ import print_function
from netCDF4 import Dataset,num2date,date2num
from matplotlib.colors import from_levels_and_colors
from cartopy import crs
from cartopy.feature import NaturalEarthFeature, COLORS
from metpy.units import units
from datetime import datetime
from metpy.plots import StationPlot
from metpy.interpolate import cross_section
from PIL import Image, ImageDraw
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import ImageGrid
#
import os
import netCDF4 as nc
import scipy.spatial as ss
import matplotlib.pyplot as plt
import metpy.calc as mpcalc
import xarray as xr
import cartopy.crs as ccrs
import matplotlib
import numpy as np
import datetime
import cartopy.feature as cfeature
import matplotlib.gridspec as gridspec
import scipy.ndimage as ndimage
#
root_dir = '/users/.../my_var/'
#
data = xr.open_dataset('omega.nc')
data = data.metpy.parse_cf().squeeze()
print(data)
#

The data are the following:

<xarray.Dataset>
Dimensions:    (bnds: 2, latitude: 161, level: 16, longitude: 137)
Coordinates:
    time       datetime64[ns] 2005-02-01
  * longitude  (longitude) float32 8.0 8.25 8.5 8.75 ... 41.25 41.5 41.75 42.0
  * latitude   (latitude) float32 -5.0 -5.25 -5.5 -5.75 ... -44.5 -44.75 -45.0
  * level      (level) float64 100.0 150.0 200.0 250.0 ... 750.0 800.0 850.0
    metpy_crs  object Projection: latitude_longitude
Dimensions without coordinates: bnds
Data variables:
    time_bnds  (bnds) datetime64[ns] 2005-01-01 2005-12-01
    w          (level, latitude, longitude) float32 -0.001213 ... 0.003879
Attributes:
    CDI:          Climate Data Interface version 1.9.7.1 (http://mpimet.mpg.d...
    Conventions:  CF-1.6
    history:      Fri Feb 04 23:12:51 2022: cdo yearmean omega_SAF_DJF_1986_2...
    frequency:    year
    CDO:          Climate Data Operators version 1.9.7.1 (http://mpimet.mpg.d...

Proceeding on with providing the coordinates over which I want the cross section:

start = (-7, 17)
end = (-28, 17)
cross = cross_section(data, start, end).set_coords(('latitude', 'longitude'))

Eventually the cross_section function yields the following error message, related to missing information about ERA5's coordinates:

Traceback (most recent call last):
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/metpy/interpolate/slices.py", line 167, in cross_section
    x = data.metpy.x
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/metpy/xarray.py", line 431, in x
    return self._axis('x')
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/metpy/xarray.py", line 385, in _axis
    raise AttributeError(axis + ' attribute is not available.')
AttributeError: x attribute is not available.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/metpy/interpolate/slices.py", line 157, in cross_section
    interp_type=interp_type)
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/xarray/core/dataset.py", line 4592, in map
    for k, v in self.data_vars.items()
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/xarray/core/dataset.py", line 4592, in <dictcomp>
    for k, v in self.data_vars.items()
  File "/users/pr007/mkaryp/.local/lib/python3.6/site-packages/metpy/interpolate/slices.py", line 169, in cross_section
    raise ValueError('Data missing required coordinate information. Verify that '
ValueError: Data missing required coordinate information. Verify that your data have been parsed by MetPy with proper x and y dimension coordinates and added crs coordinate of the correct projection for each variable.

Any help on what a work around would be?

Thank you!

Answer 1

As can be seen in your dataset output, there is a variable (time_bnds) which does not have horizontal coordinate variables, which are required for every data variable in a Dataset passed to metpy.interpolate.cross_section (as can be seen in the documentation). So, we have to exclude such data variables before using cross_section somehow.

A simple fix for such a "coordinate-like" variable as time_bnds would be to just promote it to an auxiliary coordinate, like so:

data = xr.open_dataset('omega.nc')
data = data.metpy.parse_cf().squeeze().set_coords('time_bnds')

#####

start = (-7, 17)
end = (-28, 17)
cross = cross_section(data, start, end).set_coords(('latitude', 'longitude'))

This is also the same issue as https://github.com/Unidata/MetPy/discussions/2033, so feel free to read that discussion for more information and ideas on how this could be improved in MetPy in the future!