简介：本资源是一套面向气象与环境科研人员、大气科学研究生及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 python=3.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安装：
注意，wrf-python依赖于netCDF4和xarray，但通常用pip安装时会自动处理。pip install wrf-python - cartopy：专业的地图绘图库。相比Basemap（已停止维护），Cartopy更现代、功能更强大，与matplotlib集成度极高，用于绘制带有地图投影的天气图、绘制模式区域、叠加海岸线、国界等。
- netCDF4：提供底层API，当需要极高性能或操作xarray未覆盖的底层属性时使用。
2.3 集成开发环境（IDE）配置
对于此类涉及大量数据探索和脚本编写的工作，一个好的IDE至关重要。Visual Studio Code (VSCode)是目前的最佳选择之一。
- 安装VSCode及Python扩展：从官网安装VSCode，然后在扩展商店搜索并安装“Python”扩展（由Microsoft发布）。
- 选择解释器：在VSCode中，按
Ctrl+Shift+P，输入“Python: Select Interpreter”，选择我们刚才创建的wrf_env环境下的Python路径（通常类似~/miniconda3/envs/wrf_env/bin/python）。 - 配置Jupyter内核：在
wrf_env环境中运行python -m ipykernel install --user --name=wrf_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(days=i) for i in range((end_date - start_date).days + 1)] target_dates_np = np.array(target_dates, dtype='datetime64[D]') # 使用线性插值（对于排放，更复杂的方法可能考虑周末/工作日因子） # 这里需要将时间转换为数值进行插值，简单演示思路 # 实际应用中，建议对每个网格点的时间序列进行插值 nox_daily = nox_orig.interp(time=target_dates_np, method='linear') # 3. 空间重映射（从经纬度网格到WRF网格） # 首先，需要你的WRF域定义文件（通常是geo_em.d01.nc）中的经纬度变量 ds_wrf = xr.open_dataset('geo_em.d01.nc') wrf_lat = ds_wrf['XLAT_M'].isel(Time=0).values # 网格中心纬度 wrf_lon = ds_wrf['XLONG_M'].isel(Time=0).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, method='linear', bounds_error=False, fill_value=0.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(time=0) 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)}, axis=0) * 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, combine='by_coords', engine='netcdf4') # 提取海平面气压、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, R=287 J/kg/K, cp=1004 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', combine='by_coords', parallel=True, 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, dask='parallelized', 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(Time=0) # 获取该变量的经纬度坐标和地图投影信息 lats, lons = latlon_coords(slp_plot) cart_proj = get_cartopy(slp_plot) # 创建图形和坐标系 fig = plt.figure(figsize=(12, 8)) ax = plt.axes(projection=cart_proj) # 设置地图范围（通常就是模式域的范围） ax.set_xlim(wrf.cartopy_xlim(slp_plot)) ax.set_ylim(wrf.cartopy_ylim(slp_plot)) # 添加地理特征 ax.add_feature(cfeature.COASTLINE, linewidth=0.8) ax.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':') ax.add_feature(cfeature.STATES, linewidth=0.3, linestyle=':') ax.add_feature(cfeature.LAKES, alpha=0.5) ax.add_feature(cfeature.RIVERS, linewidth=0.5) # 绘制海平面气压填色图 contourf = ax.contourf(lons, lats, slp_plot, levels=20, transform=ccrs.PlateCarree(), cmap='RdBu_r') # 叠加等值线 contour = ax.contour(lons, lats, slp_plot, levels=10, colors='k', linewidths=0.5, transform=ccrs.PlateCarree()) ax.clabel(contour, inline=True, fontsize=8, fmt='%1.0f') # 绘制风矢（以10米风为例） u10_plot = u10.isel(Time=0) v10_plot = v10.isel(Time=0) # 每隔N个点取一个风矢，避免过于密集 stride = 10 ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride], u10_plot.values[::stride, ::stride], v10_plot.values[::stride, ::stride], transform=ccrs.PlateCarree(), scale=300, color='green') # 添加色标、标题 plt.colorbar(contourf, ax=ax, orientation='horizontal', pad=0.05, aspect=40, label='Sea Level Pressure (hPa)') ax.set_title(f'Sea 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', dpi=300, bbox_inches='tight') 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_point=cross_start, wrf_end_point=cross_end, latlon=True, meta=True) tk_cross = wrf.vertcross(tk, p, wrf_start_point=cross_start, wrf_end_point=cross_end, latlon=True, meta=True) # 获取剖面线的水平距离坐标 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), cmap='Spectral_r', shading='auto') # 叠加位势高度等值线 contour = ax.contour(cross_dist, wrf.to_np(z_cross.coords['vertical']), wrf.to_np(z_cross), colors='k', linewidths=0.8) ax.clabel(contour, inline=True, fontsize=8, fmt='%1.0f') ax.set_xlabel('Distance from Start (km)') ax.set_ylabel('Height (m)') ax.set_title(f'Temperature Cross Section\nStart {cross_start} -> End {cross_end}') plt.colorbar(pc, ax=ax, label='Temperature (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, fontsize=8) ax_btm.set_xlabel('Latitude, Longitude') plt.tight_layout() plt.savefig('vertical_cross_section.png', dpi=300)6. 实战问题排查与脚本优化经验
在实际操作中，你会遇到各种各样的问题。下面是一些常见坑点和解决技巧。
6.1 常见错误与解决方法
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
ImportError: libnetcdf.so.xx: cannot open shared object file | NetCDF库的底层C依赖未正确安装或路径不对。 | 使用conda install -c conda-forge netcdf4重新安装。确保conda环境已激活。检查LD_LIBRARY_PATH环境变量。 |
wrf.getvar返回None或报错 | 变量名错误，或该变量在当前WRF输出中不存在。 | 使用wrf.variables()查看文件中所有变量。WRF-Chem的化学变量名取决于机制，如PM2_5_DRY（MADE/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, timeidx=0)来获取正确的时间信息。或者用xr.decode_cf=True打开数据集。 |
| 化学物种浓度单位混乱 | 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()是最佳实践。可以设置parallel=True加速文件读取。 - 选择性读取：如果只需要少数几个变量或特定时间步，使用
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, combine='by_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(f'Saved 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')这套脚本体系的核心思想是自动化和可重复。一旦搭建完成，对于新的模拟任务，你只需要修改输入文件路径和少数参数，就能快速得到分析就绪的数据和图表，将精力从繁琐的数据处理中解放出来，更多地投入到科学问题的探索上。
本文还有配套的精品资源，点击获取
