简介本资源是一套面向气象与环境科研人员、大气科学方向研究生及Python气象数据处理初学者的WRF/WRF-Chem全流程自动化脚本集聚焦模型前处理地形/土地利用生成、初始场插值、namelist配置与后处理NetCDF解析、时空可视化、特征导出两大核心环节。压缩包共95个文件含62个Python主脚本覆盖ERA5数据下载、TsinghuaLU/VIIRS/Landsat8遥感数据预处理、垂直剖面绘制、泰勒图生成、UCM城市参数构建等、15个备份脚本.zbak、5张示例图表含geogrid、plaindraw、泰勒图等关键输出效果以及字体、许可证、README等配套文件总大小17.17MB。已有37人学习下载。用户可直接复用近百个模块化脚本快速完成从原始遥感数据→WRF输入文件→模拟结果分析→AI建模特征提取的完整链路尤其适配Scikit-learn特征工程需求显著降低WRF-Chem数据处理门槛与重复劳动成本。1. 项目概述从“跑通”到“跑好”的必经之路如果你正在或即将使用WRFWeather Research and Forecasting或WRF-Chem模型进行气象或大气化学模拟那么你大概率会遇到一个比模型本身运行更耗时、更令人头疼的环节数据的前后处理。模型本身像一台精密的发动机但给它喂什么油输入数据以及如何解读它排出的尾气输出数据才是决定你科研或业务工作成败的关键。这就是Python脚本大显身手的地方。我接触WRF系列模型超过十年从最初的手动修改namelist、用NCL画图到后来全面转向Python自动化流程深刻体会到一套高效、可靠的数据处理脚本是多么重要。它不仅仅是节省时间更是保证数据质量、实现可重复性研究、以及进行复杂分析比如化学物种追踪、源解析的基石。很多人把大量精力花在调试模型物理参数上却忽略了前后处理流程的规范化和自动化导致结果难以复现或者因为某个手动步骤的疏漏而前功尽弃。这个项目就是聚焦于用Python构建一套针对WRF/WRF-Chem模型的“预处理”与“后处理”脚本体系。预处理负责将原始再分析数据如FNL、GFS、ERA5或排放清单转换成WRF能“消化”的格式如WPS处理的中间文件或直接准备化学初始边界条件后处理则是将WRF输出的NetCDF文件转换为我们需要的分析结果如特定变量的时空分布图、垂直剖面、站点时间序列、轨迹计算、化学物种收支分析等。我们将深入每个环节不仅告诉你“怎么做”更会解释“为什么这么做”并分享那些在官方文档里找不到的实战经验和避坑指南。2. 核心需求与工具栈解析为什么是Python在开始动手之前我们需要明确核心需求并理解为什么Python是目前处理WRF数据生态中的最优选。2.1 核心需求拆解一个完整的WRF/WRF-Chem模拟流程对数据处理脚本的需求可以归纳为以下几点自动化与批处理模拟往往是多案例、多参数、长时间序列的。手动点击操作不现实脚本必须能自动循环处理不同日期、不同区域、不同排放情景的数据。格式转换与重映射这是预处理的核心。需要将全球经纬度网格的再分析数据通过WPSWRF Preprocessing System插值到我们定义的高分辨率区域网格上。对于WRF-Chem还需要处理排放清单如EDGAR、MEIC将其从原始格式如NetCDF、文本转换为WRF-Chem所需的wrfchemi_*文件网格化、时间解析的排放源。数据提取与子集化WRF输出文件通常包含所有变量、所有层次、所有时间步文件巨大。后处理脚本需要能快速、准确地提取感兴趣的变量如温度、降水、PM2.5、O3、层次如地面、850hPa、500hPa和区域如某个城市群。计算与诊断很多需要的量不是直接输出的需要基于原始变量进行计算。例如计算位温、相当位温、水汽通量、大气稳定度指数如CAPE或者对于WRF-Chem计算气溶胶光学厚度AOD、二次有机气溶胶SOA产量、物种的干湿沉降通量等。可视化与制图生成符合学术出版或业务报告要求的高质量图表包括空间填色图、等值线图、垂直剖面图、时间序列图、散点图等。质量控制与错误处理流程中任何一个环节出错如下载数据不完整、WPS运行失败、磁盘空间不足脚本应能及时发现、记录日志并尽可能优雅地处理或中断避免浪费计算资源。2.2 为什么选择Python及其工具栈早年WRF社区后处理多依赖NCLNCAR Command Language和GrADS。它们功能强大但在代码可维护性、生态系统丰富度、学习曲线和社区活跃度上已逐渐被Python超越。丰富的科学计算库NumPy提供高效的数组计算xarray是处理NetCDF数据的“神器”它引入了带标签的多维数组概念操作NetCDF就像操作Pandas DataFrame一样直观完美替代了NCL的大部分功能pandas处理时间序列和表格数据得心应手。强大的可视化能力Matplotlib是基石高度可定制Cartopy是新一代的地图绘图库替代了古老的Basemap支持更多投影和更现代的地理数据处理对于交互式探索Plotly或HvPlot是很好的选择。专门的WRF处理包wrf-python是UCAR官方维护的Python工具包提供了一系列函数用于从WRF输出文件中提取变量、插值到气压层、计算诊断量如涡度、散度、降水类型等极大简化了后处理。卓越的自动化与集成能力Python可以轻松调用系统命令subprocess模块来运行WPS、real.exe、wrf.exe可以编写爬虫自动下载FNL/GFS数据也可以与任务调度系统如Slurm、PBS结合实现端到端的自动化模拟流程。活跃的社区与可复用性代码易于模块化、函数化方便团队共享和积累。Jupyter Notebook为交互式分析和教学提供了绝佳环境。注意虽然Python是主力但并不意味着完全抛弃WPS。WPS中的geogrid、ungrib、metgrid程序仍然是预处理不可或缺的核心。我们的Python脚本角色是“编排者”和“增强者”负责准备输入、控制流程、处理WPS生成的数据以及进行后续深度分析。3. 预处理脚本详解为模型准备“食材”预处理是模拟的起点目标是为real.exeWRF或real.execonvert_emiss.exeWRF-Chem生成正确的输入文件。我们将流程分解并实现自动化。3.1 驱动数据下载与时间管理WRF需要初始场和边界条件通常来自NCEP FNL、GFS或ECMWF ERA5再分析数据。手动下载费时费力。脚本核心功能自动生成下载列表根据模拟的起止时间start_date,end_date和时效如GFS每6小时一次生成所有需要下载的文件时间列表。调用下载工具使用subprocess调用wget、curl或专门的客户端如cdsapifor ERA5rda-apifor FNL进行批量下载。这里需要处理认证API key和网络错误重试。完整性校验下载后检查文件大小或使用ncdump -h快速检查文件头确保数据完整。import subprocess from datetime import datetime, timedelta def download_gfs_data(start_dt, end_dt, interval_hours6, output_dir./gfs_data): 自动下载GFS数据示例需预先配置好wget或ftp base_url ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.{date:%Y%m%d/%H}/ os.makedirs(output_dir, exist_okTrue) current_dt start_dt while current_dt end_dt: # 构造文件名例如 gfs.t00z.pgrb2.0p25.f000 # 实际文件名更复杂需要根据分辨率确定 file_pattern fgfs.t{current_dt.hour:02d}z.pgrb2.0p25.f000 full_url base_url.format(datecurrent_dt) file_pattern cmd [wget, -c, -N, -P, output_dir, full_url] try: subprocess.run(cmd, checkTrue, timeout300) print(f成功下载: {file_pattern}) except subprocess.CalledProcessError as e: print(f下载失败 {file_pattern}: {e}) # 这里可以加入重试逻辑 except subprocess.TimeoutExpired: print(f下载超时 {file_pattern}) current_dt timedelta(hoursinterval_hours)实操心得使用-c断点续传和-N只下载比本地新的文件是下载大文件时的好习惯避免重复下载。ERA5数据通过ECMWF的CDS API下载是最佳途径需要先注册获取API key。使用cdsapi库时注意请求的数据量过大可能导致排队脚本中需要加入等待和状态查询逻辑。本地缓存对于公共数据建议在服务器或工作站建立本地缓存目录所有项目共享避免重复下载占用带宽和磁盘。3.2 WPS流程自动化封装运行WPSgeogrid - ungrib - metgrid通常涉及多次修改namelist.wps和执行命令。我们可以用Python脚本统一管理。脚本核心功能动态生成namelist.wps将域设置、时间范围、路径参数等写入一个模板字典然后用f90nml库专门读写Fortran namelist的Python库或简单的字符串替换生成最终的namelist.wps文件。这比手动编辑更不易出错尤其在做参数敏感性试验时。顺序执行与错误监控用subprocess依次运行geogrid.exe、ungrib.exe、metgrid.exe。关键是要检查每个步骤的日志文件如geogrid.log或退出状态码一旦出错就立即停止并报错而不是继续执行下去产生无用的中间文件。清理中间文件ungrib会生成大量的FILE:*文件在确认metgrid成功后脚本可以自动清理这些临时文件释放磁盘空间。import f90nml import subprocess import shutil def run_wps(domain_config, time_config, wps_dir./WPS): 自动化运行WPS流程 domain_config: 字典包含max_dom, dx, dy, parent_id等域信息 time_config: 字典包含start/end_date, interval_seconds等时间信息 # 1. 准备namelist.wps namelist_template { share: { wrf_core: ARW, max_dom: domain_config[max_dom], start_date: time_config[start_date], # 格式: 2023-01-01_00:00:00 end_date: time_config[end_date], interval_seconds: time_config[interval_seconds], }, geogrid: { ... }, # 地理数据路径、地图参数等 ungrib: { ... }, metgrid: { ... }, } namelist_path os.path.join(wps_dir, namelist.wps) f90nml.write(namelist_template, namelist_path) # 2. 链接Vtable (根据数据源选择) vtable_src fVtable.{data_source} # 如 Vtable.GFS vtable_dst os.path.join(wps_dir, Vtable) if os.path.exists(vtable_src): if os.path.lexists(vtable_dst): os.unlink(vtable_dst) os.symlink(vtable_src, vtable_dst) else: raise FileNotFoundError(fVtable文件 {vtable_src} 不存在) # 3. 运行geogrid os.chdir(wps_dir) try: subprocess.run([./geogrid.exe], checkTrue, stdoutsubprocess.PIPE, stderrsubprocess.PIPE) print(geogrid.exe 执行成功) except subprocess.CalledProcessError as e: print(fgeogrid.exe 执行失败: {e}) with open(geogrid.log, r) as f: print(f.read()) return False # 4. 运行ungrib # ... 类似逻辑先链接GRIB数据文件 # 5. 运行metgrid # ... # 6. 可选清理FILE:*文件 # if keep_intermediate is False: # for f in glob.glob(FILE:*): # os.remove(f) return True注意事项路径问题geogrid.exe需要访问静态地理数据GEOG目录确保路径在namelist中设置正确并且该目录有足够权限和空间。Vtable匹配不同的驱动数据GFS, FNL, ERA5需要不同的Vtable文件链接错误会导致ungrib读不出变量。内存与磁盘高分辨率、多层次的模拟metgrid可能会生成巨大的中间文件确保/tmp或工作目录有足够空间TB级别。3.3 WRF-Chem排放清单处理专项这是WRF-Chem预处理中最具挑战性的部分之一。排放清单通常以年度全球网格如0.1°x0.1°或区域列表形式提供需要将其处理成WRF-Chem能识别的、时间分辨率更高如逐小时、空间上与模拟域匹配的NetCDF文件。核心步骤与脚本实现清单读取与合并使用xarray打开NetCDF格式的排放清单如EDGAR或用pandas读取CSV文本清单。可能需要合并人为源、生物源、海盐、沙尘等不同来源。时空重映射空间使用xesmf库一个基于ESMF的Python重映射工具或scipy插值将清单从原始网格重映射到WRF模拟域通过geo_em.d0x.nc文件定义的不规则网格。这比简单的最近邻插值更准确特别是对化学物种。时间排放具有日变化、周变化和季节变化。脚本需要根据预设的时间变化因子如交通排放的早晚高峰、工业排放的昼夜差异将年均或月均排放量分解到每个模拟小时。这通常通过一个时间剖面文件来实现。物种映射与分割排放清单中的物种如NOx, SO2, NMVOC需要映射到WRF-Chem化学机制如RADM2, MOZART, SAPRC99中的具体物种。同时需要将VOC按照化学机制要求的比例分割成具体的模型物种如PAR, OLE, TOL等。生成wrfchemi文件调用convert_emiss.exe需要编译进WRF-Chem或使用anthro_emis等预处理程序但更灵活的方式是直接用xarray和netCDF4库按照wrfchemi的格式要求从头创建NetCDF文件。这要求你深刻理解该文件的维度Time, emissions_zdim, south_north, west_east和变量属性。import xarray as xr import xesmf as xe def process_emission_to_wrfchem(orig_emis_nc, geo_em_file, chem_mechanismRADM2, output_prefixwrfchemi): 将原始排放清单处理为WRF-Chem输入文件简化示例 # 1. 读取原始排放数据和WRF地理网格文件 ds_emis xr.open_dataset(orig_emis_nc) ds_geo xr.open_dataset(geo_em_file) # 2. 创建从原始网格到WRF网格的重映射器 (使用双线性插值) # 注意排放重映射需谨慎质量守恒是关键。xesmf提供守恒插值选项。 regridder xe.Regridder(ds_emis, ds_geo, bilinear, reuse_weightsTrue) # 3. 对每种排放物种进行重映射 emis_vars [NOx, SO2, CO, VOC] # 示例变量名 emis_regridded {} for var in emis_vars: if var in ds_emis: emis_regridded[var] regridder(ds_emis[var]) # 4. 应用时间变化因子此处简化实际需根据小时、星期几计算 # 假设我们有一个24小时的时间变化因子数组 time_profile [24] # 需要将年均排放量根据模拟的每个小时进行分配 # 这是一个复杂的步骤涉及日期时间计算和乘法运算 # 5. 物种映射与分割 (以VOC为例) # 假设 chem_mechanism_mapping 是一个字典定义了原始VOC到机制物种的分配系数 # 例如{PAR: 0.4, OLE: 0.3, TOL: 0.1, ...} # emis_regridded[PAR] emis_regridded[VOC] * 0.4 # ... # 6. 创建符合wrfchemi格式的Dataset # 需要严格按照WRF-Chem要求的维度、变量名、属性和单位来构建 # 此处省略详细的构建过程它涉及创建时间维度、垂直层次通常为1层等 # 7. 写入NetCDF文件 # ds_wrfchemi.to_netcdf(f{output_prefix}_00z_d01) print(排放处理流程完成示例)避坑技巧质量守恒检查重映射前后计算全球或区域的总排放量确保误差在可接受范围内1%。xesmf的conservative插值方法可以更好地保持总量。垂直分配大部分清单只提供地面排放。对于高空点源如电厂烟囱需要根据烟囱参数高度、温度、速度进行垂直分配这通常需要额外的处理脚本或使用PREP-CHEM-SRC等工具。背景化学场除了排放WRF-Chem还需要化学初始边界条件如wrfinput中的化学物种浓度。对于全球化学模拟如MOZART提供的数据处理方式与气象场类似但变量更多。确保ungrib使用的Vtable能识别化学变量或者使用专门的工具如mozbc来处理。4. 后处理脚本详解从海量输出中提炼“洞察”WRF模拟完成后输出是庞大的NetCDF文件集wrfout_d0*_*。后处理的目标是将其转化为有科学或业务价值的图形和数据集。4.1 高效数据提取与子集化直接使用NetCDF4库或xarray打开所有wrfout文件进行全局操作极易导致内存溢出。必须采用惰性加载和分块处理的策略。脚本策略使用xarray的延迟计算xarray.open_mfdataset函数可以同时打开多个文件并创建一个逻辑上统一的数据集而不会立即将数据读入内存。结合chunks参数指定分块大小可以与Dask库集成实现并行计算。选择性提取在打开数据集时或之后立即使用.sel()选择和.isel()索引选择方法提取需要的Time、bottom_top或bottom_top_stag、south_north、west_east维度子集。变量计算后再持久化对于需要计算的诊断变量利用xarray的向量化运算在“虚拟”数据集上进行最后对结果调用.compute()或.persist()才真正触发计算并加载到内存。对于最终需要反复使用的中间结果可以保存到新的、更小的NetCDF文件中。import xarray as xr import numpy as np def extract_and_calc_pm25(wrfout_files, domain1, save_path./processed/pm25.nc): 从一系列wrfout文件中提取并计算PM2.5浓度适用于包含气溶胶物种的机制 假设使用MOSAIC或GOCART气溶胶方案 # 1. 延迟打开多个文件 # 使用 chunks{Time: 10, south_north: 100, west_east: 100} 进行分块 ds xr.open_mfdataset(wrfout_files, chunks{Time: 24}, combineby_coords, parallelTrue) # 2. 选择特定域 (假设变量有‘XTIME’等维度标识) # 更常见的是直接打开对应domain的文件如 wrfout_d01_* # 这里我们假设ds已经只包含一个域的数据 # 3. 计算PM2.5 (以GOCART机制为例PM2.5近似为部分物种之和) # 注意不同化学机制PM2.5的算法不同此处仅为示例。 # 例如可能包括硫酸盐(so4)、硝酸盐(no3)、铵盐(nh4)、有机碳(oc)、黑碳(bc)、海盐(ss)、沙尘(du)的细模态部分 # 需要根据实际输出变量名调整 species_list [so4_a01, no3_a01, nh4_a01, oc_a01, bc_a01] # 细模态气溶胶 pm25 None for species in species_list: if species in ds: if pm25 is None: pm25 ds[species] else: pm25 pm25 ds[species] if pm25 is not None: # 添加单位等属性 pm25.attrs[long_name] Fine Particulate Matter (PM2.5) pm25.attrs[units] ug m-3 # 4. 将pm25作为一个新变量放回数据集或创建新数据集 ds_pm25 pm25.to_dataset(namePM2.5) # 5. 保存到新文件触发实际计算和写入 ds_pm25.to_netcdf(save_path) print(fPM2.5数据已保存至 {save_path}) return ds_pm25 else: print(未找到计算PM2.5所需的物种变量。) return None注意事项变量名差异WRF输出变量名与所选物理参数化方案、化学机制紧密相关。务必先用ncdump -h或xarray查看变量列表。wrf-python的getvar函数可以屏蔽部分差异通过通用名称如‘temp’ ‘uvmet’获取变量。内存管理对于超大型数据集如高分辨率、长时间模拟即使分块也可能在计算聚合操作如时间平均、空间平均时内存暴涨。考虑分步处理先按时间片处理保存中间结果再对中间结果进行聚合。时间处理WRF输出时间可能是模拟分钟数或日期字符串。使用wrf-python的getvar函数配合time关键字或使用xarray的cftime支持可以将其转换为Python的datetime对象方便后续分析。4.2 诊断量与可视化计算除了直接输出变量许多重要物理量需要计算。wrf-python库封装了大部分常用诊断量的计算。常用诊断量计算示例import wrf import numpy as np def calculate_meteorological_diagnostics(ds): 计算一些常用的气象诊断量 ds: 通过wrf-python或xarray打开的wrfout数据集 # 注意wrf.getvar 需要原始文件路径或已加载的netCDF4 Dataset对象 # 这里假设ds是netCDF4.Dataset对象 # 如果ds是xarray Dataset可能需要from netCDF4 import Dataset; ncfile Dataset(wrfout_file) # 获取经纬度坐标用于绘图 lats, lons wrf.latlon_coords(ds) # 获取地面气压、温度、露点温度等 slp wrf.getvar(ds, slp) # 海平面气压 t2 wrf.getvar(ds, T2) # 2米温度 td2 wrf.getvar(ds, td2) # 2米露点温度 u10, v10 wrf.getvar(ds, uvmet10) # 10米风地图投影坐标 # 计算相对湿度 rh2 wrf.getvar(ds, rh2) # 2米相对湿度 # 计算降水量累计降水转瞬时或逐小时 rainc wrf.getvar(ds, RAINC) # 对流降水 rainnc wrf.getvar(ds, RAINNC) # 非对流降水 total_rain rainc rainnc # 总累计降水 # 计算逐小时降水需要时间维度 # hour_rain total_rain.diff(dimTime) # 注意单位换算 # 计算垂直剖面例如沿一条线的位温剖面 # 假设我们想画一条从点(start_lat, start_lon)到(end_lat, end_lon)的垂直剖面 start_point (35.0, 115.0) # (lat, lon) end_point (40.0, 120.0) # 提取剖面线上的数据 # 注意wrf.vertcross 需要等压面或高度面上的变量需要先插值 p wrf.getvar(ds, pressure) # 获取气压场 tk wrf.getvar(ds, tk) # 温度 (K) # 将温度插值到气压面上 # 这里只是示意实际调用较复杂 # theta_cross wrf.vertcross(theta, p, wrf_line, latlonTrue) diagnostics { slp: slp, t2: t2, td2: td2, u10: u10, v10: v10, rh2: rh2, total_rain: total_rain, lats: lats, lons: lons } return diagnostics可视化脚本核心 计算完成后使用Cartopy和Matplotlib绘图。关键是要正确处理WRF的地图投影信息。import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt from wrf import get_cartopy, cartopy_xlim, cartopy_ylim def plot_surface_map(wrf_file, var_name, time_idx0, level0, domain1): 绘制WRF输出变量的地面填色图 # 从文件获取数据和地图投影信息 ncfile Dataset(wrf_file) var_data wrf.getvar(ncfile, var_name, timeidxtime_idx) # 获取地图投影 cart_proj get_cartopy(var_data) # 创建图形 fig plt.figure(figsize(12, 8)) ax plt.axes(projectioncart_proj) # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale(50m), linewidth0.8) ax.add_feature(cfeature.BORDERS.with_scale(50m), linestyle:, linewidth0.5) ax.add_feature(cfeature.STATES.with_scale(50m), linewidth0.3) # 添加河流、湖泊等... # 设置图形范围 ax.set_xlim(cartopy_xlim(var_data)) ax.set_ylim(cartopy_ylim(var_data)) # 绘制填色图 # 需要将数据坐标转换为地图投影坐标 lats, lons wrf.latlon_coords(var_data) # 使用pcolormesh或contourf # 注意如果var_data是三维带bottom_top需要切片 var_data[level, :, :] if len(var_data.shape) 3: plot_data var_data[level, :, :] else: plot_data var_data # 这里使用pcolormesh对于不规则WRF网格需使用经纬度坐标 mesh ax.pcolormesh(lons, lats, plot_data, transformccrs.PlateCarree(), shadingauto, cmapviridis) # 添加色标 plt.colorbar(mesh, axax, orientationhorizontal, pad0.05, shrink0.8) # 添加标题 plt.title(f{var_data.description} ({var_data.units}) at Time Index {time_idx}) plt.tight_layout() plt.savefig(f{var_name}_d0{domain}_t{time_idx}.png, dpi300) plt.close(fig) ncfile.close()实操心得投影是关键WRF使用自己的地图投影通常在namelist中定义如Lambert Conformal, Mercator等。wrf-python的get_cartopy函数能自动提取这个投影信息确保绘制的图形在空间上是正确的。直接使用PlateCarree投影绘制会导致严重变形。性能优化绘制大量时间步的动画或批量出图时避免在循环中重复创建地图投影和地理特征可以将其移到循环外。使用matplotlib的FuncAnimation或批量生成图片后用ffmpeg合成视频。图形美化Cartopy提供了丰富的地图元素。合理使用cfeature添加海岸线、国界、河流、湖泊并使用gridliner添加经纬度网格和标签能让图更专业。5. 实战中的常见问题与排查技巧即使脚本写得再完善在实际运行中也会遇到各种问题。以下是一些典型问题及其解决思路。5.1 预处理阶段问题问题现象可能原因排查步骤与解决方案ungrib.exe运行失败提示“Not a GRIB file”或“No fields found”1. GRIB文件损坏或下载不完整。2. 链接的Vtable与GRIB数据源不匹配。3. GRIB数据版本如GFS 0.25度 vs 0.5度与Vtable预期不符。1. 用grib_dump或wgrib2检查GRIB文件头。2. 核对数据源GFS/FNL/ERA5并链接正确的Vtable。3. 查看WPS日志ungrib.log看具体是哪个变量读不出。有时需要手动修改Vtable以适配新数据格式。metgrid.exe运行缓慢或内存溢出1. 模拟区域太大或分辨率太高。2. 垂直层数过多。3. 输入的气象场时间频率太高。1. 检查namelist.wps中的e_we/e_sn和dx/dy。2. 减少垂直层数需与WRF的namelist.input一致。3. 增加interval_seconds降低驱动数据的时间频率如从3小时降到6小时。4. 在运行metgrid时尝试增加进程数如果支持MPI。生成的wrfinput文件中某些变量全是缺省值1.metgrid插值过程中某些变量在边界处丢失。2. 静态数据如土地利用类型、土壤类型在geogrid步骤未正确读取。1. 检查metgrid.log是否有警告。可以尝试在namelist.wps的metgrid部分调整interp_option如从average_gcell改为nearest_neighbor。2. 检查GEOG目录路径及权限确保所需的地理数据文件存在。5.2 后处理阶段问题问题现象可能原因排查步骤与解决方案用xarray打开wrfout文件报错或变量找不到1. 文件损坏。2. 变量名因物理方案不同而异。3. 文件格式问题如NetCDF4 classic vs 64-bit offset。1. 用ncdump -h快速查看文件结构和变量列表这是最可靠的方法。2. 使用wrf-python的getvar函数它提供了统一的接口来获取物理量如‘temp’ ‘uvmet’内部处理了变量名差异。3. 确保安装的netCDF4库支持对应的NetCDF版本。计算诊断量如涡度、散度结果异常NaN或极大值1. 数据在网格边界或地形附近存在无效值。2.wrf-python的getvar函数参数使用不当如未指定投影坐标。3. 输入变量单位不匹配。1. 使用xarray的.where()或numpy的np.isnan过滤无效值后再计算。2. 仔细阅读wrf-python官方文档确认函数所需的输入数据格式是原始变量还是getvar提取后的。例如计算相对涡度需要风场在地图投影坐标上。3. 打印输入变量的单位和值范围进行验证。绘图时图形扭曲或坐标错误1. 使用了错误的地图投影。2. 数据经纬度坐标与绘图投影不匹配。3. Cartopy版本与WRF-Python兼容性问题。1.始终坚持使用wrf.get_cartopy()获取投影并用transformccrs.PlateCarree()参数将数据lons, lats转换到地理坐标进行绘制。2. 使用wrf.cartopy_xlim/ylim设置图形范围而不是手动指定。3. 检查wrf-python和Cartopy的版本尽量使用较新且兼容的版本组合。处理大量文件时内存不足1. 一次性加载所有数据。2. 中间变量未及时释放。1.使用xarray.open_mfdataset并设置chunks参数启用Dask进行惰性加载和分块计算。2. 对于循环处理显式关闭不再需要的文件句柄ds.close()。3. 将大任务分解为多个小任务分步处理并保存中间结果。5.3 脚本调试与优化心得日志是生命线在脚本的关键步骤开始、结束、出错点添加详细的日志记录包括时间、操作、关键参数和状态。使用Python的logging模块可以方便地控制日志级别DEBUG, INFO, WARNING, ERROR并输出到文件这对于在后台长时间运行的批处理任务至关重要。增量式开发与测试不要试图一次性写完处理整个流程的巨型脚本。应该模块化开发每个函数只负责一个明确的任务如下载、运行WPS、计算PM2.5、绘图。然后针对一个小的、代表性的案例如单日、低分辨率模拟进行完整流程测试确保每个模块都工作正常。参数配置文件将模拟的域设置、时间范围、路径、化学机制选项等所有可配置参数写进一个单独的配置文件如config.yaml或config.ini。主脚本读取这个配置文件。这样当需要更换模拟案例时只需修改配置文件而无需改动脚本代码大大提高了复用性和可维护性。利用并行加速对于可并行的任务如处理不同时间段的文件、绘制多幅图使用multiprocessing或concurrent.futures模块实现多进程并行。xarray结合Dask可以自动对分块数据进行并行计算。注意平衡任务粒度和进程数避免进程间通信开销过大。6. 从脚本到流程构建自动化模拟系统当预处理和后处理脚本都稳定可靠后我们可以将它们整合构建一个端到端的自动化模拟与后处理系统。这个系统通常由一个主控脚本驱动它协调以下任务参数解析与配置加载读取用户提供的案例配置文件。数据准备调用下载脚本获取驱动数据调用排放处理脚本生成化学输入。WPS/WRF执行依次运行geogrid,ungrib,metgrid,real.exe,wrf.exe可能还有convert_emiss.exe。每一步都进行状态检查失败则发送通知如邮件并中止。后处理流水线WRF运行成功后自动触发后处理脚本进行数据提取、诊断量计算和批量绘图。结果打包与归档将关键的输入文件namelist、日志、核心输出数据和图表打包转移到归档存储并清理临时工作目录以释放空间。这样的系统可以部署在高性能计算集群上通过作业调度系统如Slurm提交实现无人值守的批量模拟实验极大提升科研和生产效率。我个人在构建这类系统时最后通常会写一个简单的Makefile或使用Snakemake/Nextflow这样的工作流管理工具来定义任务之间的依赖关系这比纯Python脚本在管理复杂流水线时更加清晰和强大。例如可以定义“绘图”任务依赖于“计算PM2.5”任务而“计算PM2.5”又依赖于“WRF运行完成”。当某个中间文件更新时工作流引擎可以自动识别并重新运行下游任务。无论采用哪种方式核心思想是一致的将重复、易错的手工操作固化为可靠、可复用的代码把宝贵的精力留给真正的科学问题分析和模型物理机制的探索。这套基于Python的WRF数据处理脚本就是你从模型“使用者”进阶为模型“驾驭者”的关键工具。本文还有配套的精品资源点击获取
Python自动化WRF/WRF-Chem前后处理：从数据准备到结果分析实战指南
简介本资源是一套面向气象与环境科研人员、大气科学方向研究生及Python气象数据处理初学者的WRF/WRF-Chem全流程自动化脚本集聚焦模型前处理地形/土地利用生成、初始场插值、namelist配置与后处理NetCDF解析、时空可视化、特征导出两大核心环节。压缩包共95个文件含62个Python主脚本覆盖ERA5数据下载、TsinghuaLU/VIIRS/Landsat8遥感数据预处理、垂直剖面绘制、泰勒图生成、UCM城市参数构建等、15个备份脚本.zbak、5张示例图表含geogrid、plaindraw、泰勒图等关键输出效果以及字体、许可证、README等配套文件总大小17.17MB。已有37人学习下载。用户可直接复用近百个模块化脚本快速完成从原始遥感数据→WRF输入文件→模拟结果分析→AI建模特征提取的完整链路尤其适配Scikit-learn特征工程需求显著降低WRF-Chem数据处理门槛与重复劳动成本。1. 项目概述从“跑通”到“跑好”的必经之路如果你正在或即将使用WRFWeather Research and Forecasting或WRF-Chem模型进行气象或大气化学模拟那么你大概率会遇到一个比模型本身运行更耗时、更令人头疼的环节数据的前后处理。模型本身像一台精密的发动机但给它喂什么油输入数据以及如何解读它排出的尾气输出数据才是决定你科研或业务工作成败的关键。这就是Python脚本大显身手的地方。我接触WRF系列模型超过十年从最初的手动修改namelist、用NCL画图到后来全面转向Python自动化流程深刻体会到一套高效、可靠的数据处理脚本是多么重要。它不仅仅是节省时间更是保证数据质量、实现可重复性研究、以及进行复杂分析比如化学物种追踪、源解析的基石。很多人把大量精力花在调试模型物理参数上却忽略了前后处理流程的规范化和自动化导致结果难以复现或者因为某个手动步骤的疏漏而前功尽弃。这个项目就是聚焦于用Python构建一套针对WRF/WRF-Chem模型的“预处理”与“后处理”脚本体系。预处理负责将原始再分析数据如FNL、GFS、ERA5或排放清单转换成WRF能“消化”的格式如WPS处理的中间文件或直接准备化学初始边界条件后处理则是将WRF输出的NetCDF文件转换为我们需要的分析结果如特定变量的时空分布图、垂直剖面、站点时间序列、轨迹计算、化学物种收支分析等。我们将深入每个环节不仅告诉你“怎么做”更会解释“为什么这么做”并分享那些在官方文档里找不到的实战经验和避坑指南。2. 核心需求与工具栈解析为什么是Python在开始动手之前我们需要明确核心需求并理解为什么Python是目前处理WRF数据生态中的最优选。2.1 核心需求拆解一个完整的WRF/WRF-Chem模拟流程对数据处理脚本的需求可以归纳为以下几点自动化与批处理模拟往往是多案例、多参数、长时间序列的。手动点击操作不现实脚本必须能自动循环处理不同日期、不同区域、不同排放情景的数据。格式转换与重映射这是预处理的核心。需要将全球经纬度网格的再分析数据通过WPSWRF Preprocessing System插值到我们定义的高分辨率区域网格上。对于WRF-Chem还需要处理排放清单如EDGAR、MEIC将其从原始格式如NetCDF、文本转换为WRF-Chem所需的wrfchemi_*文件网格化、时间解析的排放源。数据提取与子集化WRF输出文件通常包含所有变量、所有层次、所有时间步文件巨大。后处理脚本需要能快速、准确地提取感兴趣的变量如温度、降水、PM2.5、O3、层次如地面、850hPa、500hPa和区域如某个城市群。计算与诊断很多需要的量不是直接输出的需要基于原始变量进行计算。例如计算位温、相当位温、水汽通量、大气稳定度指数如CAPE或者对于WRF-Chem计算气溶胶光学厚度AOD、二次有机气溶胶SOA产量、物种的干湿沉降通量等。可视化与制图生成符合学术出版或业务报告要求的高质量图表包括空间填色图、等值线图、垂直剖面图、时间序列图、散点图等。质量控制与错误处理流程中任何一个环节出错如下载数据不完整、WPS运行失败、磁盘空间不足脚本应能及时发现、记录日志并尽可能优雅地处理或中断避免浪费计算资源。2.2 为什么选择Python及其工具栈早年WRF社区后处理多依赖NCLNCAR Command Language和GrADS。它们功能强大但在代码可维护性、生态系统丰富度、学习曲线和社区活跃度上已逐渐被Python超越。丰富的科学计算库NumPy提供高效的数组计算xarray是处理NetCDF数据的“神器”它引入了带标签的多维数组概念操作NetCDF就像操作Pandas DataFrame一样直观完美替代了NCL的大部分功能pandas处理时间序列和表格数据得心应手。强大的可视化能力Matplotlib是基石高度可定制Cartopy是新一代的地图绘图库替代了古老的Basemap支持更多投影和更现代的地理数据处理对于交互式探索Plotly或HvPlot是很好的选择。专门的WRF处理包wrf-python是UCAR官方维护的Python工具包提供了一系列函数用于从WRF输出文件中提取变量、插值到气压层、计算诊断量如涡度、散度、降水类型等极大简化了后处理。卓越的自动化与集成能力Python可以轻松调用系统命令subprocess模块来运行WPS、real.exe、wrf.exe可以编写爬虫自动下载FNL/GFS数据也可以与任务调度系统如Slurm、PBS结合实现端到端的自动化模拟流程。活跃的社区与可复用性代码易于模块化、函数化方便团队共享和积累。Jupyter Notebook为交互式分析和教学提供了绝佳环境。注意虽然Python是主力但并不意味着完全抛弃WPS。WPS中的geogrid、ungrib、metgrid程序仍然是预处理不可或缺的核心。我们的Python脚本角色是“编排者”和“增强者”负责准备输入、控制流程、处理WPS生成的数据以及进行后续深度分析。3. 预处理脚本详解为模型准备“食材”预处理是模拟的起点目标是为real.exeWRF或real.execonvert_emiss.exeWRF-Chem生成正确的输入文件。我们将流程分解并实现自动化。3.1 驱动数据下载与时间管理WRF需要初始场和边界条件通常来自NCEP FNL、GFS或ECMWF ERA5再分析数据。手动下载费时费力。脚本核心功能自动生成下载列表根据模拟的起止时间start_date,end_date和时效如GFS每6小时一次生成所有需要下载的文件时间列表。调用下载工具使用subprocess调用wget、curl或专门的客户端如cdsapifor ERA5rda-apifor FNL进行批量下载。这里需要处理认证API key和网络错误重试。完整性校验下载后检查文件大小或使用ncdump -h快速检查文件头确保数据完整。import subprocess from datetime import datetime, timedelta def download_gfs_data(start_dt, end_dt, interval_hours6, output_dir./gfs_data): 自动下载GFS数据示例需预先配置好wget或ftp base_url ftp://ftp.ncep.noaa.gov/pub/data/nccf/com/gfs/prod/gfs.{date:%Y%m%d/%H}/ os.makedirs(output_dir, exist_okTrue) current_dt start_dt while current_dt end_dt: # 构造文件名例如 gfs.t00z.pgrb2.0p25.f000 # 实际文件名更复杂需要根据分辨率确定 file_pattern fgfs.t{current_dt.hour:02d}z.pgrb2.0p25.f000 full_url base_url.format(datecurrent_dt) file_pattern cmd [wget, -c, -N, -P, output_dir, full_url] try: subprocess.run(cmd, checkTrue, timeout300) print(f成功下载: {file_pattern}) except subprocess.CalledProcessError as e: print(f下载失败 {file_pattern}: {e}) # 这里可以加入重试逻辑 except subprocess.TimeoutExpired: print(f下载超时 {file_pattern}) current_dt timedelta(hoursinterval_hours)实操心得使用-c断点续传和-N只下载比本地新的文件是下载大文件时的好习惯避免重复下载。ERA5数据通过ECMWF的CDS API下载是最佳途径需要先注册获取API key。使用cdsapi库时注意请求的数据量过大可能导致排队脚本中需要加入等待和状态查询逻辑。本地缓存对于公共数据建议在服务器或工作站建立本地缓存目录所有项目共享避免重复下载占用带宽和磁盘。3.2 WPS流程自动化封装运行WPSgeogrid - ungrib - metgrid通常涉及多次修改namelist.wps和执行命令。我们可以用Python脚本统一管理。脚本核心功能动态生成namelist.wps将域设置、时间范围、路径参数等写入一个模板字典然后用f90nml库专门读写Fortran namelist的Python库或简单的字符串替换生成最终的namelist.wps文件。这比手动编辑更不易出错尤其在做参数敏感性试验时。顺序执行与错误监控用subprocess依次运行geogrid.exe、ungrib.exe、metgrid.exe。关键是要检查每个步骤的日志文件如geogrid.log或退出状态码一旦出错就立即停止并报错而不是继续执行下去产生无用的中间文件。清理中间文件ungrib会生成大量的FILE:*文件在确认metgrid成功后脚本可以自动清理这些临时文件释放磁盘空间。import f90nml import subprocess import shutil def run_wps(domain_config, time_config, wps_dir./WPS): 自动化运行WPS流程 domain_config: 字典包含max_dom, dx, dy, parent_id等域信息 time_config: 字典包含start/end_date, interval_seconds等时间信息 # 1. 准备namelist.wps namelist_template { share: { wrf_core: ARW, max_dom: domain_config[max_dom], start_date: time_config[start_date], # 格式: 2023-01-01_00:00:00 end_date: time_config[end_date], interval_seconds: time_config[interval_seconds], }, geogrid: { ... }, # 地理数据路径、地图参数等 ungrib: { ... }, metgrid: { ... }, } namelist_path os.path.join(wps_dir, namelist.wps) f90nml.write(namelist_template, namelist_path) # 2. 链接Vtable (根据数据源选择) vtable_src fVtable.{data_source} # 如 Vtable.GFS vtable_dst os.path.join(wps_dir, Vtable) if os.path.exists(vtable_src): if os.path.lexists(vtable_dst): os.unlink(vtable_dst) os.symlink(vtable_src, vtable_dst) else: raise FileNotFoundError(fVtable文件 {vtable_src} 不存在) # 3. 运行geogrid os.chdir(wps_dir) try: subprocess.run([./geogrid.exe], checkTrue, stdoutsubprocess.PIPE, stderrsubprocess.PIPE) print(geogrid.exe 执行成功) except subprocess.CalledProcessError as e: print(fgeogrid.exe 执行失败: {e}) with open(geogrid.log, r) as f: print(f.read()) return False # 4. 运行ungrib # ... 类似逻辑先链接GRIB数据文件 # 5. 运行metgrid # ... # 6. 可选清理FILE:*文件 # if keep_intermediate is False: # for f in glob.glob(FILE:*): # os.remove(f) return True注意事项路径问题geogrid.exe需要访问静态地理数据GEOG目录确保路径在namelist中设置正确并且该目录有足够权限和空间。Vtable匹配不同的驱动数据GFS, FNL, ERA5需要不同的Vtable文件链接错误会导致ungrib读不出变量。内存与磁盘高分辨率、多层次的模拟metgrid可能会生成巨大的中间文件确保/tmp或工作目录有足够空间TB级别。3.3 WRF-Chem排放清单处理专项这是WRF-Chem预处理中最具挑战性的部分之一。排放清单通常以年度全球网格如0.1°x0.1°或区域列表形式提供需要将其处理成WRF-Chem能识别的、时间分辨率更高如逐小时、空间上与模拟域匹配的NetCDF文件。核心步骤与脚本实现清单读取与合并使用xarray打开NetCDF格式的排放清单如EDGAR或用pandas读取CSV文本清单。可能需要合并人为源、生物源、海盐、沙尘等不同来源。时空重映射空间使用xesmf库一个基于ESMF的Python重映射工具或scipy插值将清单从原始网格重映射到WRF模拟域通过geo_em.d0x.nc文件定义的不规则网格。这比简单的最近邻插值更准确特别是对化学物种。时间排放具有日变化、周变化和季节变化。脚本需要根据预设的时间变化因子如交通排放的早晚高峰、工业排放的昼夜差异将年均或月均排放量分解到每个模拟小时。这通常通过一个时间剖面文件来实现。物种映射与分割排放清单中的物种如NOx, SO2, NMVOC需要映射到WRF-Chem化学机制如RADM2, MOZART, SAPRC99中的具体物种。同时需要将VOC按照化学机制要求的比例分割成具体的模型物种如PAR, OLE, TOL等。生成wrfchemi文件调用convert_emiss.exe需要编译进WRF-Chem或使用anthro_emis等预处理程序但更灵活的方式是直接用xarray和netCDF4库按照wrfchemi的格式要求从头创建NetCDF文件。这要求你深刻理解该文件的维度Time, emissions_zdim, south_north, west_east和变量属性。import xarray as xr import xesmf as xe def process_emission_to_wrfchem(orig_emis_nc, geo_em_file, chem_mechanismRADM2, output_prefixwrfchemi): 将原始排放清单处理为WRF-Chem输入文件简化示例 # 1. 读取原始排放数据和WRF地理网格文件 ds_emis xr.open_dataset(orig_emis_nc) ds_geo xr.open_dataset(geo_em_file) # 2. 创建从原始网格到WRF网格的重映射器 (使用双线性插值) # 注意排放重映射需谨慎质量守恒是关键。xesmf提供守恒插值选项。 regridder xe.Regridder(ds_emis, ds_geo, bilinear, reuse_weightsTrue) # 3. 对每种排放物种进行重映射 emis_vars [NOx, SO2, CO, VOC] # 示例变量名 emis_regridded {} for var in emis_vars: if var in ds_emis: emis_regridded[var] regridder(ds_emis[var]) # 4. 应用时间变化因子此处简化实际需根据小时、星期几计算 # 假设我们有一个24小时的时间变化因子数组 time_profile [24] # 需要将年均排放量根据模拟的每个小时进行分配 # 这是一个复杂的步骤涉及日期时间计算和乘法运算 # 5. 物种映射与分割 (以VOC为例) # 假设 chem_mechanism_mapping 是一个字典定义了原始VOC到机制物种的分配系数 # 例如{PAR: 0.4, OLE: 0.3, TOL: 0.1, ...} # emis_regridded[PAR] emis_regridded[VOC] * 0.4 # ... # 6. 创建符合wrfchemi格式的Dataset # 需要严格按照WRF-Chem要求的维度、变量名、属性和单位来构建 # 此处省略详细的构建过程它涉及创建时间维度、垂直层次通常为1层等 # 7. 写入NetCDF文件 # ds_wrfchemi.to_netcdf(f{output_prefix}_00z_d01) print(排放处理流程完成示例)避坑技巧质量守恒检查重映射前后计算全球或区域的总排放量确保误差在可接受范围内1%。xesmf的conservative插值方法可以更好地保持总量。垂直分配大部分清单只提供地面排放。对于高空点源如电厂烟囱需要根据烟囱参数高度、温度、速度进行垂直分配这通常需要额外的处理脚本或使用PREP-CHEM-SRC等工具。背景化学场除了排放WRF-Chem还需要化学初始边界条件如wrfinput中的化学物种浓度。对于全球化学模拟如MOZART提供的数据处理方式与气象场类似但变量更多。确保ungrib使用的Vtable能识别化学变量或者使用专门的工具如mozbc来处理。4. 后处理脚本详解从海量输出中提炼“洞察”WRF模拟完成后输出是庞大的NetCDF文件集wrfout_d0*_*。后处理的目标是将其转化为有科学或业务价值的图形和数据集。4.1 高效数据提取与子集化直接使用NetCDF4库或xarray打开所有wrfout文件进行全局操作极易导致内存溢出。必须采用惰性加载和分块处理的策略。脚本策略使用xarray的延迟计算xarray.open_mfdataset函数可以同时打开多个文件并创建一个逻辑上统一的数据集而不会立即将数据读入内存。结合chunks参数指定分块大小可以与Dask库集成实现并行计算。选择性提取在打开数据集时或之后立即使用.sel()选择和.isel()索引选择方法提取需要的Time、bottom_top或bottom_top_stag、south_north、west_east维度子集。变量计算后再持久化对于需要计算的诊断变量利用xarray的向量化运算在“虚拟”数据集上进行最后对结果调用.compute()或.persist()才真正触发计算并加载到内存。对于最终需要反复使用的中间结果可以保存到新的、更小的NetCDF文件中。import xarray as xr import numpy as np def extract_and_calc_pm25(wrfout_files, domain1, save_path./processed/pm25.nc): 从一系列wrfout文件中提取并计算PM2.5浓度适用于包含气溶胶物种的机制 假设使用MOSAIC或GOCART气溶胶方案 # 1. 延迟打开多个文件 # 使用 chunks{Time: 10, south_north: 100, west_east: 100} 进行分块 ds xr.open_mfdataset(wrfout_files, chunks{Time: 24}, combineby_coords, parallelTrue) # 2. 选择特定域 (假设变量有‘XTIME’等维度标识) # 更常见的是直接打开对应domain的文件如 wrfout_d01_* # 这里我们假设ds已经只包含一个域的数据 # 3. 计算PM2.5 (以GOCART机制为例PM2.5近似为部分物种之和) # 注意不同化学机制PM2.5的算法不同此处仅为示例。 # 例如可能包括硫酸盐(so4)、硝酸盐(no3)、铵盐(nh4)、有机碳(oc)、黑碳(bc)、海盐(ss)、沙尘(du)的细模态部分 # 需要根据实际输出变量名调整 species_list [so4_a01, no3_a01, nh4_a01, oc_a01, bc_a01] # 细模态气溶胶 pm25 None for species in species_list: if species in ds: if pm25 is None: pm25 ds[species] else: pm25 pm25 ds[species] if pm25 is not None: # 添加单位等属性 pm25.attrs[long_name] Fine Particulate Matter (PM2.5) pm25.attrs[units] ug m-3 # 4. 将pm25作为一个新变量放回数据集或创建新数据集 ds_pm25 pm25.to_dataset(namePM2.5) # 5. 保存到新文件触发实际计算和写入 ds_pm25.to_netcdf(save_path) print(fPM2.5数据已保存至 {save_path}) return ds_pm25 else: print(未找到计算PM2.5所需的物种变量。) return None注意事项变量名差异WRF输出变量名与所选物理参数化方案、化学机制紧密相关。务必先用ncdump -h或xarray查看变量列表。wrf-python的getvar函数可以屏蔽部分差异通过通用名称如‘temp’ ‘uvmet’获取变量。内存管理对于超大型数据集如高分辨率、长时间模拟即使分块也可能在计算聚合操作如时间平均、空间平均时内存暴涨。考虑分步处理先按时间片处理保存中间结果再对中间结果进行聚合。时间处理WRF输出时间可能是模拟分钟数或日期字符串。使用wrf-python的getvar函数配合time关键字或使用xarray的cftime支持可以将其转换为Python的datetime对象方便后续分析。4.2 诊断量与可视化计算除了直接输出变量许多重要物理量需要计算。wrf-python库封装了大部分常用诊断量的计算。常用诊断量计算示例import wrf import numpy as np def calculate_meteorological_diagnostics(ds): 计算一些常用的气象诊断量 ds: 通过wrf-python或xarray打开的wrfout数据集 # 注意wrf.getvar 需要原始文件路径或已加载的netCDF4 Dataset对象 # 这里假设ds是netCDF4.Dataset对象 # 如果ds是xarray Dataset可能需要from netCDF4 import Dataset; ncfile Dataset(wrfout_file) # 获取经纬度坐标用于绘图 lats, lons wrf.latlon_coords(ds) # 获取地面气压、温度、露点温度等 slp wrf.getvar(ds, slp) # 海平面气压 t2 wrf.getvar(ds, T2) # 2米温度 td2 wrf.getvar(ds, td2) # 2米露点温度 u10, v10 wrf.getvar(ds, uvmet10) # 10米风地图投影坐标 # 计算相对湿度 rh2 wrf.getvar(ds, rh2) # 2米相对湿度 # 计算降水量累计降水转瞬时或逐小时 rainc wrf.getvar(ds, RAINC) # 对流降水 rainnc wrf.getvar(ds, RAINNC) # 非对流降水 total_rain rainc rainnc # 总累计降水 # 计算逐小时降水需要时间维度 # hour_rain total_rain.diff(dimTime) # 注意单位换算 # 计算垂直剖面例如沿一条线的位温剖面 # 假设我们想画一条从点(start_lat, start_lon)到(end_lat, end_lon)的垂直剖面 start_point (35.0, 115.0) # (lat, lon) end_point (40.0, 120.0) # 提取剖面线上的数据 # 注意wrf.vertcross 需要等压面或高度面上的变量需要先插值 p wrf.getvar(ds, pressure) # 获取气压场 tk wrf.getvar(ds, tk) # 温度 (K) # 将温度插值到气压面上 # 这里只是示意实际调用较复杂 # theta_cross wrf.vertcross(theta, p, wrf_line, latlonTrue) diagnostics { slp: slp, t2: t2, td2: td2, u10: u10, v10: v10, rh2: rh2, total_rain: total_rain, lats: lats, lons: lons } return diagnostics可视化脚本核心 计算完成后使用Cartopy和Matplotlib绘图。关键是要正确处理WRF的地图投影信息。import cartopy.crs as ccrs import cartopy.feature as cfeature import matplotlib.pyplot as plt from wrf import get_cartopy, cartopy_xlim, cartopy_ylim def plot_surface_map(wrf_file, var_name, time_idx0, level0, domain1): 绘制WRF输出变量的地面填色图 # 从文件获取数据和地图投影信息 ncfile Dataset(wrf_file) var_data wrf.getvar(ncfile, var_name, timeidxtime_idx) # 获取地图投影 cart_proj get_cartopy(var_data) # 创建图形 fig plt.figure(figsize(12, 8)) ax plt.axes(projectioncart_proj) # 添加地理特征 ax.add_feature(cfeature.COASTLINE.with_scale(50m), linewidth0.8) ax.add_feature(cfeature.BORDERS.with_scale(50m), linestyle:, linewidth0.5) ax.add_feature(cfeature.STATES.with_scale(50m), linewidth0.3) # 添加河流、湖泊等... # 设置图形范围 ax.set_xlim(cartopy_xlim(var_data)) ax.set_ylim(cartopy_ylim(var_data)) # 绘制填色图 # 需要将数据坐标转换为地图投影坐标 lats, lons wrf.latlon_coords(var_data) # 使用pcolormesh或contourf # 注意如果var_data是三维带bottom_top需要切片 var_data[level, :, :] if len(var_data.shape) 3: plot_data var_data[level, :, :] else: plot_data var_data # 这里使用pcolormesh对于不规则WRF网格需使用经纬度坐标 mesh ax.pcolormesh(lons, lats, plot_data, transformccrs.PlateCarree(), shadingauto, cmapviridis) # 添加色标 plt.colorbar(mesh, axax, orientationhorizontal, pad0.05, shrink0.8) # 添加标题 plt.title(f{var_data.description} ({var_data.units}) at Time Index {time_idx}) plt.tight_layout() plt.savefig(f{var_name}_d0{domain}_t{time_idx}.png, dpi300) plt.close(fig) ncfile.close()实操心得投影是关键WRF使用自己的地图投影通常在namelist中定义如Lambert Conformal, Mercator等。wrf-python的get_cartopy函数能自动提取这个投影信息确保绘制的图形在空间上是正确的。直接使用PlateCarree投影绘制会导致严重变形。性能优化绘制大量时间步的动画或批量出图时避免在循环中重复创建地图投影和地理特征可以将其移到循环外。使用matplotlib的FuncAnimation或批量生成图片后用ffmpeg合成视频。图形美化Cartopy提供了丰富的地图元素。合理使用cfeature添加海岸线、国界、河流、湖泊并使用gridliner添加经纬度网格和标签能让图更专业。5. 实战中的常见问题与排查技巧即使脚本写得再完善在实际运行中也会遇到各种问题。以下是一些典型问题及其解决思路。5.1 预处理阶段问题问题现象可能原因排查步骤与解决方案ungrib.exe运行失败提示“Not a GRIB file”或“No fields found”1. GRIB文件损坏或下载不完整。2. 链接的Vtable与GRIB数据源不匹配。3. GRIB数据版本如GFS 0.25度 vs 0.5度与Vtable预期不符。1. 用grib_dump或wgrib2检查GRIB文件头。2. 核对数据源GFS/FNL/ERA5并链接正确的Vtable。3. 查看WPS日志ungrib.log看具体是哪个变量读不出。有时需要手动修改Vtable以适配新数据格式。metgrid.exe运行缓慢或内存溢出1. 模拟区域太大或分辨率太高。2. 垂直层数过多。3. 输入的气象场时间频率太高。1. 检查namelist.wps中的e_we/e_sn和dx/dy。2. 减少垂直层数需与WRF的namelist.input一致。3. 增加interval_seconds降低驱动数据的时间频率如从3小时降到6小时。4. 在运行metgrid时尝试增加进程数如果支持MPI。生成的wrfinput文件中某些变量全是缺省值1.metgrid插值过程中某些变量在边界处丢失。2. 静态数据如土地利用类型、土壤类型在geogrid步骤未正确读取。1. 检查metgrid.log是否有警告。可以尝试在namelist.wps的metgrid部分调整interp_option如从average_gcell改为nearest_neighbor。2. 检查GEOG目录路径及权限确保所需的地理数据文件存在。5.2 后处理阶段问题问题现象可能原因排查步骤与解决方案用xarray打开wrfout文件报错或变量找不到1. 文件损坏。2. 变量名因物理方案不同而异。3. 文件格式问题如NetCDF4 classic vs 64-bit offset。1. 用ncdump -h快速查看文件结构和变量列表这是最可靠的方法。2. 使用wrf-python的getvar函数它提供了统一的接口来获取物理量如‘temp’ ‘uvmet’内部处理了变量名差异。3. 确保安装的netCDF4库支持对应的NetCDF版本。计算诊断量如涡度、散度结果异常NaN或极大值1. 数据在网格边界或地形附近存在无效值。2.wrf-python的getvar函数参数使用不当如未指定投影坐标。3. 输入变量单位不匹配。1. 使用xarray的.where()或numpy的np.isnan过滤无效值后再计算。2. 仔细阅读wrf-python官方文档确认函数所需的输入数据格式是原始变量还是getvar提取后的。例如计算相对涡度需要风场在地图投影坐标上。3. 打印输入变量的单位和值范围进行验证。绘图时图形扭曲或坐标错误1. 使用了错误的地图投影。2. 数据经纬度坐标与绘图投影不匹配。3. Cartopy版本与WRF-Python兼容性问题。1.始终坚持使用wrf.get_cartopy()获取投影并用transformccrs.PlateCarree()参数将数据lons, lats转换到地理坐标进行绘制。2. 使用wrf.cartopy_xlim/ylim设置图形范围而不是手动指定。3. 检查wrf-python和Cartopy的版本尽量使用较新且兼容的版本组合。处理大量文件时内存不足1. 一次性加载所有数据。2. 中间变量未及时释放。1.使用xarray.open_mfdataset并设置chunks参数启用Dask进行惰性加载和分块计算。2. 对于循环处理显式关闭不再需要的文件句柄ds.close()。3. 将大任务分解为多个小任务分步处理并保存中间结果。5.3 脚本调试与优化心得日志是生命线在脚本的关键步骤开始、结束、出错点添加详细的日志记录包括时间、操作、关键参数和状态。使用Python的logging模块可以方便地控制日志级别DEBUG, INFO, WARNING, ERROR并输出到文件这对于在后台长时间运行的批处理任务至关重要。增量式开发与测试不要试图一次性写完处理整个流程的巨型脚本。应该模块化开发每个函数只负责一个明确的任务如下载、运行WPS、计算PM2.5、绘图。然后针对一个小的、代表性的案例如单日、低分辨率模拟进行完整流程测试确保每个模块都工作正常。参数配置文件将模拟的域设置、时间范围、路径、化学机制选项等所有可配置参数写进一个单独的配置文件如config.yaml或config.ini。主脚本读取这个配置文件。这样当需要更换模拟案例时只需修改配置文件而无需改动脚本代码大大提高了复用性和可维护性。利用并行加速对于可并行的任务如处理不同时间段的文件、绘制多幅图使用multiprocessing或concurrent.futures模块实现多进程并行。xarray结合Dask可以自动对分块数据进行并行计算。注意平衡任务粒度和进程数避免进程间通信开销过大。6. 从脚本到流程构建自动化模拟系统当预处理和后处理脚本都稳定可靠后我们可以将它们整合构建一个端到端的自动化模拟与后处理系统。这个系统通常由一个主控脚本驱动它协调以下任务参数解析与配置加载读取用户提供的案例配置文件。数据准备调用下载脚本获取驱动数据调用排放处理脚本生成化学输入。WPS/WRF执行依次运行geogrid,ungrib,metgrid,real.exe,wrf.exe可能还有convert_emiss.exe。每一步都进行状态检查失败则发送通知如邮件并中止。后处理流水线WRF运行成功后自动触发后处理脚本进行数据提取、诊断量计算和批量绘图。结果打包与归档将关键的输入文件namelist、日志、核心输出数据和图表打包转移到归档存储并清理临时工作目录以释放空间。这样的系统可以部署在高性能计算集群上通过作业调度系统如Slurm提交实现无人值守的批量模拟实验极大提升科研和生产效率。我个人在构建这类系统时最后通常会写一个简单的Makefile或使用Snakemake/Nextflow这样的工作流管理工具来定义任务之间的依赖关系这比纯Python脚本在管理复杂流水线时更加清晰和强大。例如可以定义“绘图”任务依赖于“计算PM2.5”任务而“计算PM2.5”又依赖于“WRF运行完成”。当某个中间文件更新时工作流引擎可以自动识别并重新运行下游任务。无论采用哪种方式核心思想是一致的将重复、易错的手工操作固化为可靠、可复用的代码把宝贵的精力留给真正的科学问题分析和模型物理机制的探索。这套基于Python的WRF数据处理脚本就是你从模型“使用者”进阶为模型“驾驭者”的关键工具。本文还有配套的精品资源点击获取
