【STILT模型实操第2期】运行 STILT 模型提供 WRF 示例
目录1. 项目初始化与数据准备1.1 初始化 STILT 项目1.2 下载示例数据1.3 准备转换工具2. 核心操作步骤步骤 1批量转换 WRF 数据步骤 2配置 STILT 运行脚本 (r/run_stilt.r)步骤 3运行 STILT 模型步骤 4检查输出结果可视化轨迹文件 (..._traj.rds)可视化足迹文件 (..._foot.nc)3. 报错总结错误1R库包缺失raster错误2object simulation_id not found错误3object met_file_tres not found参考WRF-STILT模型用户教程可参考另一博客-【STILT模型第4期】WRF-STILT模型用户教程。STILT模型的安装可参考另一博客-【STILT模型实操第1期】STILT模型安装环境配置及安装本博客结合 GitHUb提供的 WRF 示例完整实现 STILT 模型模拟。1. 项目初始化与数据准备在开始之前我们需要创建一个新的 STILT 项目并下载教程所需的示例数据。1.1 初始化 STILT 项目在终端中使用 R 语言的uataq包初始化一个名为wrf-tutorial的新项目并进入该目录Rscript-euataq::stilt_init(wrf-tutorial)cdwrf-tutorial也可按照个人习惯修改文件名称stilt_init(1_stilt_test_project)cd/data/wanzhou/STILT/1_stilt_test_project/1.2 下载示例数据在项目的根目录下克隆 GitHub 上的教程数据仓库包含 WRF 示例输出文件gitclone--depth1https://github.com/uataq/stilt-tutorials# 清空文件夹rm-r/data/wanzhou/STILT/1_stilt_test_project/stilt-tutorials/1.3 准备转换工具将之前提到的arw2arl转换工具复制到教程工作目录中cpexe/arw2arl stilt-tutorials/03-wrf/完成以上步骤后项目目录结构应该如文档中所示包含了exe/、r/以及存放了 WRF 示例数据的stilt-tutorials/03-wrf/等文件夹。PROJECTbin/... docs/... exe/ arw2arl... r/ run_stilt.r... stilt-tutorials/ 03-wrf/ wrfout/ wrfout_*.nc... arw2arl arw2arl_batch.sh... test/...2. 核心操作步骤步骤 1批量转换 WRF 数据切换到 WRF 教程目录并运行提供的 Shell 批处理脚本。sourceactivate stilt_env(cd stilt-tutorials/03-wrf./arw2arl_batch.sh)注意请确保配置文件WRFDATA.CFG与脚本在同一目录下。结果该脚本运行完毕后会在stilt-tutorials/03-wrf/arlout目录下生成按天分组的 ARL 格式气象数据文件。步骤 2配置 STILT 运行脚本 (r/run_stilt.r)成功转换数据后需要编辑 STILT 的主配置文件r/run_stilt.r。# 拷贝备份cp/data/wanzhou/STILT/1_stilt_test_project/r/run_stilt.r/data/wanzhou/STILT/1_stilt_test_project/r/run_stilt_template.r由于本次模拟的中心位于美国犹他州盐湖城 (Salt Lake City)需要根据提供的代码块修改以下关键参数# Simulation timing, yyyy-mm-dd HH:MM:SS (UTC)t_start-2015-09-06 00:00:00t_end-2015-09-06 01:00:00run_times- seq(fromas.POSIXct(t_start, tzUTC), toas.POSIXct(t_end, tzUTC), byhour)# Receptor location(s)lati-40.740long--111.860zagl-5# Footprint grid settings, must set at least xmn, xmx, ymn, ymx belowhnf_plume- T projection-projlonglatsmooth_factor-1time_integrate- F xmn--114.0xmx--110.99ymn-39.0ymx-42.0xres-0.01yres- xres# Meteorological data inputmet_path- file.path(stilt_wd,stilt-tutorials,03-wrf,arlout)met_file_format-%Y-%m-%dmet_subgrid_buffer-0.1met_subgrid_enable- F met_subgrid_levels- NA n_met_min-1# Model controln_hours--6numpar-500rm_dat- T run_foot- T run_trajec- Ttimeout-3600varsiwant- c(time,indx,long,lati,zagl,foot,mlht,dens,samt,sigw,tlgr,dmas)模拟时间 (Simulation timing)设置为2015-09-06 00:00:00到01:00:00。接收点位置 (Receptor location)纬度40.740经度-111.860离地高度5米。足迹网格设置 (Footprint grid settings)设定经纬度边界 (xmn,xmx,ymn,ymx) 及分辨率 (xres,yres)。气象数据路径 (Meteorological data input)met_path - file.path(stilt_wd, stilt-tutorials, 03-wrf, arlout) met_file_format - %Y-%m-%d模型控制 (Model control)例如后向追踪时间n_hours - -6向后追踪 6 小时粒子数numpar - 500等。(注请直接参考您提供的代码块将这些参数替换到run_stilt.r文件中。)步骤 3运行 STILT 模型保存修改后的r/run_stilt.r脚本。然后回到终端在项目根目录下执行该脚本cd/data/wanzhou/STILT/1_stilt_test_project/ Rscript r/run_stilt.r步骤 4检查输出结果如果模型运行成功STILT 生成的轨迹文件 (trajectory) 和足迹文件 (footprint) 将会被保存在项目根目录下的./out文件夹中。轨迹文件 (…_traj.rds) —— “粒子去哪了”本质记录了数千个拉格朗日粒子在回溯时间内如 24 小时或 72 小时每一时刻的三维坐标经纬度、高度、混合层高度、气温等物理量。作用让你看到空气团是从哪个方向飘过来的经过了哪些地区。足迹文件 (…_foot.nc) —— “地表贡献了多少”本质这是轨迹数据经过空间网格化处理后的结果通常是 NetCDF 格式。单位通常是p p m / ( μ m o l / ( m 2 ⋅ s ) ) ppm / (\mu mol / (m^2 \cdot s))ppm/(μmol/(m2⋅s))。作用它定义了地表排放对观测点浓度的敏感性。如果某个格点的足迹值很高说明那个地方的排放对你观测到的浓度贡献巨大。可以进入该目录查看最终的模拟结果。可视化轨迹文件 (…_traj.rds)cd/data/wanzhou/STILT/1_stilt_test_project/out/by-id/201509060000_-111.86_40.74_5/sourceactivate stilt_env绘图相关库包安装如下install.packages(c(ggplot2, terra, ncdf4, maps))R 绘图代码如下# # 1. 绘制并保存轨迹 (Trajectory) 图片# # 读取 RDS 文件traj_data_raw-readRDS(traj_file)# 【新增处理】检查并转换数据格式if(is.data.frame(traj_data_raw)){traj_data-traj_data_raw}elseif(is.list(traj_data_raw)){# 如果是 list通常说明它是多个轨迹数据框的集合将其按行合并 (rbind)cat(检测到轨迹数据为 list 格式正在尝试合并为 data.frame...\n)# 使用 do.call 将 list 中的所有元素合并为一个 data.frametraj_data-do.call(rbind,traj_data_raw)# 注意如果你安装了 dplyr 包使用 dplyr::bind_rows(traj_data_raw) 会更安全且高效}else{stop(未知的轨迹数据格式请运行 str(readRDS(traj_file)) 检查数据结构。)}# 提取地图背景states_map-map_data(state)# 绘制轨迹图p_traj-ggplot()# 绘制地图背景geom_polygon(datastates_map,aes(xlong,ylat,groupgroup),fillwhite,colorgray70)# 绘制粒子轨迹geom_point(datatraj_data,aes(xlong,ylati,colorzagl),size0.5,alpha0.6)scale_color_viridis_c(nameHeight (m))coord_quickmap(xlimrange(traj_data$long,na.rmTRUE)c(-2,2),ylimrange(traj_data$lati,na.rmTRUE)c(-2,2))theme_minimal()labs(titleSTILT Particle Trajectories,xLongitude,yLatitude)# 保存轨迹图片ggsave(filenamepaste0(output_dir,trajectory_plot.png),plotp_traj,width8,height6,dpi300)cat(轨迹图已保存至:,paste0(output_dir,trajectory_plot.png\n))可视化足迹文件 (…_foot.nc)使用 ncview 进行结果可视化R 绘图如下R 绘图代码如下# # 2. 绘制并保存足迹 (Footprint) 图片# # 读取 NetCDF 文件# STILT 的 foot.nc 可能包含多个时间步长 (layers)这里默认读取并对所有层求时间积分求和foot_raster-rast(foot_file)# 如果有多个时间维度通常将它们相加得到总的 footprintif(nlyr(foot_raster)1){foot_total-sum(foot_raster,na.rmTRUE)}else{foot_total-foot_raster}# 将 0 值替换为 NA 以便在图中透明显示背景foot_total[foot_total0]-NA# 转换为 data.frame 以便使用 ggplot2 绘图foot_df-as.data.frame(foot_total,xyTRUE)colnames(foot_df)[3]-footprintp_foot-ggplot()geom_polygon(datastates_map,aes(xlong,ylat,groupgroup),fillwhite,colorgray70)geom_raster(datafoot_df,aes(xx,yy,fillfootprint))scale_fill_viridis_c(optionplasma,translog10,na.valuetransparent,nameFootprint\n(log10))coord_quickmap(xlimrange(foot_df$x,na.rmTRUE),ylimrange(foot_df$y,na.rmTRUE))theme_minimal()labs(titleSTILT Surface Influence Footprint,xLongitude,yLatitude)# 保存足迹图片ggsave(filenamepaste0(output_dir,footprint_plot.png),plotp_foot,width8,height6,dpi300)cat(足迹图已保存至:,paste0(output_dir,footprint_plot.png\n))3. 报错总结错误1R库包缺失raster运行 r/run_stilt.r R命令报错如下在运行 run_stilt.r 时程序发现缺少 raster 包并尝试自动安装它。但在安装 raster 的依赖包 terra 时系统卡在了底层地理空间库的缺失上。在 Linux 终端保证在 stilt_env 环境下运行以下命令使用 conda 安装所需库包condainstall-cconda-forge r-raster condainstall-cconda-forge r-terra r-raster r-rgdal错误2object ‘simulation_id’ not found错误如下修改 run_stilt.r 文件增加 simulation_id - NA 配置错误3object ‘met_file_tres’ not found错误如下修改 run_stilt.r 文件增加 ‘met_file_tres’ 配置参考1、GitHub-WRF-STILT tutorial