做环境数据分析的人手机里十有八九存着这样一批数据某个站点的逐时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。数据量少、某个方向频率极高时适当设offset0.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提供了polarPlotpolarPlot(mydata, pollutant pm25, ws ws, wd wd)这张图的极径方向是风速圆心为静风越往外风速越大角度是风向颜色代表该风速-风向组合下的平均浓度。本质上它相当于把pollutionRose的每个方位扇区按风速再细分形成一个二维浓度场。4.2 典型浓度场怎么解读看polarPlot有经验之后基本可以快速分类下面是我常用的判读思路浓度场类型高值位置可能的解释典型污染物中心热区型图中心静风区本地源主导静稳条件下累积NO2、CO、PM2.5方向热区型某个方位、中等风速带上风向工业点源或特定排放源SO2、PM10外围环带型图外围高风速区远距离输送贡献显著O3、硫酸盐第一次拿到站点数据时建议把PM2.5、PM10、NO2、SO2、O3各画一张polarPlot横向对比。不同污染物的空间分布格局差异很大NO2多呈中心热区型交通源SO2如果呈方向热区型附近大概率有燃煤设施O3的polarPlot往往中心低、外围高因为O3是二次污染物更强风速和更充分的混合条件反而有利于生成。4.3 polarCluster用聚类自动识别污染情景polarPlot看多了之后大概率会想能不能自动把这些“风速-风向-浓度”组合分类openair提供了polarClusterclusters - 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转换时没指定tzR按系统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 常见问题速查现象大概率原因处理建议所有图时间整体偏移tz时区设置错误或缺失显式指定tz Asia/Shanghai后重新转换图是出来了但时间轴乱date列不是POSIXct或列名大小写问题用str(mydata$date)检查类型用tolower(names())规范列名玫瑰图有异常扇区风速、风向野值未清洗先用subset过滤合理范围中文标注变方块系统缺中文字体安装CJK字体后用theme设置family数据量大跑图很慢逐小时甚至分钟级全国数据按站点分批处理或先聚合再画关于openair的更多进阶功能比如后向轨迹分析、气象标准化weather normalization、源解析相关模块官方文档和Carslaw发表的论文里都有详细介绍。但先把上面这些基础流程跑通日常80%的分析需求都能覆盖。我在实际使用中还有一个体会openair最大的价值不是单个函数而是它把“时间-气象-浓度”三者统一处理的框架。气象数据本身只是一堆数字但在这个框架里风向风速被转化成污染溯源的语言时间被拆解成可以解释的尺度。工具是固定的但看数据的视角一旦打开后面很多分析都会变得顺理成章。建议刚开始接触这个包的朋友拿一个站点的数据把我上面写的整套流程跑一遍从windRose到timeVariation一两个小时就能把所有基础图摸熟。等这些函数都玩顺了再回去看之前费半天劲做出来的Excel图表你会觉得整个世界清爽了很多。