WRF风场后处理实战：用Python+Cartopy绘制500hPa风场矢量图（附完整代码）

WRF风场后处理实战用PythonCartopy绘制500hPa风场矢量图附完整代码气象数据分析中风场可视化是理解大气环流特征的关键环节。WRFWeather Research and Forecasting模式输出的数据包含丰富的三维风场信息但如何高效提取特定高度层如500hPa数据并实现专业级可视化是许多初学者面临的挑战。本文将手把手教你用Python的Cartopy地理信息库从WRF输出文件中提取500hPa风场数据绘制带中国地图的矢量风场图解决中文标注、坐标轴设置等实际问题。1. 环境准备与数据加载1.1 必需库安装确保已安装以下Python库可通过pip install命令安装pip install netCDF4 wrf-python cartopy matplotlib numpy1.2 关键库功能说明wrf-python专用于WRF数据提取的官方工具包cartopy地理空间数据处理与可视化netCDF4处理WRF输出的NetCDF格式文件1.3 数据文件结构典型的WRF输出文件如wrfout_d02_2023-06-01_18_00_00.nc包含以下关键变量import netCDF4 as nc ds nc.Dataset(wrfout_d02_2023-06-01_18_00_00.nc) print(ds.variables.keys()) # 查看所有变量2. 500hPa风场数据提取2.1 变量获取与插值WRF模式输出的是σ坐标下的三维数据需插值到等压面from wrf import getvar, interplevel # 获取基础变量 p getvar(ds, pressure) # 气压场hPa u getvar(ds, ua, unitsm/s) # 纬向风 v getvar(ds, va, unitsm/s) # 经向风 # 500hPa插值 u_500 interplevel(u, p, 500) v_500 interplevel(v, p, 500) lats, lons getvar(ds, lat), getvar(ds, lon)2.2 数据降采样技巧原始数据分辨率过高会导致箭头过密建议每5个网格点取一个样本skip 5 u_sub u_500[::skip, ::skip] v_sub v_500[::skip, ::skip] lons_sub lons[::skip, ::skip] lats_sub lats[::skip, ::skip]3. 地理信息可视化配置3.1 地图投影与范围设置Cartopy支持多种投影方式东亚区域推荐使用PlateCarree投影import cartopy.crs as ccrs import matplotlib.pyplot as plt proj ccrs.PlateCarree() fig plt.figure(figsize(12, 8)) ax fig.add_subplot(111, projectionproj) # 设置中国区域范围 extent [100, 125, 15, 40] # 经度西东纬度南北 ax.set_extent(extent, crsproj)3.2 地理要素添加# 添加海岸线、国界线等 ax.coastlines(resolution50m, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) # 中国省级行政区划需下载shapefile china cfeature.ShapelyFeature( Reader(Province_9.shp).geometries(), proj, edgecolork, facecolornone ) ax.add_feature(china, linewidth0.6)4. 矢量风场绘制进阶技巧4.1 箭头参数优化quiver ax.quiver( lons_sub, lats_sub, u_sub, v_sub, scale400, # 控制箭头长度 width0.003, # 箭头杆宽度 headwidth4, # 箭头头部宽度 headlength5, # 箭头头部长度 transformproj ) # 添加参考标尺 qk ax.quiverkey( quiver, 0.85, 0.05, 20, 20 m/s, labelposE, coordinatesaxes )4.2 色标映射风速大小将风速大小映射到颜色维度import numpy as np speed np.sqrt(u_sub**2 v_sub**2) quiv ax.quiver( lons_sub, lats_sub, u_sub, v_sub, speed, cmapjet, # 使用jet色带 scale400, width0.003, transformproj ) # 添加色标 cbar plt.colorbar(quiv, axax, shrink0.7) cbar.set_label(Wind Speed (m/s))5. 专业级地图修饰5.1 经纬度网格设置from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter gl ax.gridlines( draw_labelsTrue, linewidth1.2, colorgray, alpha0.5, linestyle-- ) gl.top_labels False # 关闭顶部标签 gl.right_labels False # 关闭右侧标签 # 设置刻度格式 ax.xaxis.set_major_formatter(LongitudeFormatter()) ax.yaxis.set_major_formatter(LatitudeFormatter())5.2 中文显示解决方案plt.rcParams[font.sans-serif] [SimHei] # 设置中文字体 plt.rcParams[axes.unicode_minus] False # 解决负号显示问题 ax.set_title(500hPa风场 - 2023年6月1日18时, fontsize14)6. 完整代码整合import numpy as np import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.io.shapereader import Reader from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter import matplotlib.pyplot as plt from wrf import getvar, interplevel from netCDF4 import Dataset # 数据加载 ds Dataset(wrfout_d02_2023-06-01_18_00_00.nc) p getvar(ds, pressure) u getvar(ds, ua, unitsm/s) v getvar(ds, va, unitsm/s) # 500hPa插值 u_500 interplevel(u, p, 500) v_500 interplevel(v, p, 500) lats, lons getvar(ds, lat), getvar(ds, lon) # 降采样 skip 5 u_sub u_500[::skip, ::skip] v_sub v_500[::skip, ::skip] lons_sub lons[::skip, ::skip] lats_sub lats[::skip, ::skip] speed np.sqrt(u_sub**2 v_sub**2) # 绘图设置 plt.rcParams[font.sans-serif] [SimHei] proj ccrs.PlateCarree() fig plt.figure(figsize(15, 10)) ax fig.add_subplot(111, projectionproj) # 地图范围与要素 extent [100, 125, 15, 40] ax.set_extent(extent, crsproj) ax.coastlines(resolution50m, linewidth0.8) ax.add_feature(cfeature.BORDERS, linestyle:, linewidth0.5) # 添加中国省界 china cfeature.ShapelyFeature( Reader(Province_9.shp).geometries(), proj, edgecolork, facecolornone ) ax.add_feature(china, linewidth0.6) # 绘制风场 quiv ax.quiver( lons_sub, lats_sub, u_sub, v_sub, speed, cmapjet, scale400, width0.003, transformproj ) # 色标与图例 cbar plt.colorbar(quiv, axax, shrink0.7) cbar.set_label(风速 (m/s)) qk ax.quiverkey( quiv, 0.85, 0.05, 20, 20 m/s, labelposE, coordinatesaxes ) # 网格与标题 gl ax.gridlines( draw_labelsTrue, linewidth1.2, colorgray, alpha0.5, linestyle-- ) gl.top_labels gl.right_labels False ax.xaxis.set_major_formatter(LongitudeFormatter()) ax.yaxis.set_major_formatter(LatitudeFormatter()) ax.set_title(500hPa风场 - 2023年6月1日18时, fontsize16) plt.tight_layout() plt.savefig(wrf_500hPa_wind.png, dpi300, bbox_inchestight) plt.show()7. 常见问题解决方案7.1 中文显示异常若中文仍显示为方框确认系统已安装中文字体如SimHei指定字体路径import matplotlib.font_manager as fm font_path /path/to/SimHei.ttf font_prop fm.FontProperties(fnamefont_path) ax.set_title(标题, fontpropertiesfont_prop)7.2 箭头密度调整根据区域大小调整skip参数大区域分析skip10小区域精细分析skip37.3 性能优化建议对于大型WRF输出文件使用xarray替代netCDF4加速数据读取先计算风速再降采样speed np.sqrt(u**2 v**2)[::skip, ::skip]通过本教程你将掌握WRF风场数据可视化的核心技能。实际应用中可进一步尝试叠加等高线显示位势高度场添加天气系统标注如低压中心制作时间序列动画