元cross_section函数不识别ERA5中的坐标

Question

我使用ERA5数据集来绘制南部非洲垂直速度的横截面。数据的一个示例是这里。我尝试复制给定的这里示例，但是，由于数据集的坐标中缺少信息，cross_section函数无法工作。我使用的代码如下：

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)
#

这些数据如下：

<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...

继续提供我想要横截面的坐标：

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

最终，cross_section函数会产生以下错误消息，与ERA5 5的坐标缺失信息有关：

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.

对周围的工作有什么帮助吗？

谢谢!

Answer 1

从您的数据集输出中可以看到，有一个变量(time_bnds)没有水平坐标变量，对于传递给metpy.interpolate.cross_section (从文档中可以看到)的数据集中的每个数据变量，都需要这些变量。因此，在使用cross_section之前，我们必须排除这些数据变量。

对于time_bnds这样的“类坐标”变量，一个简单的修复方法就是将它提升到一个辅助坐标，如下所示：

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'))

这也是与https://github.com/Unidata/MetPy/discussions/2033相同的问题，因此，请阅读该讨论，以获得更多关于如何在未来的MetPy中改进这一点的信息和想法！