## 1. 从零开始：为什么选择Python、Cartopy和Pysplit？ 如果你正在研究大气污染物的来源，或者想搞清楚某次雾霾天气的“元凶”到底来自何方，那么“反向轨迹分析”绝对是你绕不开的技术。简单来说，它就像给空气粒子装上了一个“时光倒流GPS”，我们可以模拟出在过去几天里，到达我们监测点的空气团是从哪里飘过来的。传统的做法往往依赖一些专门的商业软件或者命令行工具，操作繁琐，可视化效果也一言难尽。今天，我就来分享一套我用了好几年的“组合拳”：用Python生态里的Cartopy和Pysplit，来打造一个高效、美观且完全免费的反向轨迹分析与可视化工作流。 几年前，大家做这个可能首选Basemap。但Basemap已经停止维护很久了，在Python 3的新版本上安装就是个头疼事，各种依赖冲突。Cartopy是它的“官方继任者”，基于更现代的PROJ和GEOS库，地图投影更准确，性能也更好，关键是安装起来顺滑多了。而Pysplit这个库，它本质上是一个Python接口，帮你封装了调用NOAA的HYSPLIT模型（大气轨迹和扩散模型的标准）的复杂过程。你不用再去手动写那些晦涩的HYSPLIT控制文件，用Python几行代码就能批量生成轨迹。 这套组合的优势在哪呢？第一是**全流程自动化**。从下载气象数据、计算轨迹到绘图，全部可以用Python脚本搞定，特别适合需要处理大量时间序列数据的科研或业务场景。第二是**可视化自由度极高**。Cartopy提供了强大的地理绘图能力，你可以轻松定制地图范围、添加省界、河流、标注站点，做出出版级质量的图。第三是**可复现性强**。所有步骤都是代码，换个时间、换个地点，改几个参数重新跑一遍就行，完全避免了手动操作带来的误差和不一致。 我自己的项目里，从研究华北平原秋冬季霾的输送通道，到分析沿海城市臭氧污染的区域贡献，这套工具链都立下了汗马功劳。接下来，我就手把手带你走通整个流程，从环境搭建到出图，把踩过的坑和总结的技巧都告诉你。 ## 2. 搭建你的分析环境：告别Basemap，拥抱Cartopy 工欲善其事，必先利其器。第一步就是把我们的Python环境配置好。我强烈建议使用Conda来管理环境，它能很好地处理地理信息库复杂的二进制依赖。 ### 2.1 创建并激活虚拟环境 打开你的终端（Windows用Anaconda Prompt，Mac/Linux用终端），执行以下命令。这里我用Python 3.9，因为它和各个库的兼容性非常稳定。 ```bash conda create -n air_trajectory python=3.9 conda activate air_trajectory ``` ### 2.2 安装核心绘图库Cartopy 这是最关键的一步，也是和旧教程最大的区别。**直接用Conda安装Cartopy**，它会自动处理好PROJ、GEOS、Shapely等一系列底层依赖，这是最省心的方式。 ```bash conda install -c conda-forge cartopy ``` `-c conda-forge` 指定从conda-forge这个频道安装，这里的库通常更新更及时。安装完成后，你可以在Python里试试 `import cartopy`，没有报错就成功了。 ### 2.3 安装其他必需依赖库 除了Cartopy，我们还需要一些帮手。 ```bash # 使用conda或pip安装均可，conda优先 conda install -c conda-forge xarray netCDF4 matplotlib scipy pip install pysplit openpyxl ``` - **xarray, netCDF4**: 用来读取和处理气象数据（通常是.nc格式的NetCDF文件），这是GDAS1数据的标准格式。 - **matplotlib**: Python绘图的基础，Cartopy是基于它扩展的。 - **scipy**: 提供一些数学工具，可能在后续数据处理中用到。 - **pysplit**: 我们的主角之一，用于生成和操作HYSPLIT轨迹。 - **openpyxl**: 有时候污染源数据是Excel格式，用它来读取。 **关于Basemap和Pysplit的兼容性问题**：Pysplit的早期版本确实依赖Basemap来绘图。但我们的策略是**只用Pysplit来计算轨迹，用Cartopy来绘图**。因此，我们完全不需要安装Basemap。如果你发现导入Pysplit时提示缺少Basemap，不用担心，Pysplit的核心计算功能不受影响。如果实在想消除警告，可以找到Pysplit安装目录下的 `mapdesigner.py` 文件，将其中的 `from mpl_toolkits.basemap import Basemap` 以及 `MapDesign`类中相关的绘图方法注释掉即可，因为我们不会调用它们。 ### 2.4 安装HYSPLIT模型核心 Pysplit只是一个Python外壳，真正的轨迹计算引擎是NOAA的HYSPLIT模型。我们需要安装它。 1. **访问NOAA ARL官网**，找到HYSPLIT下载页面。 2. 选择对应你操作系统的版本（Windows或Linux）。对于Windows用户，通常下载一个预编译的安装包，直接运行安装程序，记住安装路径（例如 `C:\hysplit`）。 3. 对于Linux用户，可以下载预编译的压缩包，解压后就能用。如果需要图形界面，还得额外安装Tcl/Tk和ImageMagick，但我们的自动化流程一般用命令行模式，可以不装。 安装好后，记下HYSPLIT的工作目录路径，后面配置Pysplit时会用到。 ### 2.5 获取地形数据 为了让轨迹计算更准确（尤其是考虑地形对气流的影响），我们需要一个地形高程数据文件 `ETOPO2v2c_f4.nc`。这个文件可以在NOAA官网找到。下载后，把它放在一个固定的目录，比如你项目下的 `./conf/altitude_data/` 文件夹里。在后续的轨迹生成配置中，需要指定这个文件的路径。 至此，我们的软件环境就准备好了。可以看到，核心就是Cartopy+Pysplit+HYSPLIT，思路清晰，没有历史包袱。 ## 3. 获取“燃料”：GDAS1气象数据详解与自动下载 计算空气轨迹，就像做一道物理计算题，除了公式（HYSPLIT模型），还需要初始条件（监测点位置、时间）和驱动数据（气象场）。GDAS1数据就是那个最重要的“驱动数据”，它提供了全球范围的气压、温度、风速、风向等网格化气象数据。 ### 3.1 认识GDAS1数据 GDAS1是NOAA全球数据同化系统的1度分辨率产品，覆盖全球，时间分辨率是6小时（00， 06， 12， 18 UTC）。数据格式是NetCDF，文件名有固定规律：`gdas1.[月份英文缩写3位][年份后2位].w[周数]`。例如： - `gdas1.jan23.w1`：2023年1月第1周的数据。 - `gdas1.apr05.w4`：2005年4月第4周的数据。 你需要为你想要模拟的轨迹的**每个时间点**，准备对应的GDAS1数据文件。比如你要模拟过去72小时内每小时的轨迹，那么你可能需要包含这72小时在内的多个6小时间隔的数据文件。 ### 3.2 手动与自动下载 **手动下载**：你可以直接访问NOAA ARL的数据存档网站，找到对应的年份、月份和周数，点击下载。这对于一次性获取少量数据很方便。 **自动下载（推荐）**：对于长期监测或需要大量历史数据的情况，写个Python脚本自动抓取是必须的。NOAA提供了FTP服务器。下面是我常用的一个简化版下载函数，你可以把它保存为 `download_gdas1.py`： ```python import ftplib import os from datetime import datetime, timedelta def download_gdas1_week(year, month, week, local_dir='./gdas1_data'): """ 下载指定年、月、周的GDAS1数据。 参数示例：year='23', month='jan', week='1' """ os.makedirs(local_dir, exist_ok=True) ftp = ftplib.FTP('arlftp.arlhq.noaa.gov') ftp.login() # 匿名登录 ftp.cwd('archives/gdas1') filename = f'gdas1.{month}{year}.w{week}' local_path = os.path.join(local_dir, filename) print(f'正在下载 {filename} ...') with open(local_path, 'wb') as f: ftp.retrbinary(f'RETR {filename}', f.write) print(f'下载完成: {local_path}') ftp.quit() # 示例：下载最近7天的数据（通常对应最新一周） # 需要先根据当前日期计算出对应的年、月、周 ``` 更实用的脚本会包含日期循环、错误重试、以及将“最近7天”这种表述转换为标准文件名格式的逻辑。网上也有很多现成的脚本可以参考，核心就是利用 `ftplib` 连接 `arlftp.arlhq.noaa.gov`，进入 `/archives/gdas1` 目录进行下载。 ## 4. 生成轨迹文件：用Pysplit驱动HYSPLIT 有了气象数据，我们就可以开始计算轨迹了。Pysplit的核心价值在这里体现：它把HYSPLIT复杂的参数配置和控制文件生成工作，变成了简单的Python函数调用。 ### 4.1 准备配置文件 我习惯把所有的运行参数写在一个JSON配置文件里，这样管理和修改起来非常清晰。创建一个 `config_backtraj.json` 文件： ```json { "storage_dir": "./trajectory_results", "working_dir": "C:/hysplit/working", "meteo_dir": "./gdas1_data", "altitude_data": "./conf/altitude_data/ETOPO2v2c_f4.nc", "location": [38.03, 114.48], "altitudes": [500, 1000, 1500], "runtime": -72, "basename": "sjz", "years": ["23"], "months": ["jun"], "hours": [0, 6, 12, 18], "dates": [15, 16, 17, 18] } ``` 我来解释一下关键参数： - `working_dir`: HYSPLIT模型的工作目录，需要你有写入权限。 - `meteo_dir`: 你存放GDAS1数据的文件夹路径。 - `location`: 起始点的经纬度 `[纬度, 经度]`，例如石家庄。 - `altitudes`: 起始高度列表（米AGL）。计算多个高度的轨迹可以分析不同垂直层的气流。 - `runtime`: 轨迹运行时间（小时）。负数表示反向（向后追溯），正数表示正向（向前预测）。分析污染来源用负数，比如 `-72`。 - `basename`: 输出轨迹文件的前缀。 - `years`, `months`, `hours`, `dates`: 定义了你要在哪些具体时刻启动轨迹计算。这个例子是计算2023年6月15、16、17、18日，每天UTC时间0、6、12、18时，共4天*4时次=16组轨迹，每组轨迹又包含500、1000、1500米3个高度，最终会生成48条轨迹。 ### 4.2 批量生成轨迹 使用Pysplit的 `bulk_trajgen` 模块可以轻松实现批量生成。下面是一个核心脚本： ```python import json from pysplit import bulk_trajgen # 加载配置 with open('config_backtraj.json', 'r') as f: config = json.load(f) # 调用批量生成函数 traj_files = bulk_trajgen.get_trajectories( working_dir=config['working_dir'], storage_dir=config['storage_dir'], meteo_dir=config['meteo_dir'], altitude_data_path=config['altitude_data'], years=config['years'], months=config['months'], hours=config['hours'], altitudes=config['altitudes'], coordinates=config['location'], runtime=config['runtime'], basename=config['basename'], dates=config['dates'] ) print(f"轨迹生成完成！文件保存在：{config['storage_dir']}") ``` 运行这个脚本，Pysplit就会在后台为你调用HYSPLIT可执行程序，依次生成所有轨迹。生成的轨迹文件是文本格式（.txt或.traj），里面记录了轨迹点上每个时间步的经纬度、高度、气压等信息。这个过程可能会花点时间，取决于你计算轨迹的数量。在我的经验里，批量处理成百上千条轨迹是常事，自动化脚本节省了大量手动操作时间。 ## 5. 可视化轨迹：用Cartopy绘制专业地图 计算出的轨迹数据是文本，不直观。接下来就是Cartopy大显身手的时候了。我们将把一条条轨迹线画到地图上。 ### 5.1 读取轨迹数据 Pysplit提供了读取轨迹文件的工具函数。 ```python from pysplit import traj_io # 假设我们生成了一条轨迹文件路径 traj_file = './trajectory_results/sjz2023061500_500.txt' # 加载轨迹 trajectory = traj_io.read_hysplit_trajectory(traj_file) # trajectory对象包含了很多信息，比如： lats = trajectory['latitudes'] # 纬度序列 lons = trajectory['longitudes'] # 经度序列 alts = trajectory['altitudes'] # 高度序列 times = trajectory['timestamps'] # 时间序列 ``` ### 5.2 使用Cartopy绘制单条轨迹 现在，我们用Cartopy来画图。Cartopy的强大在于它内置了各种地图投影和地理特征。 ```python import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature # 1. 创建地图和坐标系 fig = plt.figure(figsize=(12, 8)) # 使用PlateCarree投影（等经纬度），这是最常用的 ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) # 2. 设置地图范围（聚焦于东亚区域） ax.set_extent([70, 140, 15, 55], crs=ccrs.PlateCarree()) # [lon_min, lon_max, lat_min, lat_max] # 3. 添加地理特征，让地图更丰富 ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.3) # 陆地 ax.add_feature(cfeature.OCEAN, facecolor='lightblue', alpha=0.3) # 海洋 ax.add_feature(cfeature.COASTLINE, linewidth=0.5) # 海岸线 ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5) # 国界线 ax.add_feature(cfeature.RIVERS, edgecolor='blue', linewidth=0.5, alpha=0.5) # 河流 # 添加省界（需要额外数据，如中国省界shp文件，这里用Cartopy自带的州省界示意） # ax.add_feature(cfeature.STATES, linewidth=0.3) # 4. 绘制轨迹线 # 将经纬度数据转换为地图投影坐标 ax.plot(lons, lats, transform=ccrs.PlateCarree(), color='red', linewidth=2, marker='o', markersize=4, label='500m AGL') # 5. 标记起始点 ax.scatter(lons[0], lats[0], transform=ccrs.PlateCarree(), color='green', s=100, zorder=5, label='Start Point (Shijiazhuang)') # 6. 添加网格和标签 ax.gridlines(draw_labels=True, dms=True, x_inline=False, y_inline=False, linestyle='--', alpha=0.5) ax.set_title('72-Hour Backward Trajectory (2023-06-15 00:00 UTC)', fontsize=16) ax.legend(loc='upper right') plt.tight_layout() plt.savefig('./trajectory_plot.png', dpi=300, bbox_inches='tight') plt.show() ``` 这段代码生成了一张专业的轨迹图。你可以自由调整颜色、线型、标记符号，也可以轻松地在一张图上叠加多条不同高度或不同时间的轨迹，用颜色或线型区分它们，直观地对比气流路径的差异。 ### 5.3 进阶可视化：轨迹簇与频率分析 在实际研究中，我们往往不是看单条轨迹，而是分析一段时间内（比如一个月）的所有轨迹，形成“轨迹簇”，或者计算轨迹的端点频率（轨迹最终从哪里来），这能更好地揭示主要的输送通道。 Pysplit和Cartopy结合，可以方便地实现。思路是：循环读取所有轨迹文件，将每条轨迹的路径点（或端点）收集起来，然后用散点图或密度图的形式画在地图上。Cartopy的 `ax.scatter` 和 `ax.hist2d` 函数在这里非常有用。你可以用颜色深浅或点的大小来表示频率高低，一张图就能清晰展示污染物的主要潜在来源区域。 ## 6. 深入分析：潜在源贡献（PSCF）计算与绘图 轨迹图告诉我们“气团从哪来”，而潜在源贡献值分析则进一步回答“哪些来源区域对高污染浓度的贡献更大”。PSCF的基本思想是：如果一个网格区域有越多的轨迹在污染时段经过，且经过时的污染物浓度越高，那么这个区域的PSCF值就越高，它作为污染源的可能性就越大。 ### 6.1 PSCF计算原理简述 1. **网格化研究区域**：将地图划分成许多小网格（比如0.5°×0.5°）。 2. **统计轨迹端点**：对于每一条在污染事件发生时计算的反向轨迹，记录其每个时间步（比如每小时）的端点落在哪个网格里。 3. **区分污染与非污染轨迹**：根据监测点的污染物浓度是否超过某个阈值，将轨迹分为“污染轨迹”和“清洁轨迹”。 4. **计算比值**：对于每个网格，计算 `PSCF_ij = M_ij / N_ij`。其中，`N_ij` 是所有轨迹端点落在该网格的次数，`M_ij` 是“污染轨迹”端点落在该网格的次数。这个比值越高，说明该网格是潜在污染源的可能性越大。 5. **权重函数**：为了减少由于`N_ij`过小（统计样本少）带来的不确定性，通常会引入一个权重函数 `W_ij` 来修正PSCF值，得到 `WPSCF_ij = PSCF_ij * W_ij`。 ### 6.2 利用修改的pyPSCF进行计算 完全从零实现PSCF算法涉及大量的轨迹数据处理和网格统计。幸运的是，开源社区有 `pyPSCF` 这样的项目，它提供了一个GUI工具，内部集成了轨迹生成和PSCF计算。我们可以借鉴或修改其PSCF计算的核心部分。 你需要准备两份数据： 1. **轨迹文件集合**：就是上一节我们批量生成的所有轨迹文件。 2. **污染浓度时间序列**：一个包含时间戳和污染物浓度（如PM2.5）的表格（CSV或Excel）。时间需要与你生成轨迹的启动时间对应（注意时区转换，通常轨迹用UTC，监测数据用本地时间）。 一个简化的工作流程如下： ```python # 伪代码，展示逻辑 def calculate_pscf(traj_files_dir, pollution_data_path, threshold): # 1. 读取所有轨迹，提取端点 all_endpoints = [] polluted_endpoints = [] for traj_file in traj_files: traj = read_trajectory(traj_file) traj_time = traj.start_time # 轨迹启动时间 # 找到对应时间的污染浓度 conc = get_pollution_concentration(pollution_data_path, traj_time) # 提取该轨迹所有端点（或特定时间段的端点） endpoints = extract_endpoints(traj) all_endpoints.extend(endpoints) if conc > threshold: polluted_endpoints.extend(endpoints) # 2. 定义网格 grid_lons = np.arange(70, 140, 0.5) grid_lats = np.arange(15, 55, 0.5) # 3. 统计每个网格的N和M N_grid, _, _ = np.histogram2d([e[0] for e in all_endpoints], [e[1] for e in all_endpoints], bins=[grid_lons, grid_lats]) M_grid, _, _ = np.histogram2d([e[0] for e in polluted_endpoints], [e[1] for e in polluted_endpoints], bins=[grid_lons, grid_lats]) # 4. 计算PSCF并应用权重 pscf_grid = M_grid / N_grid pscf_grid[N_grid < sample_threshold] = np.nan # 样本太少设为NaN # 应用权重函数W... wpscf_grid = apply_weight(pscf_grid, N_grid) return wpscf_grid, grid_lons, grid_lats ``` ### 6.3 绘制PSCF分布图 计算出WPSCF网格数据后，用Cartopy绘制填色图是展示结果的最佳方式。 ```python import numpy as np import matplotlib.pyplot as plt import cartopy.crs as ccrs import cartopy.feature as cfeature from matplotlib import colors from matplotlib.cm import ScalarMappable # 假设我们已经得到了 wpscf_grid, grid_lons, grid_lats # grid_lons和grid_lats是网格中心点的坐标向量 fig = plt.figure(figsize=(14, 10)) ax = fig.add_subplot(1, 1, 1, projection=ccrs.PlateCarree()) ax.set_extent([70, 140, 15, 55], crs=ccrs.PlateCarree()) # 添加基础地理特征 ax.add_feature(cfeature.LAND, facecolor='lightgray', alpha=0.2) ax.add_feature(cfeature.COASTLINE, linewidth=0.8) ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.8) # 创建网格（用于pcolormesh） lon_grid, lat_grid = np.meshgrid(grid_lons, grid_lats) # 定义颜色映射和标准化 # PSCF值范围通常在0-1之间 cmap = plt.cm.YlOrRd # 黄-橙-红渐变色，适合表示风险程度 norm = colors.Normalize(vmin=0, vmax=0.8) # 根据你的数据调整vmax # 绘制填色图 # 注意：pcolormesh要求网格是二维的，且数据点与网格角点对应。 # 这里我们假设wpscf_grid是(lat, lon)形状，且grid_lons, grid_lats是中心点。 # 实际使用时可能需要调整，或使用contourf。 cf = ax.pcolormesh(lon_grid, lat_grid, wpscf_grid.T, # 注意转置以适应Cartopy的维度顺序 cmap=cmap, norm=norm, transform=ccrs.PlateCarree(), alpha=0.8) # 添加颜色条 cbar = plt.colorbar(ScalarMappable(norm=norm, cmap=cmap), ax=ax, orientation='horizontal', pad=0.05, aspect=40, shrink=0.8) cbar.set_label('Weighted Potential Source Contribution Function (WPSCF)', fontsize=12) ax.set_title('Potential Source Contribution of PM2.5 for Shijiazhuang (Summer 2023)', fontsize=16, pad=20) ax.gridlines(draw_labels=True, linestyle='--', alpha=0.5) plt.tight_layout() plt.savefig('./pscf_map.png', dpi=300, bbox_inches='tight') plt.show() ``` 这张图上的暖色区域（红色）就是需要重点关注的可能污染源区。结合地理信息，你就能解读出污染物主要的区域输送路径和潜在贡献源了。 ## 7. 实战经验与避坑指南 最后，分享几个我在实际项目中积累的经验，希望能帮你少走弯路。 **第一，时间对齐是重中之重。** HYSPLIT模型使用的气象数据（GDAS1）是UTC时间，而你本地的监测数据很可能是北京时间（UTC+8）。在配置轨迹生成参数（`config_backtraj.json`中的`hours`）和匹配污染浓度时，一定要统一时间基准。我建议**全部使用UTC时间**来思考和配置，这样可以避免很多混乱。在最终绘图标注时，再根据需求转换为本地时间。 **第二，气象数据要完整覆盖轨迹时段。** 计算一条72小时的反向轨迹，不仅需要轨迹启动时刻的气象数据，还需要之前72小时的数据。因为模型是积分运算。所以，你需要确保你的GDAS1数据文件在时间上是连续的，能够覆盖你所有轨迹计算的时间窗口。下载数据时最好多下载前后几天的数据作为缓冲。 **第三，计算资源与效率。** 批量计算成千上万条轨迹是CPU密集型任务。如果你的`working_dir`放在机械硬盘上，频繁的I/O可能会成为瓶颈。尽量将工作目录设置在固态硬盘上。另外，Pysplit的 `bulk_trajgen` 是顺序执行的，对于超大规模计算，可以考虑用Python的 `multiprocessing` 模块进行并行化，但要注意HYSPLIT本身是否支持并行以及你的许可证限制。 **第四，Cartopy绘图优化。** 绘制包含大量轨迹或高分辨率PSCF网格的地图时，可能会比较慢。可以尝试以下方法：1) 在调试阶段使用较低分辨率的地图数据（Cartopy的 `cfeature` 有 `scale` 参数控制）；2) 对于最终出图，可以先保存为矢量格式（如PDF、SVG），再转换为位图，这样清晰度更高；3) 使用 `ax.contourf` 代替 `ax.pcolormesh` 绘制PSCF，有时效率更高且自带平滑效果。 **第五，结果的解读要结合气象学知识。** 轨迹和PSCF是强大的工具，但它们给出的是统计和模拟结果。一个高PSCF值的区域，不一定就是真正的污染源。需要结合当地的地形（山脉、河谷）、盛行风场、污染源排放清单（哪些地方有工厂、城市群）进行综合判断。工具帮你缩小了范围，但最终的结论需要严谨的科学分析。 这套从数据下载、轨迹计算到高级分析和可视化的完整流程，我已经在多个环境空气质量分析项目中验证过其稳定性和有效性。它最大的魅力在于将复杂的专业模型和精美的科学绘图，都整合到了灵活的Python脚本中，让你能更专注于科学问题本身，而不是繁琐的软件操作。希望这篇文章能帮你顺利上手，开启你的大气轨迹分析之旅。如果在实践中遇到具体问题，多查阅Cartopy和Pysplit的官方文档，以及HYSPLIT的用户手册，里面有很多细节参数值得深入探索。