### 处理ERA5再分析数据并绘制相对涡度平流图 为了实现对ERA5再分析数据的后处理以及绘制相对涡度平流图，可以按照以下方式操作： #### 数据准备 首先需要确保已经成功下载了所需的ERA5数据。如果采用的是单次下载的方式，则可以通过脚本优化为批量下载模式[^1]。 对于ERA5数据的读取和处理，通常会使用 `xarray` 和 `netCDF4` 库来加载 `.nc` 文件格式的数据集。这些库能够高效地解析 NetCDF 格式的气象数据。 #### 后处理流程 以下是具体的操作步骤说明及其对应的 Python 实现代码片段： ##### 加载所需库 ```python import xarray as xr import numpy as np from scipy.ndimage import uniform_filter import matplotlib.pyplot as plt ``` ##### 读取 ERA5 数据 假设已有的文件名为 `era5_data.nc`，可以直接利用 `xarray.open_dataset()` 方法打开该文件。 ```python data = xr.open_dataset('era5_data.nc') u_component_wind = data['u'] # 假设 u 风分量变量名是 'u' v_component_wind = data['v'] # 假设 v 风分量变量名是 'v' geopotential_height = data['z'] # 地转位势高度 (可选用于计算其他参数) ``` ##### 计算相对涡度 相对涡度 ζ 的定义涉及风场的空间导数项 ∂v/∂x − ∂u/∂y 。这里我们借助中心差分法近似求解偏微分方程。 ```python def calculate_relative_vorticity(u, v, dx=7e3, dy=7e3): """dx 和 dy 是网格间距，默认单位米""" dudy = (u.shift(lat=-1) - u.shift(lat=1)) / (2 * dy) dvdx = (v.shift(lon=-1) - v.shift(lon=1)) / (2 * dx) relative_vorticity = dvdx - dudy return relative_vorticity.fillna(0).values # 替换可能存在的 NaN 值以防绘图错误 relative_vorticity = calculate_relative_vorticity(u_component_wind.mean(dim='time'), v_component_wind.mean(dim='time')) ``` ##### 平滑滤波 由于原始格点分辨率较高可能导致噪声干扰图像质量，因此建议应用简单的均匀卷积核进行空间平均过滤。 ```python smoothed_rvor = uniform_filter(relative_vorticity, size=(3, 3)) ``` ##### 绘制矢量叠加地图背景下的相对涡度分布 最后一步就是把结果可视化出来啦！下面展示了一个基本的地图配置方案。 ```python plt.figure(figsize=(8,6)) lon = u_component_wind.lon.values lat = u_component_wind.lat.values X,Y=np.meshgrid(lon,lat) cmap=plt.cm.RdBu_r levels=np.linspace(-1.e-4,1.e-4,num=9) cf=plt.contourf(X,Y,smoothed_rvor*1e5,cmap=cmap, levels=levels, extend="both") cb=plt.colorbar(cf); cb.set_label(r"$\times 10^{-5} s^{-1}$", fontsize=12) quiver_skip={'lon':slice(None,None,4),'lat':slice(None,None,4)} QV=plt.quiver( X[::4][::4], Y[::4][::4], u_component_wind.isel(time=0)[quiver_skip].mean(dim='level').values, v_component_wind.isel(time=0)[quiver_skip].mean(dim='level').values, scale=500,scale_units='inches',color='k') plt.title("Relative Vorticity Advection with Wind Field Overlay") plt.xlabel("Longitude"); plt.ylabel("Latitude"); plt.show() ``` 以上即完成了从数据导入至最终图形输出的一整套工作链条描述。