使用 Xarray 多文件数据集创建一堆图形

Question

我有一堆 .grib 文件，我希望能够读取它们并为其创建图形。我有加载一个文件并绘制一个图形的代码，但我想找出循环它的最佳方法，这样我就不必更改每个文件的文件名和其他参数。到目前为止，我已从使用 xr.open_dataset 切换到 xr.open_mfdataset。我尝试使用 time 进行 For 循环，但在plot = ax.contourf 处出现错误，提示“TypeError：输入 Z 必须是 2D，而不是 3D。”数据集包含纬度、经度、t2m（2 米气温）和时间。请参阅下面的代码和附图。

ds=xr.open_mfdataset('./Data/ecmwf_dataset/t2m/*.grib', combine='nested', concat_dim='time',
                     engine='cfgrib',backend_kwargs={'indexpath': ''})

lons = ds['longitude'][:]
lats = ds['latitude'][:]
sfcT = ds['t2m']
time = ds['time'].dt.strftime('%m-%d-%Y')

for i in range(len(ds.time)):
    fig = plt.figure(figsize=(10, 10))
    ax = fig.add_subplot(111, projection=ccrs.NorthPolarStereo(central_longitude=-47))     

    ax.set_extent([-180, 180, 60, 90], crs=ccrs.PlateCarree())
    lon2d,lat2d = np.meshgrid(lons,lats)

    min_plot_value = -4
    max_plot_value = 4
    contour_interval = 1
    clevs = np.arange(min_plot_value,max_plot_value+1,contour_interval)

    plot = xr.plot.contourf(lon2d, lat2d, sfcT.isel(time=i), clevs, cmap=plt.cm.RdBu_r, 
                       extend='both', transform=ccrs.PlateCarree())

    resol = '50m'
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
      scale=resol, edgecolor='k', facecolor=cfeature.COLORS['land'])

    plt.savefig('./Figures/ArcticT2mAnom_'+str(time.values)+'.png', facecolor='white', dpi=300)

ValueError                                Traceback (most recent call last)
Input In [131], in <cell line: 1>()
     11 contour_interval = 1
     12 clevs = np.arange(min_plot_value,max_plot_value+1,contour_interval)
---> 14 plot = xr.plot.contourf(lon2d, lat2d, sfcT.isel(time=i),
     15                         clevs, cmap=plt.cm.RdBu_r, extend='both', transform=ccrs.PlateCarree())
     17 resol = '50m'  # use data at this scale
     18 land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
     19     scale=resol, edgecolor='k', facecolor=cfeature.COLORS['land'])

File ~/opt/anaconda3/lib/python3.8/site-packages/xarray/plot/plot.py:1132, in _plot2d.<locals>.newplotfunc(darray, x, y, figsize, size, aspect, ax, row, col, col_wrap, xincrease, yincrease, add_colorbar, add_labels, vmin, vmax, cmap, center, robust, extend, levels, infer_intervals, colors, subplot_kws, cbar_ax, cbar_kwargs, xscale, yscale, xticks, yticks, xlim, ylim, norm, **kwargs)
   1126 elif rgb is not None and not imshow_rgb:
   1127     raise ValueError(
   1128         'The "rgb" keyword is only valid for imshow()'
   1129         "with a three-dimensional array (per facet)"
   1130     )
-> 1132 xlab, ylab = _infer_xy_labels(
   1133     darray=darray, x=x, y=y, imshow=imshow_rgb, rgb=rgb
   1134 )
   1136 xval = darray[xlab]
   1137 yval = darray[ylab]

File ~/opt/anaconda3/lib/python3.8/site-packages/xarray/plot/utils.py:382, in _infer_xy_labels(darray, x, y, imshow, rgb)
    376 def _infer_xy_labels(darray, x, y, imshow=False, rgb=None):
    377     """
    378     Determine x and y labels. For use in _plot2d
    379 
    380     darray must be a 2 dimensional data array, or 3d for imshow only.
    381     """
--> 382     if (x is not None) and (x == y):
    383         raise ValueError("x and y cannot be equal.")
    385     if imshow and darray.ndim == 3:

File ~/opt/anaconda3/lib/python3.8/site-packages/xarray/core/common.py:136, in AbstractArray.__bool__(self)
    135 def __bool__(self: Any) -> bool:
--> 136     return bool(self.values)

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

Answer 1

我的代码在进行这些更改后可以工作：“for i in range(len(ds.time)):”和“sfcT.isel(time=i)”

ds=xr.open_mfdataset('./Data/ecmwf_dataset/t2m/*.grib', combine='nested', concat_dim='time',
                     engine='cfgrib',backend_kwargs={'indexpath': ''})

lons = ds['longitude'][:]
lats = ds['latitude'][:]
sfcT = ds['t2m']
time = ds['time'].dt.strftime('%m-%d-%Y')

for i in range(len(ds.time)):
    fig = plt.figure(figsize=(10,10))
    ax = fig.add_subplot(111, projection=ccrs.NorthPolarStereo(central_longitude=-47))
    ax.set_extent([-180, 180, 60, 90], crs=ccrs.PlateCarree())
    
    lon2d,lat2d = np.meshgrid(lons,lats)
 
    
    min_plot_value = -4
    max_plot_value = 4
    contour_interval = 1
    clevs = np.arange(min_plot_value,max_plot_value+1,contour_interval)

    plot = ax.contourf(lon2d, lat2d, sfcT.isel(time=i),
                            clevs, cmap=plt.cm.RdBu_r, extend='both', transform=ccrs.PlateCarree())

    resol = '50m'  # use data at this scale
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', \
        scale=resol, edgecolor='k', facecolor=cfeature.COLORS['land'])