如何使用ERA5数据绘制日平均风场及解决netCDF4幂运算报错问题

你当前触发的类型错误本质是：netCDF4.Variable是netCDF库定义的特殊数据容器，不支持直接进行幂运算（**）。你已经通过切片操作U2M[:]、V2M[:]将原始变量读取为numpy数组，且完成了缺省值替换存储在U2M_nans、V2M_nans中，计算风速时直接使用这两个numpy数组即可解决报错。
另外补充注意点：你使用的MERRA2_400.inst3_3d_asm_Nv是三维大气同化产品，U、V变量维度为[时间, 气压层, 纬度, 经度]，如果你需要计算2米风，需要改用对应二维地表产品；如果要计算特定高度层的风场，需要先从lev维度选取对应层的数据。

from netCDF4 import Dataset
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.ticker as mticker

# 读取nc数据
data = Dataset(r'C:/Users/MERRA2_400.inst3_3d_asm_Nv.20200101.nc4', mode='r')
# 经纬度处理
lons = data.variables['lon'][:]
lats = data.variables['lat'][:]
lon, lat = np.meshgrid(lons, lats)

# 读取U、V变量，示例取最低层（近地面）气压层的数据，索引可根据需求调整
U = data.variables['U']
V = data.variables['V']
U_data = U[:]
V_data = V[:]
# 替换缺省值为nan
U_data[U_data == U._FillValue] = np.nan
V_data[V_data == V._FillValue] = np.nan
# 选取近地面层（lev维度索引为-1，即最底层）
U_ground = U_data[:, -1, :, :]
V_ground = V_data[:, -1, :, :]

# 计算风速，使用numpy数组计算，避免类型错误
ws = np.sqrt(U_ground**2 + V_ground**2)
# 计算日平均
ws_daily_avg = np.nanmean(ws, axis=0)
u_daily_avg = np.nanmean(U_ground, axis=0)
v_daily_avg = np.nanmean(V_ground, axis=0)

# 可视化部分
fig = plt.figure(figsize=(12,8))
ax = plt.axes(projection=ccrs.PlateCarree())
# 绘制海岸线
ax.coastlines(resolution='50m', linewidth=0.8)
# 绘制风速填色图
contourf = ax.contourf(lon, lat, ws_daily_avg, levels=20, cmap='RdYlBu_r', transform=ccrs.PlateCarree())
# 每10个格点绘制一个风矢，避免图面过密
skip = 10
ax.quiver(lon[::skip, ::skip], lat[::skip, ::skip], 
          u_daily_avg[::skip, ::skip], v_daily_avg[::skip, ::skip],
          transform=ccrs.PlateCarree(), color='k', scale=150)
# 添加色标
plt.colorbar(contourf, shrink=0.8, label='平均风速 (m/s)')
# 添加格网和刻度
gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.5, linestyle='--', color='gray', alpha=0.7)
gl.top_labels = False
gl.right_labels = False
gl.xformatter = LONGITUDE_FORMATTER
gl.yformatter = LATITUDE_FORMATTER
gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 30))
gl.ylocator = mticker.FixedLocator(np.arange(-90, 91, 30))
# 添加标题
plt.title('2020年1月1日日平均近地面风场', fontsize=14)
plt.show()

# 关闭数据集
data.close()

• 如果你使用的是专门的2米风MERRA2产品，U、V变量维度为[时间, 纬度, 经度]，不需要选取气压层，直接使用读取的数组计算即可
• 风矢的skip参数和scale参数可根据你的研究区域和风速大小调整，避免风矢过密或者过长影响可读性

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