<think>我们首先需要明确用户的问题：FLEXPART-WRF和FLEXPART输出文件是否存在差异，以及如何对FLEXPART输出数据进行绘图处理。 根据之前提供的知识，FLEXPART-WRF是FLEXPART与WRF模式耦合的版本，其输出文件格式与传统的FLEXPART（例如基于ECMWF气象数据的版本）在数据结构上基本一致，都是NetCDF格式，包含浓度场、粒子轨迹和释放源信息。但是，由于驱动气象场来源不同（WRF提供的高分辨率气象场），输出数据在空间分辨率和时间分辨率上可能有差异，且输出变量名称和维度可能略有不同。不过，现代FLEXPART版本（包括FLEXPART-WRF）都采用了统一的NetCDF输出格式，因此处理方式可以通用。 接下来，我们重点回答两个部分： 1. FLEXPART-WRF与FLEXPART输出文件的差异（如果有的话） 2. 对FLEXPART输出数据进行绘图处理的方法（使用Python） ### 一、FLEXPART-WRF与FLEXPART输出文件的差异 1. **气象驱动源差异**： - FLEXPART-WRF使用WRF模式输出的高分辨率气象场，而传统FLEXPART可能使用全球模式（如ECMWF）的气象场。这导致输出数据在空间和时间分辨率上可能不同，但文件结构相同。 2. **输出文件结构**： - 两者都输出NetCDF格式文件，包含三类文件：浓度文件（`grid_conc_*.nc`）、粒子轨迹文件（`partposit_*.nc`）和释放源文件（`releases_*.nc`）[^2]。 - 变量和维度命名基本一致，例如浓度变量为`spec001_mr`（物质1的混合比），维度为`time`, `height`, `lat`, `lon`。 3. **可能的差异点**： - **垂直坐标**：WRF模式通常使用地形跟随坐标，而全球模式使用气压坐标，但FLEXPART输出时会将浓度插值到固定的高度层上，因此输出文件中的高度层是相同的（以米为单位的地面以上高度）[^2]。 - **地图投影**：WRF模式通常使用某种地图投影（如Lambert投影），但FLEXPART输出时统一转换为经纬度坐标（WGS84），所以用户最终得到的NetCDF文件中的经纬度都是标准的，无需担心投影问题。 因此，**处理FLEXPART-WRF和传统FLEXPART输出的NetCDF文件可以使用相同的Python方法**。 ### 二、FLEXPART输出数据绘图处理（Python） 使用`xarray`和`cartopy`库进行数据处理和绘图。 #### 1. 浓度场绘图（二维平面） 以绘制地面浓度为例： ```python import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs import numpy as np # 打开浓度文件 ds = xr.open_dataset('grid_conc_001.nc') # 将时间戳转换为易读格式（如果时间是以秒为单位的Unix时间戳） ds['time'] = xr.DataArray(pd.to_datetime(ds.time.values, unit='s'), dims='time') # 选择第一个时间步和地面层（高度=0） # 注意：如果高度层有多个，可能需要选择索引或实际高度值 # 例如：height=0 或 height=0.0，或者通过索引 ds.isel(height=0) conc_surface = ds['spec001_mr'].sel(time=ds.time[0], height=0, method='nearest') # 创建地图 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) # 绘制浓度填色图 contour = conc_surface.plot.contourf( ax=ax, transform=ccrs.PlateCarree(), cmap='viridis', levels=np.logspace(-10, -5, 20), # 根据浓度范围调整 add_colorbar=True, cbar_kwargs={'label': 'Concentration (kg/kg)'} ) # 添加海岸线、经纬度网格 ax.coastlines() ax.gridlines(draw_labels=True, linewidth=0.5, color='gray', alpha=0.5, linestyle='--') # 设置标题 plt.title(f"Surface Concentration at {str(ds.time[0].values)[:19]}") # 保存图片 plt.savefig('surface_concentration.png', dpi=300, bbox_inches='tight') ``` #### 2. 垂直剖面图 绘制经度或纬度垂直剖面： ```python # 选择某一经度的垂直剖面（例如经度=116.4） lon_slice = 116.4 conc_cross = ds['spec001_mr'].sel(time=ds.time[0], lon=lon_slice, method='nearest') # 创建二维网格（高度和纬度） X, Y = np.meshgrid(conc_cross.lat, conc_cross.height) # 绘图 plt.figure(figsize=(10, 6)) plt.contourf(X, Y, conc_cross.T, levels=100, cmap='viridis') # 转置使维度匹配 plt.colorbar(label='Concentration (kg/kg)') plt.xlabel('Latitude') plt.ylabel('Height (m)') plt.title(f"Vertical Cross-section at Lon={lon_slice}°E") plt.savefig('vertical_cross_section.png') ``` #### 3. 粒子轨迹动态图 绘制粒子轨迹随时间的变化： ```python import matplotlib.animation as animation # 打开粒子轨迹文件 ds_part = xr.open_dataset('partposit_001.nc') # 转换为DataFrame（注意：粒子轨迹文件通常很大，谨慎操作） df = ds_part.to_dataframe().reset_index() # 选择部分粒子（例如前100个粒子）以减少数据量 particle_ids = df.particle_id.unique()[:100] df_sub = df[df.particle_id.isin(particle_ids)] # 创建动画 fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection=ccrs.PlateCarree()) ax.coastlines() # 初始化散点图 scatter = ax.scatter([], [], s=1, color='red', transform=ccrs.PlateCarree()) # 更新函数 def update(frame): # 选择当前时间步（frame代表时间步索引） time_step = df_sub.time.unique()[frame] df_frame = df_sub[df_sub.time == time_step] scatter.set_offsets(np.c_[df_frame.longitude, df_frame.latitude]) ax.set_title(f"Particle Trajectories at Time: {pd.to_datetime(time_step, unit='s')}") return scatter, # 创建动画 ani = animation.FuncAnimation(fig, update, frames=len(df_sub.time.unique()), interval=100, blit=True) # 保存动画 ani.save('particle_trajectories.gif', writer='pillow', fps=10) ``` #### 4. 注意事项 - **单位转换**：通常需要将浓度单位从$kg/kg$转换为$μg/m^3$，这需要空气密度数据（可从WRF输出获取或使用标准大气估算）： ```python # 假设rho_air是空气密度（kg/m³） ds['conc_ugm3'] = ds['spec001_mr'] * rho_air * 1e9 ``` - **大文件处理**：使用`chunks`参数分块读取（如`xr.open_dataset(..., chunks={'time': 10})`）以避免内存不足。 - **时间处理**：确保时间维度已正确解析，便于时间选择。 ### 三、总结 FLEXPART-WRF和传统FLEXPART的输出文件在结构上基本一致，可以使用相同的Python工具（如`xarray`和`cartopy`）进行处理和可视化。绘图主要包括浓度场平面图、垂直剖面图以及粒子轨迹动画等。 ### §§相关问题§§ 1. 如何将FLEXPART输出浓度场数据与观测站点数据进行对比验证？ 2. FLEXPART粒子轨迹数据如何转换为GIS软件可编辑的格式（如Shapefile）？ 3. 如何计算FLEXPART模拟的污染物在特定区域的沉降总量？ </think> ### FLEXPART-WRF与FLEXPART输出文件差异及绘图处理方法 #### 一、FLEXPART-WRF与FLEXPART输出文件差异 1. **气象驱动源差异**： - **FLEXPART-WRF**：使用WRF高分辨率气象场输出，包含更精细的局地环流特征[^2] - **传统FLEXPART**：通常使用全球模式（如ECMWF）输出，分辨率较低 $$ \Delta x_{\text{WRF}} \approx 1-10\text{km}, \quad \Delta x_{\text{GLOBAL}} \approx 25-50\text{km} $$ 2. **输出结构差异**： | 特征 | FLEXPART-WRF | 传统FLEXPART | |---------------------|----------------------------------|---------------------------| | **垂直坐标** | WRF地形追随坐标 | 标准气压高度 | | **地图投影** | 保留WRF原始投影 | 经纬度网格 | | **时间分辨率** | 更高（通常≤1小时） | 较低（通常3-6小时） | | **变量命名** | 相同核心变量（`spec001_mr`等） | 相同核心变量 | 3. **数据兼容性**： - **相同点**：均输出NetCDF格式，维度结构（time, height, lat, lon）一致 - **关键区别**：WRF版本包含地形相关变量（如`HGT`）和自定义投影参数[^2] #### 二、FLEXPART数据绘图处理（Python实现） ##### 1. 基础浓度场可视化 ```python import xarray as xr import matplotlib.pyplot as plt import cartopy.crs as ccrs # 读取文件（兼容两种版本） ds = xr.open_dataset('grid_conc_001.nc') # 处理WRF特殊投影 if 'WRF' in ds.attrs.get('model_version', ''): proj = ccrs.LambertConformal( central_longitude=ds.attrs['STAND_LON'], central_latitude=ds.attrs['MOAD_CEN_LAT'] ) else: proj = ccrs.PlateCarree() # 绘制地面浓度 fig = plt.figure(figsize=(12, 8)) ax = fig.add_subplot(111, projection=proj) ds['spec001_mr'].sel(height=0, time=ds.time[0]).plot.contourf( ax=ax, transform=ccrs.PlateCarree(), cmap='viridis', levels=20, add_colorbar=True ) ax.coastlines() ax.gridlines() plt.title('Surface Concentration') plt.savefig('conc_map.png', dpi=300) ``` ##### 2. 三维垂直剖面图 ```python # 选择经度剖面 cross_section = ds['spec001_mr'].sel( lon=116.4, method='nearest' ).transpose('height', 'lat') # 创建高度-纬度剖面 plt.figure(figsize=(10, 6)) cross_section.plot.contourf( y='height', levels=30, cmap='Spectral_r', cbar_kwargs={'label': 'Concentration (kg/kg)'} ) plt.ylim(0, 5000) # 限制高度范围 plt.title(f'Vertical Cross-section at Lon={116.4}°E') plt.savefig('vertical_profile.png') ``` ##### 3. 时间序列分析 ```python # 提取特定位置浓度变化 site_conc = ds['spec001_mr'].sel( lat=39.9, lon=116.4, height=0, method='nearest' ) # 转换为μg/m³（需空气密度数据） air_density = 1.225 # kg/m³ site_conc_ugm3 = site_conc * air_density * 1e9 # 绘制时间序列 plt.figure(figsize=(12, 5)) site_conc_ugm3.plot() plt.ylabel('Concentration (μg/m³)') plt.title('Beijing Urban Area Concentration Trend') plt.grid(True) plt.savefig('time_series.png') ``` ##### 4. 粒子轨迹可视化 ```python from matplotlib.animation import FuncAnimation # 读取粒子数据 part_ds = xr.open_dataset('partposit_001.nc') # 创建动态轨迹图 fig, ax = plt.subplots(figsize=(10, 8)) scatter = ax.scatter([], [], s=2) def update(frame): """更新粒子位置""" frame_data = part_ds.isel(points=slice(frame*1000, (frame+1)*1000)) scatter.set_offsets(np.c_[frame_data.longitude, frame_data.latitude]) return scatter, ani = FuncAnimation(fig, update, frames=100, interval=50) ani.save('particle_trajectories.gif', writer='pillow') ``` #### 三、WRF特殊数据处理技巧 1. **地形叠加**： ```python # 添加地形阴影 terrain = ds['HGT'].isel(time=0) ax.contourf(terrain.lon, terrain.lat, terrain, cmap='terrain', alpha=0.3, transform=ccrs.PlateCarree()) ``` 2. **投影转换**： ```python # 将WRF坐标转为经纬度 from wrf import to_np, getvar lats, lons = getvar(ds, 'lat'), getvar(ds, 'lon') ``` 3. **垂直层处理**： ```python # 获取真实高度（WRF地形追随坐标） z_agl = getvar(ds, 'height_agl') # 地面以上高度 ``` #### 四、差异处理建议 1. **坐标统一**： ```python # 将WRF投影坐标转标准经纬度 if 'Lambert' in proj.proj4_init: ds = ds.assign_coords({ 'lon': lons.mean(dim='south_north'), 'lat': lats.mean(dim='west_east') }) ``` 2. **时间对齐**： ```python # 处理WRF更高时间分辨率 if len(ds.time) > 100: # 判断为高分辨率 ds_hourly = ds.resample(time='1H').mean() ``` 3. **单位标准化**： ```python def standardize_units(ds): """统一浓度单位到μg/m³""" if 'kg/kg' in ds['spec001_mr'].attrs.get('units', ''): ds['conc_ugm3'] = ds['spec001_mr'] * 1e9 * air_density return ds ``` > **提示**：使用`ds.attrs`检查`model_version`属性确定数据来源，WRF版本通常包含`WRF`关键字[^1][^2]。 ###