openair包实战指南：用R语言做气象数据分析与污染溯源

做环境数据分析的人，手机里十有八九存着这样一批数据：某个站点的逐时PM2.5、PM10、NO2浓度，再加上同一时次的风速、风向、温度、湿度。数据量不算大，几万行而已，但真要从中看出“污染物到底从哪来、什么时间高、这几年有没有好转”，光靠Excel和Origin折腾半天也理不出头绪。我第一次处理这类数据时，用透视表拉到怀疑人生，直到用了R语言的openair包，才算找到正路。

openair虽然定位是空气质量分析工具，但它的核心输入恰恰是气象数据——风速、风向这些要素几乎贯穿所有核心函数。这篇文章就围绕“用openair包分析气象数据”这个主题，把我实际处理监测数据的完整流程写出来，包括数据准备、风向风速可视化、污染来源方向判断、时间变化规律，以及长期趋势评估。适合正在写环境类论文、做数据分析工作的同学，也适合刚入门R但想快速出图的研究生。

1. openair包到底是什么——先搞懂它的设计逻辑

1.1 从开发背景看包的定位

openair包由英国学者David Carslaw等人开发维护，最早是围绕英国自动城市与乡村监测网络（AURN）的数据分析需求而生。它的名字拆开就是一套面向开放空气监测数据的分析工具箱。这个定位决定了它和一般的绘图包（比如ggplot2）不一样：openair不是给你一堆散装积木，而是把环境数据常见的分析场景做成了封装好的函数，你只要传入数据框、指定污染物列，就能得到一张信息量很大的图。

可能有人会问：我手里是气象数据，不是空气质量数据，用openair合适吗？答案是合适的。因为openair几乎所有核心函数都依赖风速（ws）、风向（wd）这两个气象要素做分析，污染玫瑰图、极坐标图本质上都是气象条件与污染物浓度的联合统计。哪怕你只做纯气象数据的风向风速分析，用windRose这类函数也比手动画极坐标图方便得多。

1.2 openair的数据哲学：一切围绕date列

openair对数据框只有一个硬性要求：必须有一列名为date的POSIXct时间列。注意，必须叫date，小写，类型必须是POSIXct而不是字符型或Date型。我见过不少人在这一步卡住：读入的CSV里时间是字符，直接丢给openair绘图，报错还算好的，更怕的是图能出来但时间轴完全错乱。

为什么这么强调date列？因为openair内部的cutData函数要根据时间戳自动切分出小时、星期、季节、工作日/周末这些分类变量，很多分组绘图功能（比如type = "season"）都依赖这一步。date列格式不对，后面所有按时间分组的功能全部失灵，这是我反复踩过之后的深刻体会。

1.3 安装与运行环境建议

安装很简单：

install.packages("openair") library(openair)

openair在CRAN上长期维护，直接install.packages安装即可。它依赖ggplot2、plyr、dplyr、lattice这些常用包，安装时会自动拉取。我建议在RStudio里专门建一个项目，把数据文件和脚本放在一起，跑图方便管理。R版本建议用4.x，老版本R在安装新版openair时偶尔会有依赖冲突，提前升级能省去不必要的麻烦。

2. 数据预处理——把原始监测记录变成openair认识的格式

2.1 date列：时区、格式一个都不能错

假设原始数据是CSV，里面有“时间”这一列，格式类似于"2023/1/1 0:00"。直接读进来是字符型，需要先转换：

mydata <- read.csv("station_data.csv", fileEncoding = "UTF-8") mydata$date <- as.POSIXct(mydata$date, format = "%Y-%m-%d %H:%M:%S", tz = "Asia/Shanghai")

两条命令解决，但有几个细节值得单独拎出来说。

第一，时区参数tz一定要显式指定。如果省略，R默认按系统时区解析，假如你的电脑在UTC时区或者R环境被设置成UTC，转换出来的时间会比北京时间慢8小时。所有图的时间轴都会整体平移，而且你在图上不一定能立刻察觉，等对数据的时候才发现对不上，返工成本很高。中国站点的数据，建议统一用tz = "Asia/Shanghai"。

第二，format参数必须和原始时间字符串严格匹配。%Y代表四位年份，%m是两位月份，%H:%M:%S是时分秒。如果你的数据里时间精确到分钟、或者没有秒，要相应调整格式串。实在拿不准的，可以用anytime::anytime()自动解析，但我习惯显式写format，排查问题更直接。

2.2 变量命名与单位约定

openair的绘图函数默认会去找名为ws、wd的列。如果你的列名不是这两个，要么改列名，要么在函数里显式指定：

names(mydata)[names(mydata) == "WindSpeed"] <- "ws" names(mydata)[names(mydata) == "WindDirection"] <- "wd"

也可以用参数传列名，例如windRose(mydata, ws = "WindSpeed", wd = "WindDirection")。

单位方面，openair默认风速单位是m/s，风向是0到360度的方位角（0为正北，按顺时针递增）。污染物浓度一般用µg/m³或ppb，没有硬性要求，但整个数据集要保持一致。温度、湿度这类气象要素openair也能画，比如用timePlot画温湿度时间序列，或者用scatterPlot看温湿度与浓度的关系，列名是什么都可以，调用时指定就行。

2.3 缺失值与异常值处理

监测数据里有缺失值太正常了。openair多数函数在计算均值时会自动剔除NA，所以直接保留缺失值也没关系。但有两点要注意：如果某个小时的整行都是NA，画calendarPlot时那个格子会显示空白，不影响全局；如果缺失值比例超过20%，相关时段的分析结果就要谨慎下结论。

异常值则是另一码事。比如风速出现999或者负值，风向出现-999，这种“野值”必须提前清洗，否则玫瑰图里会莫名多出一些扇区。我一般这样处理：

mydata <- subset(mydata, ws >= 0 & ws < 60) # 风速合理范围 mydata <- subset(mydata, wd >= 0 & wd <= 360) # 风向合理范围 mydata <- subset(mydata, pm25 > 0 | is.na(pm25))

第三行的写法是保留正浓度或缺失值，把负浓度去掉。负浓度在仪器观测里偶尔会出现，属于零点漂移问题，量不大直接剔除即可；如果负值比例很高，需要先做校准，不能简单删了事。

2.4 分钟数据先聚合到小时

openair很多时间分析功能（比如timeVariation）最细的粒度就是小时。如果你的原始数据是分钟级或秒级的，建议先聚合到小时再做分析。一方面很多监测网络的标况浓度本身就是小时均值，另一方面分钟级数据量太大，跑polarPlot这类密集计算函数时会明显变慢。聚合代码很简单：

library(dplyr) hourly <- mydata %>% group_by(date = cut(date, "hour")) %>% summarise(across(c(ws, wd, pm25, no2), mean, na.rm = TRUE)) %>% mutate(date = as.POSIXct(date))

风向的聚合要小心：风向是圆周变量，直接求算术平均跨越0/360度时会产生错误结果。比如350度和10度简单平均是180度，显然不对。如果风向在时段内变化剧烈，简单平均会出大问题。更稳妥的办法是用圆形平均，或者用风速的u、v分量分解后求平均风向。小时尺度内风向变化相对平稳，用算术平均多数情况下误差可控；但如果你要做日均风向，务必用circular包或自己写向量平均。

3. 风向玫瑰图与污染玫瑰图——快速锁定污染来源方位

3.1 windRose：先看风的基本面

画风向玫瑰图是openair最入门的操作：

windRose(mydata, ws = "ws", wd = "wd")

一张图出来，风速大小分布和主导风向一目了然。默认设置下，图中16个方位扇区展示了每小时内风速和风向的联合频率，颜色从蓝到绿到黄代表风速由低到高，扇区半径越长说明该风向出现频率越高。

看玫瑰图要养成一个习惯：先看主导风向，再看风速分布。比如某个城市冬季北风频率高，说明污染物容易从北方输送过来，分析重污染过程时就要重点看北方上游的排放源。如果某个方向风速普遍偏低，那么即使浓度高也不一定是外来输入，更可能是本地静稳累积。

3.2 pollutionRose：把浓度叠加到风向上

windRose只能看风本身，污染玫瑰图pollutionRose才是真正回答“污染物从哪来”的工具：

pollutionRose(mydata, pollutant = "pm25", ws = "ws", wd = "wd")

这张图上，每个方位扇区的长度代表该风向的频率，颜色代表该风向下PM2.5的平均浓度。如果一个扇区颜色偏红，说明这个方向吹来的风携带的污染物浓度高。这是典型的“源指示”信号——如果东南方向浓度明显高于其他方向，大概率东南方向有局地排放源或输入通道。

不过提醒一句：污染玫瑰图只能提示方向性，不能下因果结论。风和浓度是联合统计，不是因果关系。高浓度扇区可能真的来自上风向的源，也可能只是这个方向上风速偏低导致污染物累积。要区分这两种情形，需要用到下一节的polarPlot。

3.3 type参数：快速做分组对比

openair的分组功能几乎贯穿所有绘图函数，核心参数就是type。比如按季节分组：

windRose(mydata, ws = "ws", wd = "wd", type = "season") pollutionRose(mydata, pollutant = "pm25", ws = "ws", wd = "wd", type = "season")

type参数支持的值很多，包括"year"、"season"、"month"、"weekday"、"daylight"、"site"，以及你数据框里任意一个分类列。这个参数本质上是分面绘图，openair在内部自动按分组变量做切割。我实际用得最多的是type = "season"和type = "year"——前者看季节差异，后者看年际变化。连续三年的污染玫瑰图并排放置，某个方向的浓度贡献有没有逐年下降，一眼就能看出来。

3.4 玫瑰图的参数调节细节

默认参数能快速出图，但写论文时往往需要调整细节。常用的几个参数：

• ws.int：风速分档间隔，默认是1或2 m/s。数据风速普遍偏小时设0.5，风速大时设2或3。
• breaks：颜色断点，手动设置可以让多张图之间的色标一致。比如breaks = c(0, 1, 2, 3, 4, 6, 8)。
• offset：花瓣中心空白比例，默认0。数据量少、某个方向频率极高时，适当设offset（0.1或0.2）能让图更好看。
• paddle：花瓣形状，默认FALSE画直方条，TRUE画梭形花瓣，风格不同。
• angle：扇区角度，默认约30度，数据量少时可以加大到45度，避免扇区过碎。

每次调整参数后建议用ggsave保存PNG或PDF。openair绘图的返回值是一个列表，里面第一项就是ggplot对象：

p <- windRose(mydata, ws = "ws", wd = "wd") ggsave("wind_rose_season.png", plot = p$plot, width = 8, height = 6, dpi = 300)

pdf输出同理，把文件名后缀换成.pdf即可，投稿时矢量图是刚需。

4. polarPlot极坐标图——风速、风向与浓度的联动分析

4.1 为什么用极坐标而不是散点图

pollutionRose只能看每个方向的平均浓度，但它忽略了一个重要维度：风速。同一方向来的污染，低风速和高风速下的含义完全不同——低风速高浓度，多半是本地源排放后迅速累积；高风速高浓度，说明是上风向远距离输送过来的。为了同时展示风速和风向两个维度，openair提供了polarPlot：

polarPlot(mydata, pollutant = "pm25", ws = "ws", wd = "wd")

这张图的极径方向是风速（圆心为静风，越往外风速越大），角度是风向，颜色代表该风速-风向组合下的平均浓度。本质上，它相当于把pollutionRose的每个方位扇区按风速再细分，形成一个二维浓度场。

4.2 典型浓度场怎么解读

看polarPlot有经验之后，基本可以快速分类，下面是我常用的判读思路：

第一次拿到站点数据时，建议把PM2.5、PM10、NO2、SO2、O3各画一张polarPlot横向对比。不同污染物的空间分布格局差异很大：NO2多呈中心热区型（交通源），SO2如果呈方向热区型，附近大概率有燃煤设施；O3的polarPlot往往中心低、外围高，因为O3是二次污染物，更强风速和更充分的混合条件反而有利于生成。

4.3 polarCluster：用聚类自动识别污染情景

polarPlot看多了之后，大概率会想：能不能自动把这些“风速-风向-浓度”组合分类？openair提供了polarCluster：

clusters <- polarCluster(mydata, pollutant = "pm25", n.clusters = 4)

这个函数基于极坐标下的浓度特征做聚类，把相似的气象-浓度组合归为一类，每类对应一种典型的污染情景。输出除了聚类结果，还会给出一列时间标识，标出每个时刻属于哪一类。你可以进一步统计每一类的出现频率和平均浓度，量化“本地累积型”“远距离输送型”“静稳型”等情景各自的贡献权重。

聚类数n.clusters怎么选？我的经验是普通城市站点3到5类比较合适，太少类别含糊，太多则场景碎片化，解释起来很费力。openair会提供多聚类数方案的汇总图，结合组内差异下降的趋势来判断，不用拍脑袋硬定。

4.4 参数调整与统计方法

polarPlot默认用均值作为浓度统计量，但对异常值比较敏感。如果数据集里有少量极端高值，建议改用中位数：

polarPlot(mydata, pollutant = "pm25", statistic = "median")

statistic还可以设成"max"、"frequency"、"stdev"等。stdev模式看的是浓度变异性，有时能发现均值图上被掩盖的“脉冲式”污染过程。

另外一对参数是limits和cols，控制色标范围与配色。多张polarPlot对比时务必统一limits，否则色标范围不一致，同一颜色代表的浓度完全不同，放在一起对比容易得出错误结论：

polarPlot(mydata, pollutant = "pm25", limits = c(0, 100), cols = "YlOrRd")

多站点、多月对比时，这个习惯尤其重要。

5. 时间维度上的发现——日变化、日历热图与长期趋势

5.1 timeVariation：早高峰晚高峰一看便知

气象数据里的时间规律是重头戏，openair的timeVariation函数把时间拆成三个层级展示：

timeVariation(mydata, pollutant = "pm25")

输出是三张并列的小图：左上显示一天24小时的浓度变化，右上显示一周7天的变化，下方显示12个月的逐月变化。这张图几乎是我做任何环境数据分析时的第一步，能快速回答几个核心问题：峰值在几点，是早高峰还是夜间边界层变化；工作日与周末有没有差别，交通源特征是否明显；冬季和夏季的差异有多大，采暖与光化学反应的影响如何。

一个容易被忽略的用法是对污染物做差分分析，比如比较工作日和周末的NO2日变化曲线，看早高峰是否消失或推迟，这是判断交通源贡献的常用手段。实现方式也很简单，用type参数分组就行：

timeVariation(mydata, pollutant = "no2", type = "season")

5.2 calendarPlot与trendLevel：日历热图和月-小时热图

如果想知道“过去一年里哪些天浓度特别高”，calendarPlot是最直观的选择：

calendarPlot(mydata, pollutant = "pm25")

它把数据排成一个月一个格子的日历热图，横轴是星期，纵轴是一个月内的天数，每个格子的颜色代表当天的日均浓度。重污染天气在图上会形成明显的红块，连续几天的红块通常对应一次完整的重污染过程。这张图在答辩汇报和论文里很出效果，“一目了然”是它最大的优势。

配合trendLevel还能看更细的时间交互，它画的是“月-小时”二维热图：

trendLevel(mydata, pollutant = "pm25")

横轴是月份，纵轴是24小时，颜色代表该月该小时的浓度均值。你能直观看到污染高值集中在哪些月份、哪些时段。北方城市经常出现“冬季夜间高、夏季午后低”的格局，这背后是采暖排放在静稳边界层里累积的结果，trendLevel一张图就能说清楚。

5.3 长期趋势是否真实改善：TheilSen与smoothTrend

评估“这几年污染是不是真的在好转”，是很多报告的核心结论。直接用逐日均值画趋势线噪声太大，肉眼很难判断。openair提供了两个函数：

TheilSen(mydata, pollutant = "pm25", avg.time = "month") smoothTrend(mydata, pollutant = "pm25")

TheilSen基于Theil-Sen估计器计算趋势斜率，对异常值和离群点不敏感，输出里包含斜率的置信区间。如果置信区间不跨越0，说明趋势在统计上显著。smoothTrend画出平滑曲线和置信带，能够一眼看出浓度在哪些年份真的下降、哪些年份在反弹。

这里提醒一句：趋势分析的时间尺度很重要。做年际趋势时务必先把数据聚合到月或季度再跑TheilSen，直接拿小时数据跑会因自相关太强而产生偏差。另外，如果研究期间监测站点位置变更或仪器更换过，趋势结果要谨慎解释，最好在报告里注明这些背景信息。

6. 分组对比与多站点分析的几个实用操作

6.1 用selectByDate精准切取时间段

有时候只关心某段时间，比如采暖期：

heating <- selectByDate(mydata, start = "2023-11-15", end = "2024-03-15")

selectByDate还支持hour参数，可以只看早晚高峰时段：

peak <- selectByDate(mydata, hour = c(7, 8, 9, 17, 18, 19))

这个功能在交通源分析里很实用，可以快速把高峰时段的数据单独抽出来画污染玫瑰图，看看污染来向和非高峰时段有什么不同。

6.2 多站点数据合并与site分组

如果手里有多个站点的数据，直接纵向合并再加一列站点标识：

all_data <- rbind(site_a, site_b, site_c) all_data$site <- rep(c("A", "B", "C"), times = c(nrow(site_a), nrow(site_b), nrow(site_c)))

之后所有绘图函数都能用type = "site"分面。比如：

timeVariation(all_data, pollutant = "pm25", type = "site") polarPlot(all_data, pollutant = "pm25", type = "site")

多站点对比最大的价值在于识别空间差异。城区站和郊区站的windRose可能差不多，但polarPlot会明显不同：城区站中心热区型，郊区站方向热区型。这种空间信息对解读区域污染成因非常有帮助。

6.3 把基础作图封装成自己的函数

当分析流程固定下来后，强烈建议写一个自己的封装函数，避免每次重复调整参数。比如我做年度对比时常用：

plot_season_rose <- function(df, pollutant = "pm25") { pollutionRose(df, pollutant = pollutant, ws = "ws", wd = "wd", type = "season", cols = "YlOrRd") }

封装的意义不只是省几行代码，更重要的是保证所有图风格一致。写论文时几十张图统一色标、统一尺寸，会省掉大量后期修图时间。这是我从第二次做完整项目开始才养成的习惯，早期每张图都要单独调参，反反复复改到崩溃。

7. 踩坑记录与新手避雷指南

7.1 date列时间偏移8小时

我遇到最多的问题就是时间偏移。症状是timeVariation画出来的日变化曲线和实际逐时数据对不上，或者calendarPlot里某些天的颜色怪怪的。原因几乎都是时区问题：CSV里时间是北京时间，但as.POSIXct转换时没指定tz，R按系统UTC解析后所有时间慢了8小时。

排查方法很简单，转换后立刻检查：

range(mydata$date)

看起止时间和原始数据是否一致。不一致就重新设置时区再转，不要手动往时间戳上加8小时，那种补丁式修法后续还会埋雷。

7.2 列名冲突与大小写

openair对列名比较挑剔，尤其是date列必须是精确小写。如果你的数据里同时有"Date"和"date"两列，R在某些函数里会用错列，画出来的图时间轴完全看不懂。解决方法是数据读入后第一时间做列名规范化：

names(mydata) <- tolower(names(mydata))

7.3 中文字体显示问题

openair默认绘图字体对中文支持不好，图里的中文标题或站点名经常变成方框。两个解决办法：图里尽量用英文或拼音；画完后用ggplot2的theme设置中文字体：

p <- timeVariation(mydata, pollutant = "pm25") p$plot + theme(text = element_text(family = "Noto Sans CJK SC"))

如果系统里没有合适字体，先安装"Noto Sans CJK"或"WenQuanYi"，然后重启RStudio。这个坑在Windows上尤其常见，Mac上会好一些。

7.4 大数据量的性能处理

openair处理几十万行的数据毫不费力，但如果你有数千万行（比如全国多站点多年的分钟数据），polarPlot、timeVariation这些涉及平滑和置信区间计算的函数就会明显变慢。我的建议是：先按站点或时间范围分批跑通分析思路，思路确定后再用全量数据出最终图。

7.5 常见问题速查

关于openair的更多进阶功能，比如后向轨迹分析、气象标准化（weather normalization）、源解析相关模块，官方文档和Carslaw发表的论文里都有详细介绍。但先把上面这些基础流程跑通，日常80%的分析需求都能覆盖。

我在实际使用中还有一个体会：openair最大的价值不是单个函数，而是它把“时间-气象-浓度”三者统一处理的框架。气象数据本身只是一堆数字，但在这个框架里，风向风速被转化成污染溯源的语言，时间被拆解成可以解释的尺度。工具是固定的，但看数据的视角一旦打开，后面很多分析都会变得顺理成章。建议刚开始接触这个包的朋友，拿一个站点的数据把我上面写的整套流程跑一遍，从windRose到timeVariation，一两个小时就能把所有基础图摸熟。等这些函数都玩顺了，再回去看之前费半天劲做出来的Excel图表，你会觉得整个世界清爽了很多。