简介本资源是一套面向气象与环境科研人员、大气科学研究生及Python气象数据处理初学者的WRF/WRF-Chem全流程自动化脚本工具集聚焦模型前处理地形/土地利用准备、初始场插值、namelist生成与后处理NetCDF解析、时空可视化、特征导出两大核心环节显著提升数值模拟数据链的处理效率与复现性。压缩包含95个文件主体为62个Python脚本涵盖GDAL地理数据处理、Cartopy绘图、PyWRF封装调用、ERA5/SST/VIIRS/Landsat等多源数据接入、15个备份脚本.zbak、5张关键结果示意图如泰勒图、垂直剖面图以及nc实测数据、字体与许可证文件总容量17.17MB。已有37人学习下载资源结构按功能模块分层组织如data_merge、draw_、lcz_pre、ucm_pre等目录附带README.md说明与namelist配置范例提供从原始遥感数据到可输入WRF的二进制地形文件、再到空气质量预测所需气象特征变量的一站式Python解决方案。1. 项目概述为什么我们需要专门的WRF数据处理脚本如果你用过WRF或者WRF-Chem模型肯定对那一大堆NetCDF格式的输出文件又爱又恨。爱的是它包含了从大气环流到化学物种浓度的海量信息恨的是想把这些数据“榨”出点有用的结论过程实在繁琐。官方工具如NCL、ARWpost、RIP4虽然能用但要么学习曲线陡峭要么灵活性不足批量处理多个模拟、提取特定变量、计算自定义诊断量时效率低下。这就是为什么很多研究组和个人开发者包括我自己最终都会转向Python来搭建一套自己的数据处理流水线。这个项目要聊的就是基于Python对WRF/WRF-Chem模型数据进行预处理和后处理的整套脚本方案。所谓“预处理”主要指在模式运行前对输入数据如再分析资料、排放清单进行格式转换、插值、区域裁剪等操作确保它们能乖乖地被WRF系统“吃”进去。而“后处理”则是重头戏指的是模式跑完后对输出的wrfout文件进行解码、计算、分析和可视化把二进制数据变成我们能看懂的图表、统计量和驱动其他模型的数据。为什么用Python因为它生态强大。NetCDF4、xarray、cartopy、wrf-python这些库让读写、处理和可视化模式数据变得前所未有的直观。你可以用几行代码完成过去需要复杂脚本的任务比如计算位涡、追踪气团轨迹、或者批量提取所有模拟中某个站点上空的气溶胶垂直廓线。更重要的是Python脚本的可重复性和可定制性极高一次写好终身受益特别适合需要处理大量情景模拟如不同排放情景、不同物理参数化方案的研究工作。2. 核心工具链选型与环境搭建工欲善其事必先利其器。搭建一个稳定高效的WRF数据处理环境是后续所有工作的基础。下面是我经过多次踩坑后总结出的工具链配置。2.1 Python发行版与包管理选择首先强烈建议使用Miniconda或Anaconda来管理Python环境。WRF相关的科学计算库依赖复杂特别是NetCDF库的底层C依赖conda环境能很好地解决依赖冲突问题。不要使用系统自带的Python以免权限和依赖问题带来麻烦。创建一个专门的环境conda create -n wrf_env python3.9 conda activate wrf_env选择Python 3.9是因为它在稳定性和对新库的兼容性之间取得了很好的平衡。Python 3.10有时会遇到一些科学计算库的预编译包兼容性问题。2.2 核心依赖库详解接下来安装核心库。我推荐使用conda优先安装那些有复杂二进制依赖的包用pip作为补充。# 使用conda安装基础科学计算栈和地理处理库 conda install -c conda-forge numpy scipy pandas matplotlib jupyter conda install -c conda-forge xarray dask netcdf4 h5py conda install -c conda-forge cartopy proj geos shapely conda install -c conda-forge cfgrib eccodes # 如需处理GRIB格式数据重点库解析xarray这是处理WRF数据的核心。它完美封装了NetCDF数据模型提供了类似pandas的标签化数据操作接口。你可以通过维度名如Time,south_north,bottom_top和坐标值来切片数据比直接操作NetCDF4库直观太多。wrf-python由NCAR官方维护的Python工具包是后处理的瑞士军刀。它提供了一系列函数来从WRF输出中提取和计算标准变量如气压、温度、风速、插值到指定高度或气压层、计算诊断量如相对湿度、潜在温度以及基本的可视化功能。必须用pip安装pip install wrf-python注意wrf-python依赖于netCDF4和xarray但通常用pip安装时会自动处理。cartopy专业的地图绘图库。相比Basemap已停止维护Cartopy更现代、功能更强大与matplotlib集成度极高用于绘制带有地图投影的天气图、绘制模式区域、叠加海岸线、国界等。netCDF4提供底层API当需要极高性能或操作xarray未覆盖的底层属性时使用。2.3 集成开发环境IDE配置对于此类涉及大量数据探索和脚本编写的工作一个好的IDE至关重要。Visual Studio Code (VSCode)是目前的最佳选择之一。安装VSCode及Python扩展从官网安装VSCode然后在扩展商店搜索并安装“Python”扩展由Microsoft发布。选择解释器在VSCode中按CtrlShiftP输入“Python: Select Interpreter”选择我们刚才创建的wrf_env环境下的Python路径通常类似~/miniconda3/envs/wrf_env/bin/python。配置Jupyter内核在wrf_env环境中运行python -m ipykernel install --user --namewrf_env这样在VSCode的Jupyter Notebook中就可以选择wrf_env内核了。用Notebook做数据探索和原型开发非常方便。实用扩展推荐Rainbow CSV高亮显示CSV格式的排放清单文件。Even Better TOML如果你需要编辑WRF的namelist.input文件虽然不是TOML但语法高亮有助阅读。GitLens方便进行版本控制。注意在服务器或无图形界面的环境中熟练使用vim/nano编辑脚本配合tmux进行会话管理是必备技能。可以将常用处理流程封装成命令行脚本通过参数控制。3. WRF-Chem数据预处理实战以人为源排放清单为例WRF-Chem的预处理核心之一是准备化学初始条件和边界条件以及最重要的——排放清单。这里我们以处理全球排放清单如CEDS、EDGAR并将其转换为WRF-Chem所需的MOZART、RADM2等化学机制格式为例详解流程。3.1 排放清单预处理流程总览原始全球排放清单通常是年度、月度或每日的NetCDF或文本文件网格分辨率较低如0.5°×0.5°。我们需要将其时间插值从月均值或日均值插值到模拟所需的每小时数据。空间插值与重映射从源网格插值到WRF模拟域往往是更高分辨率的兰伯特投影或墨卡托投影。化学物种映射将清单中的排放物种如NOx, SO2, NMVOC映射到目标化学机制所定义的物种。垂直分配将地面排放量按预设的垂直剖面分配到不同的模型层。格式转换最终生成WRF-Chem预处理程序anthro_emis可读取的NetCDF格式。3.2 使用Python进行时空插值与重映射假设我们有一个CEDS的NOx月度排放文件NOx_emissions_2015_monthly_0.5x0.5.nc。import xarray as xr import numpy as np from datetime import datetime, timedelta import scipy.interpolate as spi # 1. 读取原始排放数据 ds_orig xr.open_dataset(NOx_emissions_2015_monthly_0.5x0.5.nc) # 假设变量名为 NOx维度为 (time, lat, lon) nox_orig ds_orig[NOx] # 2. 时间插值从月度到每日 # 创建原始数据的时间坐标假设是每月15日 times_orig nox_orig.time.values # 类型为 numpy.datetime64 # 创建目标时间序列例如2015年全年每日 start_date datetime(2015, 1, 1) end_date datetime(2015, 12, 31) target_dates [start_date timedelta(daysi) for i in range((end_date - start_date).days 1)] target_dates_np np.array(target_dates, dtypedatetime64[D]) # 使用线性插值对于排放更复杂的方法可能考虑周末/工作日因子 # 这里需要将时间转换为数值进行插值简单演示思路 # 实际应用中建议对每个网格点的时间序列进行插值 nox_daily nox_orig.interp(timetarget_dates_np, methodlinear) # 3. 空间重映射从经纬度网格到WRF网格 # 首先需要你的WRF域定义文件通常是geo_em.d01.nc中的经纬度变量 ds_wrf xr.open_dataset(geo_em.d01.nc) wrf_lat ds_wrf[XLAT_M].isel(Time0).values # 网格中心纬度 wrf_lon ds_wrf[XLONG_M].isel(Time0).values # 网格中心经度 # 使用scipy进行双线性插值 def regrid_to_wrf(orig_data, orig_lat, orig_lon, target_lat, target_lon): 将数据从原始网格插值到WRF网格 # orig_data: 2D数组 (lat, lon) # 创建插值函数 interp_func spi.RegularGridInterpolator((orig_lat, orig_lon), orig_data, methodlinear, bounds_errorFalse, fill_value0.0) # 准备目标网格点 target_points np.stack([target_lat.ravel(), target_lon.ravel()], axis-1) # 插值 regridded interp_func(target_points).reshape(target_lat.shape) return regridded # 假设我们处理某一天的排放数据 sample_day nox_daily.isel(time0) orig_lat sample_day.lat.values orig_lon sample_day.lon.values orig_vals sample_day.values nox_regridded regrid_to_wrf(orig_vals, orig_lat, orig_lon, wrf_lat, wrf_lon) # 将结果保存为新的xarray DataArray并添加WRF网格坐标 nox_wrf_da xr.DataArray(nox_regridded, dims[south_north, west_east], coords{XLAT: ([south_north, west_east], wrf_lat), XLONG: ([south_north, west_east], wrf_lon)})实操心得空间插值是个性能瓶颈。对于全球高分辨率数据上述循环方法会很慢。生产环境中推荐使用xESMF库它专门用于地球科学数据的重网格化支持保守插值等方法并且可以利用Dask进行并行计算效率极高。3.3 化学物种映射与垂直分配不同的化学机制对物种的定义不同。例如NMVOC非甲烷挥发性有机物在MOZART机制下被拆分成数十个具体物种。你需要一个映射文件通常是CSV格式定义源清单物种到目标机制物种的分配系数。import pandas as pd # 读取物种映射表 species_map_df pd.read_csv(MOZART_species_mapping.csv) # 假设列名为Source_Species, Target_Species, Split_Factor # 假设我们已处理完所有源物种的二维排放如NOx, SO2, CO, NMVOC等 # 并存储在字典中键为物种名值为对应的二维DataArray emissions_2d {NOx: nox_wrf_da, SO2: so2_wrf_da, ...} # 初始化一个字典用于存放目标机制的三维排放考虑了垂直分配 # 假设我们有垂直分配系数文件定义了不同物种/不同排放类型如工业、交通的垂直剖面共bottom_top层 vert_profile np.loadtxt(vertical_profile.csv, delimiter,) # 形状 (n_species, bottom_top) target_emissions_3d {} for idx, row in species_map_df.iterrows(): src_species row[Source_Species] tgt_species row[Target_Species] factor row[Split_Factor] vert_factor vert_profile[tgt_species] # 这里需要根据tgt_species索引仅为示例 if src_species in emissions_2d: # 水平排放 * 物种分配系数 emis_2d_scaled emissions_2d[src_species] * factor # 扩展到垂直维度 # 将二维排放 (south_north, west_east) 乘以垂直剖面 (bottom_top,) # 通过增加新维度并广播实现 emis_3d emis_2d_scaled.expand_dims(dim{bottom_top: len(vert_factor)}, axis0) * vert_factor[:, np.newaxis, np.newaxis] # 注意调整维度顺序以匹配WRF (Time, bottom_top, south_north, west_east) emis_3d emis_3d.transpose(bottom_top, south_north, west_east) if tgt_species not in target_emissions_3d: target_emissions_3d[tgt_species] emis_3d else: target_emissions_3d[tgt_species] emis_3d最后将target_emissions_3d字典中的所有DataArray按照WRF-Chem排放文件要求的格式变量名、属性、时间维度组装成一个xarray Dataset并写入NetCDF文件。这个文件就可以作为anthro_emis或fire_emis等预处理程序的输入。4. WRF模型后处理核心从wrfout到分析就绪数据模式跑完了拿到了一串wrfout_d01_*文件真正的分析才刚刚开始。后处理的目标是将这些文件转换为易于分析的数据集。4.1 使用wrf-python提取标准气象变量wrf-python是后处理的首选工具它封装了WRF内部计算诊断量的Fortran例程确保结果与WRF内部计算一致。import wrf import xarray as xr import numpy as np # 打开多个时间序列的wrfout文件 wrf_files [wrfout_d01_2023-07-01_00:00:00, wrfout_d01_2023-07-01_01:00:00] ds_wrf xr.open_mfdataset(wrf_files, combineby_coords, enginenetcdf4) # 提取海平面气压、10米风场、2米温湿度等近地面变量 # 这些函数直接返回xarray DataArray非常方便 slp wrf.getvar(ds_wrf, slp) # 海平面气压 u10, v10 wrf.getvar(ds_wrf, uvmet10) # 10米风地图投影坐标 t2 wrf.getvar(ds_wrf, T2) # 2米温度 rh2 wrf.getvar(ds_wrf, rh2) # 2米相对湿度 # 提取等压面变量将数据插值到指定的气压层如500hPa, 850hPa pressure_levels [1000, 850, 700, 500, 300] # 获取三维气压场 p wrf.getvar(ds_wrf, pressure) # 获取三维温度、位势高度、风场 tk wrf.getvar(ds_wrf, tk) # 温度 (K) z wrf.getvar(ds_wrf, z) # 位势高度 (m) ua, va wrf.getvar(ds_wrf, uvmet) # 水平风 (地图投影坐标) # 插值到等压面 tk_plev wrf.interplevel(tk, p, pressure_levels) z_plev wrf.interplevel(z, p, pressure_levels) ua_plev wrf.interplevel(ua, p, pressure_levels) va_plev wrf.interplevel(va, p, pressure_levels)4.2 计算自定义诊断量以潜在温度和涡度为例除了标准变量研究中经常需要计算一些诊断量。# 计算潜在温度 (Theta) # 方法1使用wrf-python内置函数推荐 theta wrf.getvar(ds_wrf, theta) # 方法2手动计算便于理解原理 # 潜在温度公式: Theta T * (P0 / P)^(R/cp) # 其中 P0 1000 hPa, R287 J/kg/K, cp1004 J/kg/K P wrf.getvar(ds_wrf, pressure) # 单位: hPa T wrf.getvar(ds_wrf, tk) # 单位: K P0 1000.0 R 287.0 cp 1004.0 theta_manual T * (P0 / P) ** (R / cp) # 计算相对涡度 (Relative Vorticity) # 首先获取地图投影下的风分量和科里奥利参数 u wrf.getvar(ds_wrf, ua) # 质量点上的东向风分量 v wrf.getvar(ds_wrf, va) # 质量点上的北向风分量 # 计算涡度 (使用wrf-python的梯度函数它考虑了地图投影变形) # 这比手动计算dx, dy要准确得多 vort wrf.vorticity(u, v)4.3 高效处理大型数据集与时间序列分析WRF输出文件动辄几十GB一次性读入内存不现实。xarray与dask的集成提供了完美的解决方案。# 使用open_mfdataset并启用dask延迟计算 # 指定chunks参数将数据在逻辑上分块不立即加载 ds_big xr.open_mfdataset(wrfout_d01_*.nc, combineby_coords, parallelTrue, chunks{Time: 10, south_north: 100, west_east: 100}) print(ds_big) # 此时数据并未真正加载只是创建了计算任务图。 # 进行一个耗时的计算例如计算整个时间序列的区域平均温度 # 这个计算会被dask分解到各个数据块上并行执行 mean_t2_series ds_big[T2].mean(dim[south_north, west_east]) # 现在才触发实际计算 mean_t2_series_computed mean_t2_series.compute() # 将结果保存到小文件中避免重复计算 mean_t2_series_computed.to_netcdf(domain_mean_T2.nc) # 对于更复杂的操作可以自定义函数并应用 def calculate_aqi(pm25_conc): 一个简化的PM2.5 AQI计算函数示例 # 根据浓度分段线性计算AQI # ... 具体计算逻辑 ... return aqi # 假设有PM2.5浓度变量 pm25 wrf.getvar(ds_big, PM2_5_DRY) # 注意变量名随化学机制而异 # 使用apply_ufunc进行逐元素计算并利用dask并行 aqi xr.apply_ufunc(calculate_aqi, pm25, daskparallelized, output_dtypes[np.float32]) aqi_computed aqi.compute()注意事项使用dask时chunks的大小设置是关键。太小会导致任务调度开销巨大太大则可能内存不足。一个经验法则是每个数据块的大小应在10MB到100MB之间。可以通过ds_big.nbytes / 1e6查看数据总大小再除以设定的块数来估算。5. 可视化与成果输出从数据到图表数据分析的最终目的是为了理解和展示。基于Cartopy和Matplotlib我们可以制作出版级别的图表。5.1 绘制二维场与叠加地图import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from wrf import get_cartopy, latlon_coords # 选择一个时间点 slp_plot slp.isel(Time0) # 获取该变量的经纬度坐标和地图投影信息 lats, lons latlon_coords(slp_plot) cart_proj get_cartopy(slp_plot) # 创建图形和坐标系 fig plt.figure(figsize(12, 8)) ax plt.axes(projectioncart_proj) # 设置地图范围通常就是模式域的范围 ax.set_xlim(wrf.cartopy_xlim(slp_plot)) ax.set_ylim(wrf.cartopy_ylim(slp_plot)) # 添加地理特征 ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5, linestyle:) ax.add_feature(cfeature.STATES, linewidth0.3, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS, linewidth0.5) # 绘制海平面气压填色图 contourf ax.contourf(lons, lats, slp_plot, levels20, transformccrs.PlateCarree(), cmapRdBu_r) # 叠加等值线 contour ax.contour(lons, lats, slp_plot, levels10, colorsk, linewidths0.5, transformccrs.PlateCarree()) ax.clabel(contour, inlineTrue, fontsize8, fmt%1.0f) # 绘制风矢以10米风为例 u10_plot u10.isel(Time0) v10_plot v10.isel(Time0) # 每隔N个点取一个风矢避免过于密集 stride 10 ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride], u10_plot.values[::stride, ::stride], v10_plot.values[::stride, ::stride], transformccrs.PlateCarree(), scale300, colorgreen) # 添加色标、标题 plt.colorbar(contourf, axax, orientationhorizontal, pad0.05, aspect40, labelSea Level Pressure (hPa)) ax.set_title(fSea Level Pressure and 10m Wind\n{slp_plot.Time.dt.strftime(%Y-%m-%d %H:%M UTC).values}) plt.tight_layout() plt.savefig(slp_wind_map.png, dpi300, bbox_inchestight) plt.show()5.2 绘制垂直剖面与时空序列图垂直剖面对于分析天气系统或污染物的垂直结构至关重要。# 假设我们想画一条从点A(lat1, lon1)到点B(lat2, lon2)的垂直剖面 cross_start (40.0, 115.0) # (lat, lon) 起点 cross_end (35.0, 120.0) # (lat, lon) 终点 # 使用wrf-python的剖面函数 z_cross wrf.vertcross(z, p, wrf_start_pointcross_start, wrf_end_pointcross_end, latlonTrue, metaTrue) tk_cross wrf.vertcross(tk, p, wrf_start_pointcross_start, wrf_end_pointcross_end, latlonTrue, metaTrue) # 获取剖面线的水平距离坐标 cross_lats wrf.to_np(z_cross.coords[xy_loc]).T[1] # 纬度 cross_lons wrf.to_np(z_cross.coords[xy_loc]).T[0] # 经度 cross_dist wrf.to_np(z_cross.coords[xy_loc]).T[2] # 距离起点距离 (km) fig, ax plt.subplots(figsize(15, 5)) # 绘制温度剖面填色图 pc ax.pcolormesh(cross_dist, wrf.to_np(z_cross.coords[vertical]), wrf.to_np(tk_cross), cmapSpectral_r, shadingauto) # 叠加位势高度等值线 contour ax.contour(cross_dist, wrf.to_np(z_cross.coords[vertical]), wrf.to_np(z_cross), colorsk, linewidths0.8) ax.clabel(contour, inlineTrue, fontsize8, fmt%1.0f) ax.set_xlabel(Distance from Start (km)) ax.set_ylabel(Height (m)) ax.set_title(fTemperature Cross Section\nStart {cross_start} - End {cross_end}) plt.colorbar(pc, axax, labelTemperature (K)) # 在x轴下方标记经纬度 ax_btm ax.twiny() ax_btm.set_xlim(ax.get_xlim()) xticks ax.get_xticks() xtick_lats np.interp(xticks, cross_dist, cross_lats) xtick_lons np.interp(xticks, cross_dist, cross_lons) xtick_labels [f({lat:.1f}N,\n{lon:.1f}E) for lat, lon in zip(xtick_lats, xtick_lons)] ax_btm.set_xticks(xticks) ax_btm.set_xticklabels(xtick_labels, fontsize8) ax_btm.set_xlabel(Latitude, Longitude) plt.tight_layout() plt.savefig(vertical_cross_section.png, dpi300)6. 实战问题排查与脚本优化经验在实际操作中你会遇到各种各样的问题。下面是一些常见坑点和解决技巧。6.1 常见错误与解决方法问题现象可能原因解决方案ImportError: libnetcdf.so.xx: cannot open shared object fileNetCDF库的底层C依赖未正确安装或路径不对。使用conda install -c conda-forge netcdf4重新安装。确保conda环境已激活。检查LD_LIBRARY_PATH环境变量。wrf.getvar返回None或报错变量名错误或该变量在当前WRF输出中不存在。使用wrf.variables()查看文件中所有变量。WRF-Chem的化学变量名取决于机制如PM2_5_DRYMADE/SORGAM、o3等。插值到等压面时出现大量NaN值请求的气压层超出了模式层的气压范围如地面气压只有980hPa却请求了1000hPa层。检查模式域的最低层气压。使用wrf.getvar(ds, ‘pressure’)查看气压范围。插值时设置missing np.nan并处理NaN。Cartopy绘图时地图变形或错位地图投影cart_proj与数据投影不匹配。务必使用wrf.get_cartopy(var)获取数据自带的地图投影对象而不是自己定义。处理大量文件时内存爆炸没有使用分块chunk或延迟加载。使用xr.open_mfdataset(..., chunks{})。对于聚合操作如mean,sum先分块计算再合并结果。避免在循环中重复打开文件。时间坐标解析错误WRF输出文件的时间属性格式不标准。使用wrf.extract_times(ds, timeidx0)来获取正确的时间信息。或者用xr.decode_cfTrue打开数据集。化学物种浓度单位混乱WRF-Chem输出单位可能是mol mol-1 dry air,ppm,μg m-3等。仔细查阅你所使用的化学机制文档如MOZART, RADM2, SAPRC99的输出说明。使用wrf.getvar提取的变量通常会自动转换为标准单位如K, hPa, m/s。6.2 脚本性能优化技巧向量化操作永远避免在Python中对大型数组使用显式循环for循环。尽量使用NumPy或xarray的向量化函数和ufunc。利用Dask进行并行计算对于多文件、多时间步的处理open_mfdataset配合chunks和.compute()是最佳实践。可以设置parallelTrue加速文件读取。选择性读取如果只需要少数几个变量或特定时间步使用xr.open_dataset(..., drop_variables[...])或preprocess函数在打开时即丢弃不需要的数据。缓存中间结果对于耗时的计算步骤如重网格化、复杂的诊断量计算将结果保存为NetCDF文件。下次直接读取中间文件避免重复计算。使用更高效的文件格式如果频繁读写可以考虑将处理后的数据保存为zarr格式它对分块存储和并行读写更友好。6.3 构建可复用的处理流水线为了提升效率建议将通用功能模块化# my_wrf_tools.py import xarray as xr import wrf import numpy as np class WRFProcessor: def __init__(self, wrfout_paths): self.ds xr.open_mfdataset(wrfout_paths, combineby_coords, chunks{Time: 24}) self.cache {} # 用于缓存计算昂贵的变量 def get_diagnostic_var(self, var_name): 获取诊断变量如果缓存中有则直接返回 if var_name in self.cache: return self.cache[var_name] if var_name pv: # 计算位涡 p wrf.getvar(self.ds, pressure) theta wrf.getvar(self.ds, theta) ua, va, wa wrf.getvar(self.ds, uvmet) pv wrf.potential_vorticity(p, theta, ua, va, wa) self.cache[pv] pv return pv elif var_name precip_total: # 计算总降水量 rainc self.ds[RAINC] rainnc self.ds[RAINNC] total rainc rainnc self.cache[precip_total] total return total else: # 尝试用wrf.getvar直接提取 var wrf.getvar(self.ds, var_name) self.cache[var_name] var return var def save_to_disk(self, var_list, output_file): 将指定变量列表保存到新的NetCDF文件 ds_out xr.Dataset() for var_name in var_list: ds_out[var_name] self.get_diagnostic_var(var_name) ds_out.to_netcdf(output_file) print(fSaved to {output_file}) # 主脚本 main.py from my_wrf_tools import WRFProcessor processor WRFProcessor([wrfout_d01_*]) # 一次性计算并保存多个常用诊断量 processor.save_to_disk([slp, T2, rh2, uvmet10, precip_total], processed_diagnostics.nc)这套脚本体系的核心思想是自动化和可重复。一旦搭建完成对于新的模拟任务你只需要修改输入文件路径和少数参数就能快速得到分析就绪的数据和图表将精力从繁琐的数据处理中解放出来更多地投入到科学问题的探索上。本文还有配套的精品资源点击获取 基于Python的WRF/WRF-Chem数据处理：从环境搭建到实战应用
简介本资源是一套面向气象与环境科研人员、大气科学研究生及Python气象数据处理初学者的WRF/WRF-Chem全流程自动化脚本工具集聚焦模型前处理地形/土地利用准备、初始场插值、namelist生成与后处理NetCDF解析、时空可视化、特征导出两大核心环节显著提升数值模拟数据链的处理效率与复现性。压缩包含95个文件主体为62个Python脚本涵盖GDAL地理数据处理、Cartopy绘图、PyWRF封装调用、ERA5/SST/VIIRS/Landsat等多源数据接入、15个备份脚本.zbak、5张关键结果示意图如泰勒图、垂直剖面图以及nc实测数据、字体与许可证文件总容量17.17MB。已有37人学习下载资源结构按功能模块分层组织如data_merge、draw_、lcz_pre、ucm_pre等目录附带README.md说明与namelist配置范例提供从原始遥感数据到可输入WRF的二进制地形文件、再到空气质量预测所需气象特征变量的一站式Python解决方案。1. 项目概述为什么我们需要专门的WRF数据处理脚本如果你用过WRF或者WRF-Chem模型肯定对那一大堆NetCDF格式的输出文件又爱又恨。爱的是它包含了从大气环流到化学物种浓度的海量信息恨的是想把这些数据“榨”出点有用的结论过程实在繁琐。官方工具如NCL、ARWpost、RIP4虽然能用但要么学习曲线陡峭要么灵活性不足批量处理多个模拟、提取特定变量、计算自定义诊断量时效率低下。这就是为什么很多研究组和个人开发者包括我自己最终都会转向Python来搭建一套自己的数据处理流水线。这个项目要聊的就是基于Python对WRF/WRF-Chem模型数据进行预处理和后处理的整套脚本方案。所谓“预处理”主要指在模式运行前对输入数据如再分析资料、排放清单进行格式转换、插值、区域裁剪等操作确保它们能乖乖地被WRF系统“吃”进去。而“后处理”则是重头戏指的是模式跑完后对输出的wrfout文件进行解码、计算、分析和可视化把二进制数据变成我们能看懂的图表、统计量和驱动其他模型的数据。为什么用Python因为它生态强大。NetCDF4、xarray、cartopy、wrf-python这些库让读写、处理和可视化模式数据变得前所未有的直观。你可以用几行代码完成过去需要复杂脚本的任务比如计算位涡、追踪气团轨迹、或者批量提取所有模拟中某个站点上空的气溶胶垂直廓线。更重要的是Python脚本的可重复性和可定制性极高一次写好终身受益特别适合需要处理大量情景模拟如不同排放情景、不同物理参数化方案的研究工作。2. 核心工具链选型与环境搭建工欲善其事必先利其器。搭建一个稳定高效的WRF数据处理环境是后续所有工作的基础。下面是我经过多次踩坑后总结出的工具链配置。2.1 Python发行版与包管理选择首先强烈建议使用Miniconda或Anaconda来管理Python环境。WRF相关的科学计算库依赖复杂特别是NetCDF库的底层C依赖conda环境能很好地解决依赖冲突问题。不要使用系统自带的Python以免权限和依赖问题带来麻烦。创建一个专门的环境conda create -n wrf_env python3.9 conda activate wrf_env选择Python 3.9是因为它在稳定性和对新库的兼容性之间取得了很好的平衡。Python 3.10有时会遇到一些科学计算库的预编译包兼容性问题。2.2 核心依赖库详解接下来安装核心库。我推荐使用conda优先安装那些有复杂二进制依赖的包用pip作为补充。# 使用conda安装基础科学计算栈和地理处理库 conda install -c conda-forge numpy scipy pandas matplotlib jupyter conda install -c conda-forge xarray dask netcdf4 h5py conda install -c conda-forge cartopy proj geos shapely conda install -c conda-forge cfgrib eccodes # 如需处理GRIB格式数据重点库解析xarray这是处理WRF数据的核心。它完美封装了NetCDF数据模型提供了类似pandas的标签化数据操作接口。你可以通过维度名如Time,south_north,bottom_top和坐标值来切片数据比直接操作NetCDF4库直观太多。wrf-python由NCAR官方维护的Python工具包是后处理的瑞士军刀。它提供了一系列函数来从WRF输出中提取和计算标准变量如气压、温度、风速、插值到指定高度或气压层、计算诊断量如相对湿度、潜在温度以及基本的可视化功能。必须用pip安装pip install wrf-python注意wrf-python依赖于netCDF4和xarray但通常用pip安装时会自动处理。cartopy专业的地图绘图库。相比Basemap已停止维护Cartopy更现代、功能更强大与matplotlib集成度极高用于绘制带有地图投影的天气图、绘制模式区域、叠加海岸线、国界等。netCDF4提供底层API当需要极高性能或操作xarray未覆盖的底层属性时使用。2.3 集成开发环境IDE配置对于此类涉及大量数据探索和脚本编写的工作一个好的IDE至关重要。Visual Studio Code (VSCode)是目前的最佳选择之一。安装VSCode及Python扩展从官网安装VSCode然后在扩展商店搜索并安装“Python”扩展由Microsoft发布。选择解释器在VSCode中按CtrlShiftP输入“Python: Select Interpreter”选择我们刚才创建的wrf_env环境下的Python路径通常类似~/miniconda3/envs/wrf_env/bin/python。配置Jupyter内核在wrf_env环境中运行python -m ipykernel install --user --namewrf_env这样在VSCode的Jupyter Notebook中就可以选择wrf_env内核了。用Notebook做数据探索和原型开发非常方便。实用扩展推荐Rainbow CSV高亮显示CSV格式的排放清单文件。Even Better TOML如果你需要编辑WRF的namelist.input文件虽然不是TOML但语法高亮有助阅读。GitLens方便进行版本控制。注意在服务器或无图形界面的环境中熟练使用vim/nano编辑脚本配合tmux进行会话管理是必备技能。可以将常用处理流程封装成命令行脚本通过参数控制。3. WRF-Chem数据预处理实战以人为源排放清单为例WRF-Chem的预处理核心之一是准备化学初始条件和边界条件以及最重要的——排放清单。这里我们以处理全球排放清单如CEDS、EDGAR并将其转换为WRF-Chem所需的MOZART、RADM2等化学机制格式为例详解流程。3.1 排放清单预处理流程总览原始全球排放清单通常是年度、月度或每日的NetCDF或文本文件网格分辨率较低如0.5°×0.5°。我们需要将其时间插值从月均值或日均值插值到模拟所需的每小时数据。空间插值与重映射从源网格插值到WRF模拟域往往是更高分辨率的兰伯特投影或墨卡托投影。化学物种映射将清单中的排放物种如NOx, SO2, NMVOC映射到目标化学机制所定义的物种。垂直分配将地面排放量按预设的垂直剖面分配到不同的模型层。格式转换最终生成WRF-Chem预处理程序anthro_emis可读取的NetCDF格式。3.2 使用Python进行时空插值与重映射假设我们有一个CEDS的NOx月度排放文件NOx_emissions_2015_monthly_0.5x0.5.nc。import xarray as xr import numpy as np from datetime import datetime, timedelta import scipy.interpolate as spi # 1. 读取原始排放数据 ds_orig xr.open_dataset(NOx_emissions_2015_monthly_0.5x0.5.nc) # 假设变量名为 NOx维度为 (time, lat, lon) nox_orig ds_orig[NOx] # 2. 时间插值从月度到每日 # 创建原始数据的时间坐标假设是每月15日 times_orig nox_orig.time.values # 类型为 numpy.datetime64 # 创建目标时间序列例如2015年全年每日 start_date datetime(2015, 1, 1) end_date datetime(2015, 12, 31) target_dates [start_date timedelta(daysi) for i in range((end_date - start_date).days 1)] target_dates_np np.array(target_dates, dtypedatetime64[D]) # 使用线性插值对于排放更复杂的方法可能考虑周末/工作日因子 # 这里需要将时间转换为数值进行插值简单演示思路 # 实际应用中建议对每个网格点的时间序列进行插值 nox_daily nox_orig.interp(timetarget_dates_np, methodlinear) # 3. 空间重映射从经纬度网格到WRF网格 # 首先需要你的WRF域定义文件通常是geo_em.d01.nc中的经纬度变量 ds_wrf xr.open_dataset(geo_em.d01.nc) wrf_lat ds_wrf[XLAT_M].isel(Time0).values # 网格中心纬度 wrf_lon ds_wrf[XLONG_M].isel(Time0).values # 网格中心经度 # 使用scipy进行双线性插值 def regrid_to_wrf(orig_data, orig_lat, orig_lon, target_lat, target_lon): 将数据从原始网格插值到WRF网格 # orig_data: 2D数组 (lat, lon) # 创建插值函数 interp_func spi.RegularGridInterpolator((orig_lat, orig_lon), orig_data, methodlinear, bounds_errorFalse, fill_value0.0) # 准备目标网格点 target_points np.stack([target_lat.ravel(), target_lon.ravel()], axis-1) # 插值 regridded interp_func(target_points).reshape(target_lat.shape) return regridded # 假设我们处理某一天的排放数据 sample_day nox_daily.isel(time0) orig_lat sample_day.lat.values orig_lon sample_day.lon.values orig_vals sample_day.values nox_regridded regrid_to_wrf(orig_vals, orig_lat, orig_lon, wrf_lat, wrf_lon) # 将结果保存为新的xarray DataArray并添加WRF网格坐标 nox_wrf_da xr.DataArray(nox_regridded, dims[south_north, west_east], coords{XLAT: ([south_north, west_east], wrf_lat), XLONG: ([south_north, west_east], wrf_lon)})实操心得空间插值是个性能瓶颈。对于全球高分辨率数据上述循环方法会很慢。生产环境中推荐使用xESMF库它专门用于地球科学数据的重网格化支持保守插值等方法并且可以利用Dask进行并行计算效率极高。3.3 化学物种映射与垂直分配不同的化学机制对物种的定义不同。例如NMVOC非甲烷挥发性有机物在MOZART机制下被拆分成数十个具体物种。你需要一个映射文件通常是CSV格式定义源清单物种到目标机制物种的分配系数。import pandas as pd # 读取物种映射表 species_map_df pd.read_csv(MOZART_species_mapping.csv) # 假设列名为Source_Species, Target_Species, Split_Factor # 假设我们已处理完所有源物种的二维排放如NOx, SO2, CO, NMVOC等 # 并存储在字典中键为物种名值为对应的二维DataArray emissions_2d {NOx: nox_wrf_da, SO2: so2_wrf_da, ...} # 初始化一个字典用于存放目标机制的三维排放考虑了垂直分配 # 假设我们有垂直分配系数文件定义了不同物种/不同排放类型如工业、交通的垂直剖面共bottom_top层 vert_profile np.loadtxt(vertical_profile.csv, delimiter,) # 形状 (n_species, bottom_top) target_emissions_3d {} for idx, row in species_map_df.iterrows(): src_species row[Source_Species] tgt_species row[Target_Species] factor row[Split_Factor] vert_factor vert_profile[tgt_species] # 这里需要根据tgt_species索引仅为示例 if src_species in emissions_2d: # 水平排放 * 物种分配系数 emis_2d_scaled emissions_2d[src_species] * factor # 扩展到垂直维度 # 将二维排放 (south_north, west_east) 乘以垂直剖面 (bottom_top,) # 通过增加新维度并广播实现 emis_3d emis_2d_scaled.expand_dims(dim{bottom_top: len(vert_factor)}, axis0) * vert_factor[:, np.newaxis, np.newaxis] # 注意调整维度顺序以匹配WRF (Time, bottom_top, south_north, west_east) emis_3d emis_3d.transpose(bottom_top, south_north, west_east) if tgt_species not in target_emissions_3d: target_emissions_3d[tgt_species] emis_3d else: target_emissions_3d[tgt_species] emis_3d最后将target_emissions_3d字典中的所有DataArray按照WRF-Chem排放文件要求的格式变量名、属性、时间维度组装成一个xarray Dataset并写入NetCDF文件。这个文件就可以作为anthro_emis或fire_emis等预处理程序的输入。4. WRF模型后处理核心从wrfout到分析就绪数据模式跑完了拿到了一串wrfout_d01_*文件真正的分析才刚刚开始。后处理的目标是将这些文件转换为易于分析的数据集。4.1 使用wrf-python提取标准气象变量wrf-python是后处理的首选工具它封装了WRF内部计算诊断量的Fortran例程确保结果与WRF内部计算一致。import wrf import xarray as xr import numpy as np # 打开多个时间序列的wrfout文件 wrf_files [wrfout_d01_2023-07-01_00:00:00, wrfout_d01_2023-07-01_01:00:00] ds_wrf xr.open_mfdataset(wrf_files, combineby_coords, enginenetcdf4) # 提取海平面气压、10米风场、2米温湿度等近地面变量 # 这些函数直接返回xarray DataArray非常方便 slp wrf.getvar(ds_wrf, slp) # 海平面气压 u10, v10 wrf.getvar(ds_wrf, uvmet10) # 10米风地图投影坐标 t2 wrf.getvar(ds_wrf, T2) # 2米温度 rh2 wrf.getvar(ds_wrf, rh2) # 2米相对湿度 # 提取等压面变量将数据插值到指定的气压层如500hPa, 850hPa pressure_levels [1000, 850, 700, 500, 300] # 获取三维气压场 p wrf.getvar(ds_wrf, pressure) # 获取三维温度、位势高度、风场 tk wrf.getvar(ds_wrf, tk) # 温度 (K) z wrf.getvar(ds_wrf, z) # 位势高度 (m) ua, va wrf.getvar(ds_wrf, uvmet) # 水平风 (地图投影坐标) # 插值到等压面 tk_plev wrf.interplevel(tk, p, pressure_levels) z_plev wrf.interplevel(z, p, pressure_levels) ua_plev wrf.interplevel(ua, p, pressure_levels) va_plev wrf.interplevel(va, p, pressure_levels)4.2 计算自定义诊断量以潜在温度和涡度为例除了标准变量研究中经常需要计算一些诊断量。# 计算潜在温度 (Theta) # 方法1使用wrf-python内置函数推荐 theta wrf.getvar(ds_wrf, theta) # 方法2手动计算便于理解原理 # 潜在温度公式: Theta T * (P0 / P)^(R/cp) # 其中 P0 1000 hPa, R287 J/kg/K, cp1004 J/kg/K P wrf.getvar(ds_wrf, pressure) # 单位: hPa T wrf.getvar(ds_wrf, tk) # 单位: K P0 1000.0 R 287.0 cp 1004.0 theta_manual T * (P0 / P) ** (R / cp) # 计算相对涡度 (Relative Vorticity) # 首先获取地图投影下的风分量和科里奥利参数 u wrf.getvar(ds_wrf, ua) # 质量点上的东向风分量 v wrf.getvar(ds_wrf, va) # 质量点上的北向风分量 # 计算涡度 (使用wrf-python的梯度函数它考虑了地图投影变形) # 这比手动计算dx, dy要准确得多 vort wrf.vorticity(u, v)4.3 高效处理大型数据集与时间序列分析WRF输出文件动辄几十GB一次性读入内存不现实。xarray与dask的集成提供了完美的解决方案。# 使用open_mfdataset并启用dask延迟计算 # 指定chunks参数将数据在逻辑上分块不立即加载 ds_big xr.open_mfdataset(wrfout_d01_*.nc, combineby_coords, parallelTrue, chunks{Time: 10, south_north: 100, west_east: 100}) print(ds_big) # 此时数据并未真正加载只是创建了计算任务图。 # 进行一个耗时的计算例如计算整个时间序列的区域平均温度 # 这个计算会被dask分解到各个数据块上并行执行 mean_t2_series ds_big[T2].mean(dim[south_north, west_east]) # 现在才触发实际计算 mean_t2_series_computed mean_t2_series.compute() # 将结果保存到小文件中避免重复计算 mean_t2_series_computed.to_netcdf(domain_mean_T2.nc) # 对于更复杂的操作可以自定义函数并应用 def calculate_aqi(pm25_conc): 一个简化的PM2.5 AQI计算函数示例 # 根据浓度分段线性计算AQI # ... 具体计算逻辑 ... return aqi # 假设有PM2.5浓度变量 pm25 wrf.getvar(ds_big, PM2_5_DRY) # 注意变量名随化学机制而异 # 使用apply_ufunc进行逐元素计算并利用dask并行 aqi xr.apply_ufunc(calculate_aqi, pm25, daskparallelized, output_dtypes[np.float32]) aqi_computed aqi.compute()注意事项使用dask时chunks的大小设置是关键。太小会导致任务调度开销巨大太大则可能内存不足。一个经验法则是每个数据块的大小应在10MB到100MB之间。可以通过ds_big.nbytes / 1e6查看数据总大小再除以设定的块数来估算。5. 可视化与成果输出从数据到图表数据分析的最终目的是为了理解和展示。基于Cartopy和Matplotlib我们可以制作出版级别的图表。5.1 绘制二维场与叠加地图import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from wrf import get_cartopy, latlon_coords # 选择一个时间点 slp_plot slp.isel(Time0) # 获取该变量的经纬度坐标和地图投影信息 lats, lons latlon_coords(slp_plot) cart_proj get_cartopy(slp_plot) # 创建图形和坐标系 fig plt.figure(figsize(12, 8)) ax plt.axes(projectioncart_proj) # 设置地图范围通常就是模式域的范围 ax.set_xlim(wrf.cartopy_xlim(slp_plot)) ax.set_ylim(wrf.cartopy_ylim(slp_plot)) # 添加地理特征 ax.add_feature(cfeature.COASTLINE, linewidth0.8) ax.add_feature(cfeature.BORDERS, linewidth0.5, linestyle:) ax.add_feature(cfeature.STATES, linewidth0.3, linestyle:) ax.add_feature(cfeature.LAKES, alpha0.5) ax.add_feature(cfeature.RIVERS, linewidth0.5) # 绘制海平面气压填色图 contourf ax.contourf(lons, lats, slp_plot, levels20, transformccrs.PlateCarree(), cmapRdBu_r) # 叠加等值线 contour ax.contour(lons, lats, slp_plot, levels10, colorsk, linewidths0.5, transformccrs.PlateCarree()) ax.clabel(contour, inlineTrue, fontsize8, fmt%1.0f) # 绘制风矢以10米风为例 u10_plot u10.isel(Time0) v10_plot v10.isel(Time0) # 每隔N个点取一个风矢避免过于密集 stride 10 ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride], u10_plot.values[::stride, ::stride], v10_plot.values[::stride, ::stride], transformccrs.PlateCarree(), scale300, colorgreen) # 添加色标、标题 plt.colorbar(contourf, axax, orientationhorizontal, pad0.05, aspect40, labelSea Level Pressure (hPa)) ax.set_title(fSea Level Pressure and 10m Wind\n{slp_plot.Time.dt.strftime(%Y-%m-%d %H:%M UTC).values}) plt.tight_layout() plt.savefig(slp_wind_map.png, dpi300, bbox_inchestight) plt.show()5.2 绘制垂直剖面与时空序列图垂直剖面对于分析天气系统或污染物的垂直结构至关重要。# 假设我们想画一条从点A(lat1, lon1)到点B(lat2, lon2)的垂直剖面 cross_start (40.0, 115.0) # (lat, lon) 起点 cross_end (35.0, 120.0) # (lat, lon) 终点 # 使用wrf-python的剖面函数 z_cross wrf.vertcross(z, p, wrf_start_pointcross_start, wrf_end_pointcross_end, latlonTrue, metaTrue) tk_cross wrf.vertcross(tk, p, wrf_start_pointcross_start, wrf_end_pointcross_end, latlonTrue, metaTrue) # 获取剖面线的水平距离坐标 cross_lats wrf.to_np(z_cross.coords[xy_loc]).T[1] # 纬度 cross_lons wrf.to_np(z_cross.coords[xy_loc]).T[0] # 经度 cross_dist wrf.to_np(z_cross.coords[xy_loc]).T[2] # 距离起点距离 (km) fig, ax plt.subplots(figsize(15, 5)) # 绘制温度剖面填色图 pc ax.pcolormesh(cross_dist, wrf.to_np(z_cross.coords[vertical]), wrf.to_np(tk_cross), cmapSpectral_r, shadingauto) # 叠加位势高度等值线 contour ax.contour(cross_dist, wrf.to_np(z_cross.coords[vertical]), wrf.to_np(z_cross), colorsk, linewidths0.8) ax.clabel(contour, inlineTrue, fontsize8, fmt%1.0f) ax.set_xlabel(Distance from Start (km)) ax.set_ylabel(Height (m)) ax.set_title(fTemperature Cross Section\nStart {cross_start} - End {cross_end}) plt.colorbar(pc, axax, labelTemperature (K)) # 在x轴下方标记经纬度 ax_btm ax.twiny() ax_btm.set_xlim(ax.get_xlim()) xticks ax.get_xticks() xtick_lats np.interp(xticks, cross_dist, cross_lats) xtick_lons np.interp(xticks, cross_dist, cross_lons) xtick_labels [f({lat:.1f}N,\n{lon:.1f}E) for lat, lon in zip(xtick_lats, xtick_lons)] ax_btm.set_xticks(xticks) ax_btm.set_xticklabels(xtick_labels, fontsize8) ax_btm.set_xlabel(Latitude, Longitude) plt.tight_layout() plt.savefig(vertical_cross_section.png, dpi300)6. 实战问题排查与脚本优化经验在实际操作中你会遇到各种各样的问题。下面是一些常见坑点和解决技巧。6.1 常见错误与解决方法问题现象可能原因解决方案ImportError: libnetcdf.so.xx: cannot open shared object fileNetCDF库的底层C依赖未正确安装或路径不对。使用conda install -c conda-forge netcdf4重新安装。确保conda环境已激活。检查LD_LIBRARY_PATH环境变量。wrf.getvar返回None或报错变量名错误或该变量在当前WRF输出中不存在。使用wrf.variables()查看文件中所有变量。WRF-Chem的化学变量名取决于机制如PM2_5_DRYMADE/SORGAM、o3等。插值到等压面时出现大量NaN值请求的气压层超出了模式层的气压范围如地面气压只有980hPa却请求了1000hPa层。检查模式域的最低层气压。使用wrf.getvar(ds, ‘pressure’)查看气压范围。插值时设置missing np.nan并处理NaN。Cartopy绘图时地图变形或错位地图投影cart_proj与数据投影不匹配。务必使用wrf.get_cartopy(var)获取数据自带的地图投影对象而不是自己定义。处理大量文件时内存爆炸没有使用分块chunk或延迟加载。使用xr.open_mfdataset(..., chunks{})。对于聚合操作如mean,sum先分块计算再合并结果。避免在循环中重复打开文件。时间坐标解析错误WRF输出文件的时间属性格式不标准。使用wrf.extract_times(ds, timeidx0)来获取正确的时间信息。或者用xr.decode_cfTrue打开数据集。化学物种浓度单位混乱WRF-Chem输出单位可能是mol mol-1 dry air,ppm,μg m-3等。仔细查阅你所使用的化学机制文档如MOZART, RADM2, SAPRC99的输出说明。使用wrf.getvar提取的变量通常会自动转换为标准单位如K, hPa, m/s。6.2 脚本性能优化技巧向量化操作永远避免在Python中对大型数组使用显式循环for循环。尽量使用NumPy或xarray的向量化函数和ufunc。利用Dask进行并行计算对于多文件、多时间步的处理open_mfdataset配合chunks和.compute()是最佳实践。可以设置parallelTrue加速文件读取。选择性读取如果只需要少数几个变量或特定时间步使用xr.open_dataset(..., drop_variables[...])或preprocess函数在打开时即丢弃不需要的数据。缓存中间结果对于耗时的计算步骤如重网格化、复杂的诊断量计算将结果保存为NetCDF文件。下次直接读取中间文件避免重复计算。使用更高效的文件格式如果频繁读写可以考虑将处理后的数据保存为zarr格式它对分块存储和并行读写更友好。6.3 构建可复用的处理流水线为了提升效率建议将通用功能模块化# my_wrf_tools.py import xarray as xr import wrf import numpy as np class WRFProcessor: def __init__(self, wrfout_paths): self.ds xr.open_mfdataset(wrfout_paths, combineby_coords, chunks{Time: 24}) self.cache {} # 用于缓存计算昂贵的变量 def get_diagnostic_var(self, var_name): 获取诊断变量如果缓存中有则直接返回 if var_name in self.cache: return self.cache[var_name] if var_name pv: # 计算位涡 p wrf.getvar(self.ds, pressure) theta wrf.getvar(self.ds, theta) ua, va, wa wrf.getvar(self.ds, uvmet) pv wrf.potential_vorticity(p, theta, ua, va, wa) self.cache[pv] pv return pv elif var_name precip_total: # 计算总降水量 rainc self.ds[RAINC] rainnc self.ds[RAINNC] total rainc rainnc self.cache[precip_total] total return total else: # 尝试用wrf.getvar直接提取 var wrf.getvar(self.ds, var_name) self.cache[var_name] var return var def save_to_disk(self, var_list, output_file): 将指定变量列表保存到新的NetCDF文件 ds_out xr.Dataset() for var_name in var_list: ds_out[var_name] self.get_diagnostic_var(var_name) ds_out.to_netcdf(output_file) print(fSaved to {output_file}) # 主脚本 main.py from my_wrf_tools import WRFProcessor processor WRFProcessor([wrfout_d01_*]) # 一次性计算并保存多个常用诊断量 processor.save_to_disk([slp, T2, rh2, uvmet10, precip_total], processed_diagnostics.nc)这套脚本体系的核心思想是自动化和可重复。一旦搭建完成对于新的模拟任务你只需要修改输入文件路径和少数参数就能快速得到分析就绪的数据和图表将精力从繁琐的数据处理中解放出来更多地投入到科学问题的探索上。本文还有配套的精品资源点击获取