从WRF输出文件到实际分析：手把手教你用Python提取并可视化温度、风场和降水数据

从WRF输出文件到专业气象分析：Python实战指南

气象模拟数据就像一座未经开采的金矿，WRF模式输出的wrfout文件包含了海量信息，但如何从中提取出有价值的温度、风场和降水数据？本文将带你用Python一步步实现从数据提取到专业可视化的完整流程。

1. 准备工作与环境配置

在开始之前，确保你的Python环境已经安装了必要的科学计算库。推荐使用Anaconda创建专用环境：

conda create -n wrf_analysis python=3.8 conda activate wrf_analysis conda install -c conda-forge netcdf4 xarray cartopy matplotlib numpy pandas

对于WRF数据处理的特殊需求，还需要安装wrf-python这个专门为WRF数据设计的工具包：

pip install wrf-python

提示：Cartopy是地理空间可视化的重要工具，如果安装遇到问题，可以先安装Proj和GEOS库

检查你的wrfout文件结构，典型的命名格式为wrfout_d01_YYYY-MM-DD_HH:MM:SS，其中d01表示域编号。建议使用以下代码快速查看文件内容：

import xarray as xr ds = xr.open_dataset('wrfout_d01_2020-07-01_00:00:00') print(ds)

2. 理解WRF数据结构和关键变量

WRF数据采用NetCDF格式存储，其多维数组结构需要特别注意维度定义：

• 时间维度(Time): 模拟输出的时间序列
• 南北维度(south_north): 网格的南北方向
• 东西维度(west_east): 网格的东西方向
• 垂直层(bottom_top): 大气垂直层次

常用气象变量及其提取方法：

提取2米温度数据的示例代码：

t2 = ds['T2'] # 获取温度变量 print(f"温度数据形状: {t2.shape}") print(f"温度范围: {t2.min().values:.1f}K 到 {t2.max().values:.1f}K")

3. 数据提取与预处理技巧

直接从WRF输出的原始数据通常需要经过处理才能用于分析。以下是常见的数据处理步骤：

温度单位转换- WRF输出的温度单位为开尔文(K)，转换为摄氏度(℃)：

t2_celsius = t2 - 273.15

风速计算- 从U10和V10分量计算10米风速和风向：

u10 = ds['U10'] v10 = ds['V10'] wind_speed = np.sqrt(u10**2 + v10**2) wind_direction = np.arctan2(-u10, -v10) * (180/np.pi) % 360

降水处理- 计算总降水量和降水强度：

total_rain = ds['RAINC'] + ds['RAINNC'] # 计算降水变化率(需要时间维度上有多个时次) rain_rate = total_rain.diff('Time')

注意：WRF的降水变量是累积值，分析时通常需要计算时段内的增量

时间处理是另一个关键点，WRF使用特殊的时间编码：

from netCDF4 import num2date times = num2date(ds['Times'][:], units=ds['Times'].units) print(f"模拟时段: {times[0]} 到 {times[-1]}")

4. 专业级气象可视化实战

气象数据的可视化需要专业的投影和绘图技术。Cartopy提供了丰富的地图投影功能：

温度场填色图：

import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=ccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle=':') ax.add_feature(cfeature.LAKES, alpha=0.5) ax.add_feature(cfeature.RIVERS) # 绘制温度场 contour = ax.contourf(ds['XLONG'][0], ds['XLAT'][0], t2_celsius[0], levels=20, cmap='jet', transform=ccrs.PlateCarree()) # 添加色标和标题 plt.colorbar(contour, label='Temperature (°C)') ax.set_title('2m Temperature at Initial Time') plt.savefig('temperature_map.png', dpi=300, bbox_inches='tight')

风场矢量图：

fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=ccrs.PlateCarree()) # 每10个点取一个风矢量(避免过于密集) stride = 10 q = ax.quiver(ds['XLONG'][0][::stride,::stride], ds['XLAT'][0][::stride,::stride], u10[0][::stride,::stride], v10[0][::stride,::stride], scale=500, color='black', transform=ccrs.PlateCarree()) ax.quiverkey(q, X=0.9, Y=1.05, U=10, label='10 m/s', labelpos='E') ax.add_feature(cfeature.COASTLINE) ax.set_title('10m Wind Vectors at Initial Time') plt.savefig('wind_vectors.png', dpi=300)

降水时间序列- 分析特定点的降水变化：

# 选择中心点附近的网格 center_x, center_y = ds.dims['west_east']//2, ds.dims['south_north']//2 point_rain = total_rain.isel(south_north=center_y, west_east=center_x) plt.figure(figsize=(10, 5)) plt.plot(times, point_rain, 'b-', linewidth=2) plt.fill_between(times, 0, point_rain, color='b', alpha=0.2) plt.title(f'Accumulated Precipitation at Grid Point ({center_x}, {center_y})') plt.ylabel('Precipitation (mm)') plt.xticks(rotation=45) plt.grid(True) plt.tight_layout() plt.savefig('precipitation_timeseries.png', dpi=300)

5. 高级分析与实用技巧

垂直剖面分析- 研究气象要素的垂直分布：

# 提取垂直速度变量 w = ds['W'] # 选择经向剖面(固定东西位置) cross_section = w.isel(west_east=center_x) fig, ax = plt.subplots(figsize=(12, 6)) contour = ax.contourf(ds['XLAT'][:,center_x,:], ds['ZNU'][0]*ds['P_TOP']/100, # 转换为高度 cross_section.mean('Time'), levels=20, cmap='coolwarm') plt.colorbar(contour, label='Vertical Velocity (m/s)') ax.set_ylabel('Pressure (hPa)') ax.set_xlabel('Latitude') ax.set_title('Meridional Cross-Section of Vertical Velocity') plt.savefig('vertical_profile.png', dpi=300)

极端天气识别- 自动检测高温区域：

heat_wave_threshold = 35 # 35℃ heat_wave_areas = t2_celsius > heat_wave_threshold # 计算每个时次的高温面积占比 heat_wave_fraction = heat_wave_areas.mean(dim=['south_north', 'west_east']) plt.figure(figsize=(10, 5)) plt.plot(times, heat_wave_fraction*100, 'r-') plt.title('Percentage of Domain Experiencing Heat Wave Conditions') plt.ylabel('Area Percentage (%)') plt.xticks(rotation=45) plt.grid(True) plt.tight_layout()

数据导出与共享- 将处理后的数据保存为通用格式：

# 创建新的数据集 processed_ds = xr.Dataset({ 'temperature': t2_celsius, 'wind_speed': wind_speed, 'wind_direction': wind_direction, 'total_precipitation': total_rain }) # 添加坐标信息 processed_ds.coords['longitude'] = ds['XLONG'][0] processed_ds.coords['latitude'] = ds['XLAT'][0] processed_ds.coords['time'] = times # 保存为NetCDF processed_ds.to_netcdf('processed_wrf_data.nc') # 也可以导出为CSV(适合小规模数据) df = processed_ds.to_dataframe() df.to_csv('wrf_data.csv')

6. 性能优化与批量处理

处理多个WRF文件或长时间序列时，性能成为关键考虑因素：

高效读取策略：

# 使用xarray的open_mfdataset处理多个文件 files = sorted(glob.glob('wrfout_d01_*')) ds = xr.open_mfdataset(files, combine='by_coords', parallel=True)

内存优化技巧：

# 只加载需要的变量 ds = xr.open_dataset('wrfout_d01_2020-07-01_00:00:00', chunks={'Time': 10}, # 使用Dask分块 decode_times=False) # 延迟时间解码 selected_vars = ds[['T2', 'U10', 'V10', 'RAINC', 'RAINNC']]

并行处理示例：

from dask.distributed import Client client = Client() # 启动本地集群 # 对多个时次并行计算 results = [] for time_idx in range(len(ds['Time'])): future = client.submit(process_single_time_step, ds.isel(Time=time_idx)) results.append(future) # 获取所有结果 processed_data = client.gather(results)

自动化脚本模板：

import xarray as xr import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from glob import glob import os def process_wrf_file(filename, output_dir='output'): """处理单个WRF输出文件""" os.makedirs(output_dir, exist_ok=True) # 读取数据 ds = xr.open_dataset(filename) # 数据处理 t2_celsius = ds['T2'] - 273.15 u10, v10 = ds['U10'], ds['V10'] wind_speed = np.sqrt(u10**2 + v10**2) total_rain = ds['RAINC'] + ds['RAINNC'] # 可视化 plot_temperature(ds, t2_celsius, output_dir) plot_wind(ds, u10, v10, output_dir) plot_precipitation(ds, total_rain, output_dir) return {'filename': filename, 'max_temp': t2_celsius.max().values} # 批量处理 results = [] for wrf_file in glob('wrfout_d01_*'): results.append(process_wrf_file(wrf_file))

在实际项目中，我发现使用wrf-python库可以显著简化某些操作，特别是涉及到WRF特殊坐标系的计算时。例如，计算相对湿度：

from wrf import getvar # 使用wrf-python提取经过处理的变量 rh = getvar(ds, 'rh2') # 2米相对湿度 slp = getvar(ds, 'slp') # 海平面气压

处理WRF数据时最常见的坑是维度对齐问题，特别是当变量定义在不同的交错网格上时。比如U分量定义在west_east_stag网格上，而V分量定义在south_north_stag网格上，直接运算会导致错误。解决方案是使用wrf-python的interplevel函数或手动插值到同一网格。