Creating a bunch of figures with a Xarray multiple file dataset

Question

I have a bunch of .grib files that I would like to be able to read in and create figures for. I have the code to load one file and plot a figure, but I'd like to figure out the best way to loop it so I don't have to change the filename and other parameters for each file. So far, I've switched from using xr.open_dataset to xr.open_mfdataset. I've tried doing a For Loop using time, but I get an error at plot = ax.contourf saying, "TypeError: Input Z must be 2D, not 3D." The dataset involves latitude, longitude, t2m (2-m air temperature), and time. See below code and attached image.

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

My code works after making these changes: "for i in range(len(ds.time)):" and "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'])