# Python+WRF风场可视化实战：从nc文件提取500hPa风速到生成专业气象图 如果你正在处理WRF模式输出的netCDF数据，并且希望将那些复杂的多维数组转化为清晰、专业、可直接用于学术报告或论文发表的气象图，那么这篇文章就是为你准备的。无论是刚接触气象数据可视化的研究生，还是需要优化现有工作流的科研人员，我们都会深入探讨如何利用Python生态中的xarray和cartopy库，将原始的格点数据转化为具有学术美感的矢量风场图。这个过程不仅仅是调用几个绘图函数，更涉及到数据插值、坐标转换、图形美学以及性能优化等一系列实战技巧。我遇到过不少同行，他们能跑出漂亮的模拟结果，却在最后的数据呈现环节卡壳，要么图形粗糙不符合出版要求，要么代码冗长难以维护。接下来，我将分享一套经过多个项目验证的、高效且优雅的解决方案。 ## 1. 理解WRF输出结构与500hPa风场数据提取 在动手写代码之前，我们必须先搞清楚WRF的netCDF文件里到底装了些什么。WRF的输出是一个典型的多维时空数据集，其结构远比普通的二维表格复杂。一个标准的`wrfout`文件通常包含数十个变量，从基本的气压、温度、湿度，到三维风场（U、V、W分量），这些变量被组织在时间、垂直层、南北格点、东西格点这四个维度上。对于风场可视化，我们最常关心的是**U（东西风分量）**和**V（南北风分量）**。 > **注意**：WRF模式中的风变量`UA`和`VA`代表的是**在模式层上的笛卡尔风分量**，而`U10`、`V10`则是10米高度的风。在进行等压面插值前，务必确认你使用的是正确的高度层风变量。 提取特定等压面（如500hPa）的数据，本质是一个**垂直插值**过程。因为WRF的输出是基于地形追随坐标（Sigma坐标）的，其垂直层是等间距的Sigma层，而非我们气象分析中常用的等压面。这就需要用到`wrf-python`这个专门为WRF数据后处理设计的工具包，它封装了复杂的坐标转换和插值算法。 首先，我们来看一个基础的数据读取和变量提取示例。这里假设你已经有了一个名为`wrfout_d01_2023-07-01_00:00:00`的典型输出文件。 ```python import xarray as xr from wrf import getvar, interplevel, to_np, latlon_coords import numpy as np # 使用xarray打开netCDF文件，它比netCDF4库提供了更友好、更强大的接口 ds = xr.open_dataset('wrfout_d01_2023-07-01_00:00:00', engine='netcdf4') # 提取气压场和风场变量 # getvar函数是wrf-python的核心，它能智能识别WRF特有的变量名和属性 pressure = getvar(ds, "pressure") # 三维气压场，单位：hPa u_wind = getvar(ds, "ua", units="m s-1") # 东西风分量，转换为米/秒 v_wind = getvar(ds, "va", units="m s-1") # 南北风分量，转换为米/秒 # 获取经纬度坐标 # latlon_coords能根据数据自动计算每个格点的经纬度，比从文件属性中读取更准确 lats, lons = latlon_coords(pressure) # 关键步骤：垂直插值到500hPa等压面 # interplevel函数执行三维线性插值，返回指定气压层上的二维数据 u_500 = interplevel(u_wind, pressure, 500) v_500 = interplevel(v_wind, pressure, 500) # 计算全风速（标量），用于后续填色或筛选 wspd_500 = np.sqrt(u_500**2 + v_500**2) print(f"数据维度: 经度={lons.shape}, 纬度={lats.shape}") print(f"500hPa U分量范围: [{u_500.min().values:.2f}, {u_500.max().values:.2f}] m/s") print(f"500hPa V分量范围: [{v_500.min().values:.2f}, {v_500.max().values:.2f}] m/s") ``` 这段代码构成了数据处理的基石。但实际操作中，你可能会遇到一些棘手的情况。例如，WRF输出可能包含多个时间步长，而`interplevel`默认只处理单个时间步。如果你的文件是时间序列，就需要增加一个时间维度循环。另一个常见问题是数据区域太大，导致全分辨率的风矢量过于密集，图形杂乱。这时就需要进行**数据抽稀**。 ```python # 处理多时间步长数据的示例 # 假设我们只处理第一个时间步（time=0） time_idx = 0 pressure_t = pressure.isel(Time=time_idx) u_wind_t = u_wind.isel(Time=time_idx) v_wind_t = v_wind.isel(Time=time_idx) # 然后进行插值 u_500_t = interplevel(u_wind_t, pressure_t, 500) v_500_t = interplevel(v_wind_t, pressure_t, 500) # 数据抽稀：每隔5个格点取一个点，使箭头密度适中 stride = 5 u_sub = u_500_t[::stride, ::stride] v_sub = v_500_t[::stride, ::stride] lons_sub = lons[::stride, ::stride] lats_sub = lats[::stride, ::stride] print(f"抽稀后格点数: {lons_sub.size} (原数据: {lons.size})") ``` 理解数据结构和掌握正确的提取方法是后续所有可视化工作的前提。错误的数据会导致错误的图形，无论你的绘图技巧多么高超。 ## 2. 构建专业气象地图：Cartopy核心配置与投影选择 有了数据，下一步就是创建地图底图。在气象领域，`cartopy`库已经成为事实上的标准，因为它原生支持多种地图投影和地理特征数据，并且与matplotlib无缝集成。选择合适的地图投影对于准确表达气象信息至关重要，不同的投影会扭曲距离、面积或方向。 对于中纬度天气系统分析，**兰勃特等角投影（Lambert Conformal Conic）** 是最常用的选择之一，因为它能较好地保持中纬度地区的形状和方向，特别适合展示西风带中的槽脊系统。而对于全球或热带地区的研究，**圆柱投影（如PlateCarree）** 或**麦卡托投影（Mercator）** 可能更合适。在WRF社区中，由于模式本身常使用兰勃特投影，因此绘图时也常保持一致。 下面我们创建一个具有专业气象地图所有要素的底图：海岸线、国界、经纬网格、以及符合出版要求的刻度标注。 ```python import cartopy.crs as ccrs import cartopy.feature as cfeature from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER import matplotlib.pyplot as plt import matplotlib.ticker as mticker # 创建图形和坐标轴，指定投影 # 这里使用PlateCarree（等经纬度投影）作为示例，实际应根据你的数据区域选择 fig = plt.figure(figsize=(14, 8), dpi=150) ax = plt.axes(projection=ccrs.PlateCarree()) # 设置图形显示范围（经度起始，经度结束，纬度起始，纬度结束） # 这个范围应该覆盖你的数据区域，并留出适当边界 extent = [105, 125, 20, 40] # 示例：东亚区域 ax.set_extent(extent, crs=ccrs.PlateCarree()) # 添加地理特征 - 这是让地图“活”起来的关键 # 分辨率选择：'110m'（粗）、'50m'（中）、'10m'（细），出版级建议用'50m'或'10m' ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8, edgecolor='black') ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5, linestyle='--', edgecolor='gray') ax.add_feature(cfeature.OCEAN, alpha=0.2) # 海洋浅色填充 ax.add_feature(cfeature.LAND, alpha=0.1) # 陆地浅色填充 # 添加河流、湖泊等（可选） ax.add_feature(cfeature.RIVERS.with_scale('50m'), linewidth=0.5, edgecolor='blue', alpha=0.5) ax.add_feature(cfeature.LAKES.with_scale('50m'), linewidth=0.5, edgecolor='blue', alpha=0.2) # 配置经纬网格线 - 学术图的精髓所在 gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=True, linewidth=0.8, color='gray', alpha=0.5, linestyle=':') # 控制哪些边显示标签 gl.top_labels = False gl.right_labels = False gl.left_labels = True gl.bottom_labels = True # 设置经纬度格式（度°格式） gl.xformatter = LONGITUDE_FORMATTER gl.yformatter = LATITUDE_FORMATTER # 自定义刻度间隔和样式 gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 5)) gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 5)) # 设置刻度标签字体大小 gl.xlabel_style = {'size': 10, 'color': 'black'} gl.ylabel_style = {'size': 10, 'color': 'black'} plt.title('专业气象地图底图示例', fontsize=14, pad=15) plt.show() ``` 对于WRF数据，特别是当你的模拟使用了特定的地图投影时，直接使用该投影绘图可以避免额外的坐标转换。`wrf-python`提供了`get_cartopy`函数来获取WRF数据自身的投影信息： ```python from wrf import get_cartopy, cartopy_xlim, cartopy_ylim # 从WRF数据中获取其使用的投影 cart_proj = get_cartopy(u_500) fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=cart_proj) # 使用WRF数据自身的投影 # 自动获取数据在投影坐标系下的范围 ax.set_xlim(cartopy_xlim(u_500)) ax.set_ylim(cartopy_ylim(u_500)) # 添加特征和网格（同上） # ... ``` 地图投影的选择和配置直接影响到图形的专业程度。一个常见的错误是使用不合适的投影导致风矢量的方向失真，或者网格线过于密集影响阅读。我建议在项目初期就确定好投影方式，并保持整个研究中的所有图形一致。 ## 3. 风场可视化进阶：矢量箭头优化与叠加等值线 现在进入最核心的部分：将风场数据可视化到地图上。单纯绘制箭头很简单，但要绘制出既美观又信息丰富的风场图，需要掌握一系列技巧。我们将重点讨论三个关键方面：**箭头样式优化**、**颜色映射选择**和**多变量叠加**。 ### 3.1 矢量箭头（Quiver）的精细控制 matplotlib的`quiver`函数是绘制矢量的基础工具，但它有十几个参数需要调整才能达到理想效果。下面是一个经过优化的配置示例： ```python # 继续使用之前的数据和地图 # 计算风速大小，用于颜色映射 wind_speed = np.sqrt(u_sub**2 + v_sub**2) # 创建带颜色的风矢量图 # 关键参数说明： # scale: 箭头缩放因子，值越大箭头越短，需要根据你的风速范围调整 # width: 箭头杆的宽度 # headwidth: 箭头头部的宽度（相对于杆宽） # headlength: 箭头头部的长度 # headaxislength: 箭头头部轴长 # angles: 指定输入数据的角度格式，'uv'表示直接使用U/V分量 # pivot: 箭头旋转支点，'tail'在尾部，'mid'在中部，'tip'在头部 # cmap: 颜色映射，这里使用'viridis'，也可用'plasma', 'inferno'等 quiv = ax.quiver(lons_sub, lats_sub, u_sub, v_sub, wind_speed, # 用风速着色 cmap='viridis', scale=200, # 这个值需要根据你的数据调整 width=0.003, headwidth=3, headlength=4, headaxislength=3.5, angles='uv', pivot='mid', # 箭头从中点开始绘制，看起来更平衡 transform=ccrs.PlateCarree()) # 添加颜色条 cbar = plt.colorbar(quiv, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8) cbar.set_label('风速 (m/s)', fontsize=11) cbar.ax.tick_params(labelsize=10) # 添加参考箭头（图例） # 在图形角落添加一个参考箭头，标注其代表的风速值 ax.quiverkey(quiv, 0.95, 0.05, 20, '20 m/s', labelpos='E', coordinates='axes', color='black', fontproperties={'size': 10}) ``` > **提示**：`scale`参数是最需要经验调整的。一个实用的技巧是：先计算你数据中最大风速对应的箭头长度，然后调整`scale`值，使得最长的箭头不超过你希望的长度。可以尝试这个公式作为起点：`scale = 最大风速 / 期望箭头长度（以度为单位）`。 ### 3.2 叠加等值线与填色图 在实际分析中，我们经常需要将风场与其他气象要素叠加，比如等压面高度场、涡度场、垂直速度等。这能提供更全面的天气系统信息。下面演示如何将500hPa风场与位势高度场叠加： ```python # 提取500hPa位势高度场 # 假设我们已经有了高度场变量z（单位：位势米） z = getvar(ds, "z", units="dm") # 转换为位势什米（常用单位） z_500 = interplevel(z, pressure, 500) # 创建图形 fig = plt.figure(figsize=(15, 10)) ax = plt.axes(projection=cart_proj) ax.set_extent(extent, crs=ccrs.PlateCarree()) # 添加地理特征（同上，略） # 绘制位势高度等值线 # 选择适当的等高线间隔，500hPa高度场通常在5520-5880位势米之间 levels = np.arange(5520, 5880, 60) # 每60位势米一条线 contours = ax.contour(lons, lats, z_500, levels=levels, colors='black', linewidths=1.2, transform=ccrs.PlateCarree()) # 在等高线上标注数值 plt.clabel(contours, inline=True, fontsize=9, fmt='%d') # 绘制风速填色图（背景） # 选择风速分级和颜色 wspd_levels = np.arange(10, 70, 5) # 从10到70 m/s，每5 m/s一级 cf = ax.contourf(lons, lats, wind_speed, levels=wspd_levels, cmap='YlOrRd', # 黄-橙-红渐变色，适合风速 alpha=0.7, # 稍微透明，避免完全遮盖等高线 extend='both', # 延伸颜色条两端 transform=ccrs.PlateCarree()) # 添加风速填色的颜色条 cbar_wspd = plt.colorbar(cf, ax=ax, orientation='horizontal', pad=0.08, shrink=0.8, aspect=30) cbar_wspd.set_label('500 hPa 风速 (m/s)', fontsize=12) # 最后绘制风矢量（在填色和等高线之上） # 这里使用黑白箭头，避免与填色图颜色冲突 quiv = ax.quiver(lons_sub, lats_sub, u_sub, v_sub, color='black', scale=250, width=0.0035, headwidth=3.5, headlength=4.5, transform=ccrs.PlateCarree()) # 添加风矢量图例 ax.quiverkey(quiv, 0.15, 0.92, 25, '25 m/s', labelpos='E', coordinates='axes', color='black', fontproperties={'size': 11, 'weight': 'bold'}) plt.title('500 hPa 位势高度（黑色等值线）与风场（箭头）叠加图\n填色表示风速大小', fontsize=14, pad=20) plt.tight_layout() ``` 这种多层叠加的图形能同时展示多个气象要素，但需要注意视觉层次的安排。通常的顺序是：先绘制填色背景，然后是等值线，最后是矢量箭头。颜色选择也要避免冲突，确保所有元素都清晰可辨。 ### 3.3 处理高密度数据的性能优化 当处理高分辨率WRF输出（如1km网格）时，直接绘制所有矢量点会导致图形过于密集且渲染缓慢。除了前面提到的数据抽稀，还有几种优化策略： **策略一：基于风速的智能抽稀** - 只在风速较大的区域密集绘制，弱风区稀疏绘制： ```python # 创建风速掩码 wind_threshold = 15 # 风速阈值，单位m/s high_wind_mask = wind_speed > wind_threshold # 对高风速区域使用更密集的抽稀 stride_high = 2 stride_low = 8 # 分别处理高风速和低风速区域 lons_high = lons[high_wind_mask][::stride_high] lats_high = lats[high_wind_mask][::stride_high] u_high = u_500.values[high_wind_mask][::stride_high] v_high = v_500.values[high_wind_mask][::stride_high] lons_low = lons[~high_wind_mask][::stride_low] lats_low = lats[~high_wind_mask][::stride_low] u_low = u_500.values[~high_wind_mask][::stride_low] v_low = v_500.values[~high_wind_mask][::stride_low] # 合并数据 lons_combined = np.concatenate([lons_high, lons_low]) lats_combined = np.concatenate([lats_high, lats_low]) u_combined = np.concatenate([u_high, u_low]) v_combined = np.concatenate([v_high, v_low]) print(f"优化后数据点: {len(lons_combined)} (原数据: {lons.size})") ``` **策略二：使用流线（Streamplot）替代箭头** - 对于某些分析，流线能更好地展示气流走向： ```python # 流线图适合展示气流整体形态，但不适合精确表示风速大小 stride = 10 # 需要更大的抽稀步长 lons_stream = lons[::stride, ::stride] lats_stream = lats[::stride, ::stride] u_stream = u_500.values[::stride, ::stride] v_stream = v_500.values[::stride, ::stride] # 计算风速用于着色 wspd_stream = np.sqrt(u_stream**2 + v_stream**2) # 绘制流线 stream = ax.streamplot(lons_stream, lats_stream, u_stream, v_stream, color=wspd_stream, # 用风速着色 cmap='coolwarm', linewidth=1.5, density=2.5, # 流线密度 arrowsize=1.2, # 箭头大小 arrowstyle='->', transform=ccrs.PlateCarree()) # 添加颜色条 cbar = plt.colorbar(stream.lines, ax=ax, orientation='vertical', pad=0.02) cbar.set_label('风速 (m/s)', fontsize=11) ``` 选择哪种可视化方式取决于你的分析目的。如果关注**精确的风速和风向**，使用矢量箭头；如果关注**气流的整体形态和路径**，流线图可能更合适。 ## 4. 学术出版级图形输出与自动化脚本 当图形在屏幕上看起来不错后，最后一步是将其输出为适合学术出版的高质量文件。这不仅仅是点击"保存"那么简单，需要考虑文件格式、分辨率、字体嵌入、颜色模式等多个细节。 ### 4.1 图形输出参数详解 ```python # 在绘图代码的最后，添加保存命令 output_filename = '500hPa_wind_field_analysis.png' # 保存图形的关键参数 plt.savefig(output_filename, dpi=600, # 出版要求通常为300-600 dpi bbox_inches='tight', # 裁剪空白边缘 pad_inches=0.05, # 保留少量内边距 facecolor='white', # 确保背景为白色 edgecolor='none', # 无边框 transparent=False) # 不透明背景 # 也可以保存为矢量格式，适合论文中的线图和地图 output_pdf = '500hPa_wind_field_analysis.pdf' plt.savefig(output_pdf, format='pdf', bbox_inches='tight', pad_inches=0.05) print(f"图形已保存为: {output_filename} 和 {output_pdf}") ``` 对于学术出版，我强烈建议同时保存为**PNG（用于预览和演示）**和**PDF/EPS（用于论文投稿）**两种格式。PDF/EPS是矢量格式，无论放大多少倍都不会失真，特别适合包含文字和线条的地图。 ### 4.2 字体配置与样式统一 学术期刊通常对图形中的字体有明确要求，最常见的是**Times New Roman**或**Arial**。在Python中配置字体需要一些技巧： ```python import matplotlib from matplotlib import rcParams # 方法1：直接指定字体文件路径（最可靠） font_path = '/usr/share/fonts/truetype/times.ttf' # Linux系统Times字体路径 # Windows: 'C:/Windows/Fonts/times.ttf' # Mac: '/System/Library/Fonts/Times.ttc' matplotlib.font_manager.fontManager.addfont(font_path) font_name = matplotlib.font_manager.FontProperties(fname=font_path).get_name() rcParams['font.family'] = font_name rcParams['font.size'] = 11 # 方法2：使用系统已安装的字体（如果知道确切名称） rcParams['font.family'] = 'Times New Roman' rcParams['mathtext.fontset'] = 'stix' # 数学字体使用STIX，与Times风格一致 # 统一图形样式 rcParams['axes.labelsize'] = 12 rcParams['axes.titlesize'] = 14 rcParams['xtick.labelsize'] = 10 rcParams['ytick.labelsize'] = 10 rcParams['legend.fontsize'] = 10 rcParams['figure.titlesize'] = 16 # 确保负号正常显示 rcParams['axes.unicode_minus'] = False ``` ### 4.3 创建可复用的自动化脚本 在实际研究中，我们通常需要处理多个时间步、多个气压层或多个模拟案例。手动调整每个图形是不现实的。下面是一个模块化的脚本框架，你可以根据自己的需求扩展： ```python """ WRF风场可视化自动化脚本 作者：[你的名字] 日期：2024年 功能：批量处理WRF输出文件，生成500hPa风场分析图 """ import os import glob import xarray as xr import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from wrf import getvar, interplevel, to_np, latlon_coords, get_cartopy from matplotlib import rcParams # 配置字体和样式 def setup_plot_style(): """配置出版级图形样式""" rcParams['font.family'] = 'Times New Roman' rcParams['font.size'] = 11 rcParams['axes.labelsize'] = 12 rcParams['axes.titlesize'] = 14 rcParams['xtick.labelsize'] = 10 rcParams['ytick.labelsize'] = 10 rcParams['legend.fontsize'] = 10 rcParams['figure.titlesize'] = 16 rcParams['axes.unicode_minus'] = False rcParams['savefig.dpi'] = 600 rcParams['savefig.bbox'] = 'tight' class WRFWindVisualizer: """WRF风场可视化器""" def __init__(self, wrf_file, output_dir='./figures'): """ 初始化可视化器 参数: wrf_file: WRF输出文件路径 output_dir: 输出目录 """ self.wrf_file = wrf_file self.output_dir = output_dir self.ds = None self.u_500 = None self.v_500 = None self.lats = None self.lons = None # 创建输出目录 os.makedirs(output_dir, exist_ok=True) def load_and_process(self, time_idx=0, pressure_level=500): """ 加载并处理WRF数据 参数: time_idx: 时间索引（对于多时间步文件） pressure_level: 气压层（hPa） """ print(f"处理文件: {self.wrf_file}") # 打开数据集 self.ds = xr.open_dataset(self.wrf_file, engine='netcdf4') # 提取变量 pressure = getvar(self.ds, "pressure") u_wind = getvar(self.ds, "ua", units="m s-1") v_wind = getvar(self.ds, "va", units="m s-1") # 如果有多时间步，选择特定时间 if 'Time' in pressure.dims: pressure = pressure.isel(Time=time_idx) u_wind = u_wind.isel(Time=time_idx) v_wind = v_wind.isel(Time=time_idx) # 获取经纬度 self.lats, self.lons = latlon_coords(pressure) # 垂直插值到指定气压层 self.u_500 = interplevel(u_wind, pressure, pressure_level) self.v_500 = interplevel(v_wind, pressure, pressure_level) print(f"数据维度: {self.u_500.shape}") print(f"U分量范围: [{self.u_500.min().values:.2f}, {self.u_500.max().values:.2f}] m/s") def create_wind_map(self, extent=None, stride=5, add_geopotential=False, add_contours=True): """ 创建风场地图 参数: extent: 地图范围 [lon_min, lon_max, lat_min, lat_max] stride: 数据抽稀步长 add_geopotential: 是否添加位势高度场 add_contours: 是否添加等高线 """ if self.u_500 is None: raise ValueError("请先调用load_and_process()方法加载数据") # 数据抽稀 u_sub = self.u_500.values[::stride, ::stride] v_sub = self.v_500.values[::stride, ::stride] lons_sub = self.lons[::stride, ::stride] lats_sub = self.lats[::stride, ::stride] # 计算风速 wind_speed = np.sqrt(u_sub**2 + v_sub**2) # 创建图形 fig = plt.figure(figsize=(14, 10)) # 获取WRF投影或使用默认投影 try: cart_proj = get_cartopy(self.u_500) ax = plt.axes(projection=cart_proj) # 自动设置范围 ax.set_xlim(cartopy_xlim(self.u_500)) ax.set_ylim(cartopy_ylim(self.u_500)) except: # 如果获取投影失败，使用等经纬度投影 ax = plt.axes(projection=ccrs.PlateCarree()) if extent: ax.set_extent(extent, crs=ccrs.PlateCarree()) else: # 自动计算数据范围 lon_min, lon_max = self.lons.min(), self.lons.max() lat_min, lat_max = self.lats.min(), self.lats.max() ax.set_extent([lon_min, lon_max, lat_min, lat_max], crs=ccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale('50m'), linewidth=0.8) ax.add_feature(cfeature.BORDERS.with_scale('50m'), linewidth=0.5, linestyle='--', alpha=0.7) ax.add_feature(cfeature.OCEAN, alpha=0.1) ax.add_feature(cfeature.LAND, alpha=0.1) # 添加位势高度场（如果要求） if add_geopotential: z = getvar(self.ds, "z", units="dm") if 'Time' in z.dims: z = z.isel(Time=0) z_500 = interplevel(z, getvar(self.ds, "pressure"), 500) # 绘制等高线 levels = np.arange(5200, 5900, 60) contours = ax.contour(self.lons, self.lats, z_500, levels=levels, colors='black', linewidths=1.0, transform=ccrs.PlateCarree()) plt.clabel(contours, inline=True, fontsize=9, fmt='%d') # 绘制风矢量 quiv = ax.quiver(lons_sub, lats_sub, u_sub, v_sub, wind_speed, cmap='viridis', scale=200, width=0.003, headwidth=3, headlength=4, pivot='mid', transform=ccrs.PlateCarree()) # 添加颜色条 cbar = plt.colorbar(quiv, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8, aspect=30) cbar.set_label('风速 (m/s)', fontsize=12) # 添加参考箭头 ax.quiverkey(quiv, 0.15, 0.92, 20, '20 m/s', labelpos='E', coordinates='axes', fontproperties={'size': 11, 'weight': 'bold'}) # 添加网格线 gl = ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle=':') gl.top_labels = False gl.right_labels = False gl.xlabel_style = {'size': 10} gl.ylabel_style = {'size': 10} # 设置标题 from datetime import datetime if 'Times' in self.ds.variables: time_str = str(self.ds.Times.values[0])[:19] title = f'500 hPa 风场 - {time_str}' else: title = '500 hPa 风场' plt.title(title, fontsize=16, pad=20) return fig, ax def save_figure(self, fig, filename_prefix='wind_field'): """保存图形""" timestamp = datetime.now().strftime('%Y%m%d_%H%M%S') filename_png = f"{filename_prefix}_{timestamp}.png" filename_pdf = f"{filename_prefix}_{timestamp}.pdf" filepath_png = os.path.join(self.output_dir, filename_png) filepath_pdf = os.path.join(self.output_dir, filename_pdf) fig.savefig(filepath_png, dpi=600, bbox_inches='tight') fig.savefig(filepath_pdf, bbox_inches='tight') print(f"图形已保存: {filepath_png}") print(f"图形已保存: {filepath_pdf}") plt.close(fig) def process_batch(self, file_pattern, pressure_levels=[500, 700, 850]): """ 批量处理多个文件 参数: file_pattern: 文件匹配模式（如 'wrfout_d01_*'） pressure_levels: 要处理的气压层列表 """ files = sorted(glob.glob(file_pattern)) print(f"找到 {len(files)} 个文件") for file_idx, wrf_file in enumerate(files): print(f"\n处理文件 {file_idx+1}/{len(files)}: {os.path.basename(wrf_file)}") # 重新初始化处理每个文件 self.__init__(wrf_file, self.output_dir) self.load_and_process(time_idx=0) for level in pressure_levels: print(f" 处理 {level} hPa 层...") # 这里可以扩展为处理不同气压层 # 实际实现需要修改load_and_process方法以接受level参数 # 创建并保存图形 fig, _ = self.create_wind_map(stride=6, add_geopotential=True) self.save_figure(fig, f'wrf_wind_{file_idx:03d}') # 使用示例 if __name__ == "__main__": # 设置图形样式 setup_plot_style() # 创建可视化器实例 visualizer = WRFWindVisualizer( wrf_file='wrfout_d01_2023-07-01_00:00:00', output_dir='./output_figures' ) # 处理单个文件 visualizer.load_and_process(time_idx=0, pressure_level=500) fig, ax = visualizer.create_wind_map( extent=[110, 120, 25, 35], # 自定义范围 stride=6, add_geopotential=True ) visualizer.save_figure(fig, '500hPa_wind_analysis') # 或者批量处理 # visualizer.process_batch('wrfout_d01_2023-07-*.nc') ``` 这个类封装了完整的处理流程，从数据加载、处理到图形生成和保存。你可以轻松地修改参数来适应不同的需求，比如改变气压层、调整图形样式或添加新的分析要素。 ### 4.4 常见问题与调试技巧 在实际使用中，你可能会遇到各种问题。下面是一些常见问题及其解决方案： **问题1：箭头太密集或太稀疏** > 调整`stride`参数和`scale`参数。可以先尝试`stride=10`和`scale=200`，然后根据效果微调。记住，`scale`值越大，箭头越短。 **问题2：图形保存后质量差** > 确保`dpi`参数足够高（至少300），并且保存为PNG或PDF格式。避免使用JPG格式，因为它是有损压缩。检查图形尺寸是否足够大（`figsize`参数）。 **问题3：投影扭曲或坐标不对** > 确认你的数据投影与地图投影一致。使用`get_cartopy`获取WRF数据自身的投影信息。如果仍有问题，尝试使用简单的`PlateCarree`投影。 **问题4：内存不足或处理缓慢** > 对于大型WRF输出文件，可以尝试以下优化： > 1. 使用`chunks`参数分块读取：`xr.open_dataset(..., chunks={'Time': 1})` > 2. 只读取需要的变量：`xr.open_dataset(..., drop_variables=['不需要的变量'])` > 3. 在处理前进行空间子集：`ds.sel(south_north=slice(100, 200), west_east=slice(100, 200))` **问题5：颜色条范围不合适** > 手动设置颜色条范围以确保不同图形间的一致性： > ```python > vmin, vmax = 0, 50 # 风速范围0-50 m/s > quiv = ax.quiver(..., clim=(vmin, vmax)) > ``` 掌握了这些技巧后，你就能高效地生成符合学术出版要求的WRF风场分析图。关键是要理解每个步骤背后的原理，而不仅仅是复制代码。随着经验的积累，你会发展出适合自己的工作流程和视觉风格。