Python绘制NOAA再分析NC文件报错‘Input z must be 2D, not 3D’求助

错误根源 #

你的omega变量是3D/4D维度（通常NOAA再分析数据的垂直速度变量维度为**[时间, 层次, 纬度, 经度]或[层次, 纬度, 经度]**），但Basemap.contourf()要求传入的变量必须是2D（仅纬度×经度的平面数据）。你当前写的omega[0,:,:]没有正确剥离多余维度，导致传入了3D数据。

修正步骤 #

1. 先确认变量维度结构 #

在读取omega后添加打印代码，明确数据结构：

print(omega.shape)
level = nc1.variables['level'][:]
print(level)  # 输出类似[1000, 850, 700, 500, 200]的气压层次值

比如输出(1,4,73,144)代表维度为[时间, 层次, 纬度, 经度]；输出(4,73,144)代表维度为[层次, 纬度, 经度]。

2. 提取500hPa对应的2D数据 #

找到500hPa在level数组中的索引，再根据维度结构提取平面数据：

• 如果是4D数据（时间+层次+经纬度）：

  level_index = np.where(level == 500)[0][0]
  omega_500 = omega[0, level_index, :, :]  # 取第1个时间步+500hPa层次
  

• 如果是3D数据（层次+经纬度）：

  level_index = np.where(level == 500)[0][0]
  omega_500 = omega[level_index, :, :]  # 直接取500hPa层次
  

3. 替换contourf的输入变量 #

将原代码中的omega[0,:,:]替换为处理后的omega_500：

cs1 = map.contourf(X,Y, omega_500, clevs1, cmap='RdBu_r')

完整修正代码 #

from netCDF4 import Dataset as NetCDFFile 
import matplotlib.pyplot as plt
import numpy as np
from mpl_toolkits.basemap import Basemap

fig = plt.figure(figsize=(10,9), constrained_layout=True ) 
ax = fig.add_subplot(111).set_title('b1' , size=10 , fontweight='bold')
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0 )

nc1 = NetCDFFile('C:/Users/user/Desktop/project/nc2024Habib/final-mjo/obs1-vv-ac.nc')
lat = nc1.variables['lat'][:]
lon = nc1.variables['lon'][:]
level = nc1.variables['level'][:]
omega = nc1.variables['omega'][:]

# 打印维度和层次值，确认数据结构
print("omega维度:", omega.shape)
print("气压层次:", level)

# 获取500hPa对应的索引
level_index = np.where(level == 500)[0][0]

# 根据维度提取2D数据
if len(omega.shape) == 4:
    omega_500 = omega[0, level_index, :, :]
elif len(omega.shape) == 3:
    omega_500 = omega[level_index, :, :]
else:
    raise ValueError("omega变量维度不符合预期，请检查NC文件结构")

map= Basemap(projection='merc',llcrnrlon=0.,llcrnrlat=10.,urcrnrlon=80.,urcrnrlat=45.,resolution='i',fix_aspect=False)  
map.drawcoastlines()
map.drawstates()
map.drawcountries()
map.drawlsmask(land_color='Linen', ocean_color='#CCFFFF')
map.drawcounties() 

parallels = np.arange(15.,45.,5.)
meridians = np.arange(30.,70.,5.)
map.drawparallels(parallels,labels=[1,0,0,0],fontsize=10 , linewidth=0.9,dashes=[4,9900], fontweight='bold' ,fontstyle='oblique')
map.drawmeridians(meridians,labels=[0,0,0,1],fontsize=10, linewidth=0.9,dashes=[4,9900] , fontweight='bold', fontstyle='oblique')

lons,lats= np.meshgrid(lon-180,lat )
X,Y = map(lons,lats)

clevs1 = np.arange(-0.02,0.01,0.005)
cs1 = map.contourf(X,Y, omega_500, clevs1,cmap='RdBu_r')
plt.show()

内容的提问来源于stack exchange，提问作者Sedighe