从WRF输出文件到专业气象分析Python实战指南气象模拟数据就像一座未经开采的金矿WRF模式输出的wrfout文件包含了海量信息但如何从中提取出有价值的温度、风场和降水数据本文将带你用Python一步步实现从数据提取到专业可视化的完整流程。1. 准备工作与环境配置在开始之前确保你的Python环境已经安装了必要的科学计算库。推荐使用Anaconda创建专用环境conda create -n wrf_analysis python3.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): 大气垂直层次常用气象变量及其提取方法变量名描述单位维度T22米温度KTime, south_north, west_eastU1010米纬向风m/sTime, south_north, west_eastV1010米经向风m/sTime, south_north, west_eastRAINC对流降水mmTime, south_north, west_eastRAINNC非对流降水mmTime, south_north, west_east提取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][:], unitsds[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(projectionccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS) # 绘制温度场 contour ax.contourf(ds[XLONG][0], ds[XLAT][0], t2_celsius[0], levels20, cmapjet, transformccrs.PlateCarree()) # 添加色标和标题 plt.colorbar(contour, labelTemperature (°C)) ax.set_title(2m Temperature at Initial Time) plt.savefig(temperature_map.png, dpi300, bbox_inchestight)风场矢量图fig plt.figure(figsize(12, 8)) ax plt.axes(projectionccrs.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], scale500, colorblack, transformccrs.PlateCarree()) ax.quiverkey(q, X0.9, Y1.05, U10, label10 m/s, labelposE) ax.add_feature(cfeature.COASTLINE) ax.set_title(10m Wind Vectors at Initial Time) plt.savefig(wind_vectors.png, dpi300)降水时间序列- 分析特定点的降水变化# 选择中心点附近的网格 center_x, center_y ds.dims[west_east]//2, ds.dims[south_north]//2 point_rain total_rain.isel(south_northcenter_y, west_eastcenter_x) plt.figure(figsize(10, 5)) plt.plot(times, point_rain, b-, linewidth2) plt.fill_between(times, 0, point_rain, colorb, alpha0.2) plt.title(fAccumulated Precipitation at Grid Point ({center_x}, {center_y})) plt.ylabel(Precipitation (mm)) plt.xticks(rotation45) plt.grid(True) plt.tight_layout() plt.savefig(precipitation_timeseries.png, dpi300)5. 高级分析与实用技巧垂直剖面分析- 研究气象要素的垂直分布# 提取垂直速度变量 w ds[W] # 选择经向剖面(固定东西位置) cross_section w.isel(west_eastcenter_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), levels20, cmapcoolwarm) plt.colorbar(contour, labelVertical 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, dpi300)极端天气识别- 自动检测高温区域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(rotation45) 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, combineby_coords, parallelTrue)内存优化技巧# 只加载需要的变量 ds xr.open_dataset(wrfout_d01_2020-07-01_00:00:00, chunks{Time: 10}, # 使用Dask分块 decode_timesFalse) # 延迟时间解码 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(Timetime_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_diroutput): 处理单个WRF输出文件 os.makedirs(output_dir, exist_okTrue) # 读取数据 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函数或手动插值到同一网格。
从WRF输出文件到实际分析：手把手教你用Python提取并可视化温度、风场和降水数据
从WRF输出文件到专业气象分析Python实战指南气象模拟数据就像一座未经开采的金矿WRF模式输出的wrfout文件包含了海量信息但如何从中提取出有价值的温度、风场和降水数据本文将带你用Python一步步实现从数据提取到专业可视化的完整流程。1. 准备工作与环境配置在开始之前确保你的Python环境已经安装了必要的科学计算库。推荐使用Anaconda创建专用环境conda create -n wrf_analysis python3.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): 大气垂直层次常用气象变量及其提取方法变量名描述单位维度T22米温度KTime, south_north, west_eastU1010米纬向风m/sTime, south_north, west_eastV1010米经向风m/sTime, south_north, west_eastRAINC对流降水mmTime, south_north, west_eastRAINNC非对流降水mmTime, south_north, west_east提取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][:], unitsds[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(projectionccrs.PlateCarree()) # 添加地理特征 ax.add_feature(cfeature.COASTLINE) ax.add_feature(cfeature.BORDERS, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS) # 绘制温度场 contour ax.contourf(ds[XLONG][0], ds[XLAT][0], t2_celsius[0], levels20, cmapjet, transformccrs.PlateCarree()) # 添加色标和标题 plt.colorbar(contour, labelTemperature (°C)) ax.set_title(2m Temperature at Initial Time) plt.savefig(temperature_map.png, dpi300, bbox_inchestight)风场矢量图fig plt.figure(figsize(12, 8)) ax plt.axes(projectionccrs.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], scale500, colorblack, transformccrs.PlateCarree()) ax.quiverkey(q, X0.9, Y1.05, U10, label10 m/s, labelposE) ax.add_feature(cfeature.COASTLINE) ax.set_title(10m Wind Vectors at Initial Time) plt.savefig(wind_vectors.png, dpi300)降水时间序列- 分析特定点的降水变化# 选择中心点附近的网格 center_x, center_y ds.dims[west_east]//2, ds.dims[south_north]//2 point_rain total_rain.isel(south_northcenter_y, west_eastcenter_x) plt.figure(figsize(10, 5)) plt.plot(times, point_rain, b-, linewidth2) plt.fill_between(times, 0, point_rain, colorb, alpha0.2) plt.title(fAccumulated Precipitation at Grid Point ({center_x}, {center_y})) plt.ylabel(Precipitation (mm)) plt.xticks(rotation45) plt.grid(True) plt.tight_layout() plt.savefig(precipitation_timeseries.png, dpi300)5. 高级分析与实用技巧垂直剖面分析- 研究气象要素的垂直分布# 提取垂直速度变量 w ds[W] # 选择经向剖面(固定东西位置) cross_section w.isel(west_eastcenter_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), levels20, cmapcoolwarm) plt.colorbar(contour, labelVertical 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, dpi300)极端天气识别- 自动检测高温区域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(rotation45) 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, combineby_coords, parallelTrue)内存优化技巧# 只加载需要的变量 ds xr.open_dataset(wrfout_d01_2020-07-01_00:00:00, chunks{Time: 10}, # 使用Dask分块 decode_timesFalse) # 延迟时间解码 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(Timetime_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_diroutput): 处理单个WRF输出文件 os.makedirs(output_dir, exist_okTrue) # 读取数据 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函数或手动插值到同一网格。
