搞过几年环境数据分析的朋友应该都绕不开一个名字openair。这个R包第一次出现在我面前是我刚接手某站点PM2.5数据、急着画风玫瑰图的时候。说实话google大法搜了很多现成代码结果不是字体别扭就是方向角算错直到同事甩过来一句“你直接用openair啊”。这个包全名就是openair英国Kings College London团队开发的R包专攻空气质量与气象数据的分析、可视化和趋势检验。你可能已经见过了它标志性的风玫瑰图、污染玫瑰图还有那个四联图一样的时间变化面板。但openair远不止“画图漂亮”这么简单它把很多气象分析中的脏活累活——比如风向风速的区间统计、浓度数据的周期分解、甚至曼肯德尔趋势检验——都封装成了几行代码就能调的“傻瓜函数”。这篇文章我打算从使用者的真实视角出发把openair从安装、数据规范化、高频画图到进阶统计完整讲一遍。不仅讲怎么用还会解释它背后的设计逻辑以及哪些环节容易踩坑。适合刚入行的环境数据分析师、做毕业论文的学生和那些想快速摸清一个站点气象与污染关系的人。1. 内容整体设计与思路拆解1.1 openair包能做什么不只是画风玫瑰图很多人认识openair是从风玫瑰图开始的这一点不奇怪因为windRose()这个函数实在太出名了。但如果你只把它当“画图工具”那真是低估了这个包。openair的能力可以归成五大块基础可视化风玫瑰、污染玫瑰、时间序列、周期与规律分析日变化、周变化、月变化、源解析辅助污染概率函数CPF、条件风玫瑰、趋势与显著性检验TheilSen斜率、Mann-Kendall检验、模型评估与校正散点图、泰勒图、统计指标计算。这五块基本覆盖了日常环境数据分析90%以上的需求。我自己用得最多的场景就是给某个污染过程做“快速体检”。比如某天NO2浓度突然飙升我需要快速知道当时风从哪来风速多大是不是静稳天气这种浓度升高是局地排放还是外来传输这一套问题用timePlot()、windRose()和pollutionRose()三件套基本十分钟就能有个初步结论。openair适合的人群也相当广。正在写毕业论文的学生可以用它快速出图环境监测站的技术人员可以用它做季度趋势汇报做大气科研的同志可以用它的趋势检验和模型评估功能。更难得的是这个包虽然脱胎于英国项目但对中国的数据格式同样友好只要你把时间列和变量列整理成规范格式就行。1.2 为什么用openair而不是自己写代码懒人哲学与统计机制用过ggplot2画风玫瑰图的人都知道这不是画一个漂亮圆盘那么简单。你得把风向切分成若干个区间比如每30度一个扇区然后统计每个扇区内风速的分布再按百分比着色。数据清洗、区间计算、坐标转换、图例设置每一步都要小心否则最后出来的图要么方向错位、要么频率比例失真。openair的做法是你只需要提供原始的风向、风速和污染物浓度列函数内部会自动完成扇区划分、统计、圆形坐标映射和图例渲染。更重要的是它不只是“画”一个图形还“计算”了一个统计结果比如pollutionRose()默认输出的是每个风向区间内污染物浓度的均值你可以通过statistic参数切换成“中位数”“分位数”“CPF概率”等。这背后是有学术源流的很多城市尺度的大气污染源方向判断靠的就是这种分扇区统计的逻辑。再说得直白一点openair的设计哲学是“让分析者把精力放在解释数据上而不是耗在写绘图代码上”。对大多数非编程出身的环保从业者来说这正是最需要的东西。哪怕你平时只会read.csv()和summary()学一点openair也能立刻上手做专业分析。2. 环境准备与数据规范化让openair兼容你的数据2.1 安装与加载新旧版本都要稳安装openair第一步在RStudio的控制台里运行install.packages(openair)正常情况下会自动装上依赖包。但R生态经常出现一个老问题你的R版本如果偏旧而openair更新版本要求更高版本的R这时候直接安装可能失败报错一般长得像“package ‘openair’ is not available for this version of R”。解决办法有两个。第一先升级R到最新版这是最省心的路第二如果因为某些原因不方便升级可以去CRAN的archive页面找旧版本的openair源码包然后用install.packages(openair_x.x.x.tar.gz, repos NULL, type source)手动安装前提是系统里有RtoolsWindows或编译工具链macOS/Linux。还有一个我常遇到的问题装完openair之后加载时报错说找不到某个依赖包比如latticeExtra、hexbin这些比较冷门的包。遇到这种情况不用慌单独补装就行install.packages(c(latticeExtra, hexbin, clustreg)) library(openair)加载成功后建议先看一眼包自带的示例数据方便我们验证函数是否正常运行data(mydata) head(mydata)这个mydata是openair内置的一个伦敦地区示例数据集里面有date、ws风速、wd风向、no2、o3、pm10等常用列。后面所有实操我建议你也拿这份数据先跑一遍熟悉了再换自己的数据。2.2 数据格式时间、风向、风速、浓度一个都不能少openair虽然好用但它对数据的“约定”很严格。核心要求就一句话数据必须是一个data.frame时间列要么叫date要么你能在调用函数时指定风向风速列的名称也要能对上。如果你直接拿原始监测导出的Excel数据丢进windRose()大概率会得到“object ws not found”这类报错。所以数据规范化的第一步是统一列名。以我自己的习惯为例一个最基本的分析数据集长这样列名含义示例date时间戳POSIXct或可转换的字符2024-06-01 00:00:00ws风速单位m/s2.3wd风向单位度0-360135no2NO2浓度单位µg/m342.5o3O3浓度单位µg/m388.1pm2.5PM2.5浓度单位µg/m335.2如果原始数据列名不是这套最简单的方法是读入R之后用rename()改掉library(dplyr) my_data - read.csv(站点数据.csv, stringsAsFactors FALSE) my_data - my_data %% rename(date 时间, ws 风速, wd 风向, no2 NO2浓度, o3 O3浓度, pm25 PM2.5_浓度)改完列名后第二个关键步骤是把时间列转成日期时间格式。很多人导入CSV后时间列是字符型肉眼看着正常但openair内部会强制转换并报错。建议统一这样处理my_data$date - as.POSIXct(my_data$date, format %Y-%m-%d %H:%M:%S, tz UTC)这里提个醒如果你的数据本来就是北京时间时区可以写tz Asia/Shanghai但不要一会儿UTC一会儿本地时间混着用尤其是做时间序列分析时时区不一致会导致按日聚合的结果出现偏移。第三个要点是异常值和缺失值的预处理。风向的合法范围是0到360度风速理论上≥0。如果原始数据里有-9999、9999这类哨兵值最好是先置成NA而不是留着让openair计算。还有重复时间戳建议先去重否则做时间聚合时报错或者结果被重复计算。缺失值这块openair最常用的windRose()会自动忽略NA但在做timePlot()时如果缺测太多会断线这时可以用interpolate()函数或自己填充。我自己整理数据的习惯流程是读CSV → 改列名 → 转时间 → 去重复 → 把哨兵值置NA → 再head()和summary()检查一遍。整个过程不超过五分钟但能帮你省掉后面两小时的调试时间。3. 核心绘图函数实操从玫瑰图到时间序列3.1 风玫瑰图与污染玫瑰图windRose和pollutionRose当数据准备好以后第一个值得画的图就是风玫瑰图。openair里的语法简洁得让人感动windRose(mydata, ws ws, wd wd, ws.int 2)上面这行代码会输出一张经典的风玫瑰图八个或更多方向的扇区每个扇区又被不同颜色分层颜色代表风速区间扇区半径长度代表该方向风出现的频率。参数ws.int 2表示风速按2 m/s为一个区间划分如果你不设置openair会自动选择间隔。对新手来说默认设置已经够用。但光看风玫瑰还不够环境分析里更常画的是污染玫瑰图pollutionRose。它跟风玫瑰的区别在于颜色和半径表达的不只是“风从哪来”而是“从某个方向吹来的风携带了多少污染物”。举个例子pollutionRose(mydata, pollutant no2, breaks c(0, 20, 40, 80, 160), col c(#4575B4, #74ADD1, #FEE090, #F46D43, #A50026))这张图一出来你马上能看出高浓度NO2主要来自哪个方向。很多时候污染玫瑰图比纯风玫瑰图更有用因为它直接反映了污染源方向。如果你发现高浓度来自某个工业园区方向那下一步就是核对那个方向上的排放源清单了。pollutionRose还有一个非常好用的参数statistic cpf。CPF全称是条件概率函数Conditional Probability Function它计算的是“当风向来自某个方向时浓度超过某个阈值的概率”。阈值默认是数据的75%分位数你可以用percentile 90改成90%分位。这个指标在识别异常高浓度事件的来向上特别有效。跑CPF的代码很直白pollutionRose(mydata, pollutant pm25, statistic cpf, percentile 90)图上扇区越大、颜色越深说明这个方向超过90%分位数浓度的概率越高。这种图在空气质量分析报告里是很有说服力的证据。顺带说一个我踩过的坑breaks参数的设置会直接改变图形信息量。如果区间设得太宽低浓度细节会被吞掉太窄又会导致颜色太碎。一般我会先用默认值跑一遍再根据浓度分布范围手动调整4到5个色阶区间兼顾美观和可读性。3.2 时间序列图timePlot快速定位污染过程玫瑰图解决的是“空间方向”问题时间序列图解决的则是“时间过程”问题。openair的timePlot()函数比R基础绘图或ggplot2更顺手的地方在于它可以同时绘制多个变量自动处理时间轴、图例和标题还支持添加平滑曲线和均值线。最基础用法timePlot(mydata, pollutant c(no2, o3), smooth TRUE)这样会得到一张包含NO2和O3两条时间曲线的图。我在这类图里最喜欢加的参数是smooth TRUE它会叠加上一条移动平均平滑曲线特别适合用来观察整体趋势而不是被逐小时的锯齿细节干扰。如果要观察特定时间段的污染过程可以用date.format或xlim来控制范围。比如去年12月有几天重污染过程我想看这七天的逐小时变化timePlot(subset(mydata, date as.POSIXct(2024-12-01) date as.POSIXct(2024-12-07)), pollutant c(pm25, no2))实际分析中timePlot最常见的一个问题是不同污染物浓度量级差异太大。比如NO2浓度在0-200 µg/m3而CO浓度可能在0-2 mg/m3如果把CO和NO2画在同一张图CO会几乎变成一条贴在底部的直线。解决办法是分图timePlot(mydata, pollutant no2) timePlot(mydata, pollutant co)或者用group TRUE把变量分面展示。openair里的分组逻辑是通过type参数实现的比如type season可以按季节拆分成四个小面板这个我在后面章节详细讲。3.3 时间变化图timeVariation看周期规律找排放特征如果说timePlot是看“哪一天出了事”那timeVariation就是看“这类事是不是经常在某个时间段发生”。这个函数把数据按小时、星期、月份三个维度拆开算平均值和置信区间一次性输出三张图timeVariation(mydata, pollutant no2)默认输出三列面板第一列是“hour of day”的24小时日变化曲线第二列是“weekday”的星期变化周一至周日第三列是“month”的逐月变化。这玩意儿看起来一目了然。我一般会先看日变化。如果一个站点的NO2日变化呈现“早晚两个高峰”那基本可以判断交通源贡献明显如果是夜间反而升高则可能是夜间残留、边界层降低或者局地排放如果O3日变化是“午后单峰”那就是典型的光化学生成特征。这些规律用timeVariation一次就能确认。如果你还想进一步对比不同站点或不同季节的差异可以加参数timeVariation(mydata, pollutant c(no2, o3), type season)这时候图会变成四行春夏秋冬每行三列时、周、月能很直观地看出不同季节污染规律的差异。比如夏天O3的午后峰值特别明显冬天NO2的早晚峰可能更突出。这种图写论文、做汇报都是加分项。timeVariation还有个优点自动生成置信区间让你知道曲线到底是信号还是噪声。如果某个时段曲线波动的置信区间很宽那可能样本量太少下结论要谨慎。这一点是很多自己用ggplot2画的均值曲线不具备的。4. 进阶分析功能趋势、关联与模型验证4.1 趋势检验TheilSen数据量小也不慌分析气象数据早晚会碰到一个灵魂拷问“这个污染物浓度到底是在升还是在降”很多人的第一反应是做个线性回归看斜率。但空气质量数据往往不满足正态分布还有季节性波动和异常值普通最小二乘回归很容易被个别极端值带偏。openair专门封装了TheilSen趋势检验这是一种基于中位数斜率的非参数回归方法在含有多个异常值和偏态分布的时间序列里非常稳健。用法如下theil - TheilSen(mydata, pollutant no2, ylab NO2 (ug/m3)) plot(theil)输出结果会包含趋势斜率、置信区间和P值。比如你看到斜率是-1.2 µg/m3/yearP值小于0.05那就能下结论说NO2呈显著下降趋势。注意TheilSen()默认是基于年度趋势的如果你的数据只有一年它会尝试计算月均趋势斜率解读的时候要看清res输出里的slope单位。我自己的经验是用TheilSen之前最好先对数据进行季节性分解或至少利用timeVariation理解季节特征否则你可能看到一个总体显著下降但其实只是春季数据多、冬季数据少导致的假象。openair里有一个配套函数叫Trend可以按站点、按季节分别输出趋势结果适合多站点对比Trend(mydata, pollutant pm25, type site)不过Trend函数的具体参数在不同版本里有些调整用之前最好?Trend看一眼帮助文档。4.2 散点图、泰勒图测点校准和模式评估的利器做环境监测的人经常会遇到两类对比需求一是把两台仪器放在一起观测看数据可比不可比二是把模型模拟值和实测值放一起评估模型准不准。这两种需求都可以用openair里的两个函数轻松解决。先看scatterPlot()它本质上就是散点图加拟合线但在细节上做了很多优化。它支持按季节、按站点分组着色还可以加入不同形式的拟合线比如线性、二次多项式甚至平滑样条scatterPlot(mydata, x no2, y pm25, modline list(linear TRUE))如果你把x轴设为观测值、y轴设为模拟值那这图就是一个快速评估模型偏差的工具。图上如果数据点都压在1:1线附近说明模型很好如果系统性偏离可以加一条拟合线看截距和斜率。openair还顺手把相关系数、均方根误差这些指标输出在控制台省去了单独计算的麻烦。泰勒图TaylorDiagram则是更专业的模型评估工具它把相关系数、标准差和中心均方根误差画在同一个极坐标图上。调用方式是这样# 假设你的数据里有obs观测和mod模拟两列 TaylorDiagram(mydata, obs obs, mod mod, group season)一张泰勒图可以同时展示春夏秋冬四个季节的模型表现点越靠近参考点代表模拟越好。模型开发、论文审稿阶段这张图几乎是“标配级”证据图。我见过不少同行用其他工具折腾半天画泰勒图其实openair几行代码就能搞定。5. 常见问题与排查技巧实录5.1 典型报错与解决思路对照表格直接查在实际使用中openair的报错信息有时候比较简略新手容易一脸懵。我把这几年来遇到过的典型报错整理成一个速查表方便你直接对号入座。报错信息常见原因解决办法object ws not found函数没找到指定列检查数据框列名或用ws 你的风速列名显式指定could not find function windRoseopenair包没加载运行library(openair)invalid type argumenttype参数写错检查type取值是否在允许集合内如season、year、weekdayError in as.POSIXct...时间列格式不对用as.POSIXct()显式转换并指定formatIn get.var(..., default) : No data to plot筛选后没有数据检查时间范围、列名或NA是否过多图形输出成空白PDF绘图设备驱动问题先dev.off()再用png()或直接保存为高分辨率图片这里特别想展开说一个最容易忽略的问题中文字符。如果数据框的列名包含中文即使你改参数指定了ws 风速在部分R版本或系统字体环境下openair内部的匹配逻辑可能出问题导致图例或变量识别失败。我的习惯是所有列名统一改成英文图中标题、标签用main、xlab参数手动加展示中文靠后面后期或文本但内部处理绝不掺中文。5.2 实操经验总结我踩过的坑和总结出来的习惯动作前面讲的都是具体函数最后这部分我聊聊自己在长期用openair时沉淀下来的几条习惯应该能帮你少踩不少坑。第一拿到数据先做三个检查时间格式、列名、缺失值。这个习惯我反复强调因为openair所有的函数都对数据格式有严格要求。你花五分钟做预处理后续所有分析丝滑顺畅反过来你急着画图大概率会卡在某个报错信息上反复调试。第二风向和风速尽量保留原始观测精度。有些原始数据里风向是整数、风速是一位小数这就很好。但你有时会遇到已经做过处理的数据比如风向被四舍五入到45度整倍数风速被分成等级这时候画出来的玫瑰图会非常“棱角分明”扇区边界很突兀。我的做法是能用原始数据就不用处理过的数据如果实在只有分级数据那就把breaks和扇区分辨率调低避免过度解释。第三画图前先跑一段“标准流水线”。我自己在实际操作中的体会是openair最适合的工作流是这样的用timePlot()快速检查数据质量看看有没有明显的野值或断测再用windRose()和pollutionRose()理解气象条件与浓度分布的方向关系然后用timeVariation()确认日、周、月尺度的周期规律最后如果需要报告或论文里的结论再用TheilSen()和scatterPlot()做趋势和相关性检验。这套流程基本能覆盖多数气象数据分析需求而且每一步的代码量都很小。第四把图形保存成矢量格式。如果你用png()保存风玫瑰图分辨率偏低时图例文字会发虚放大后没法看。我一般建议用pdf()或者RStudio的Export功能保存成PDF矢量格式这样无论插到论文里还是放大到展板上清晰度都有保障。若实在需要png()至少设置width 2400, height 1800, res 300。第五多看帮助文档。openair的帮助文档写得很详细每个参数都有详细说明和示例。别觉得英文看不下去配合一个在线翻译插件你能从帮助文档里挖出很多小技巧。比如type weekday可以直接按工作日/休息日分组这个参数我认识的大部分同行都没用过但在分析城市交通污染时特别实用。最后再分享一个小技巧openair包可以和shiny配合做成交互式分析小工具。我自己曾经把windRose()、timePlot()和timeVariation()封装进一个简单的shiny应用让不会R的同事也能自己选定站点和时间段点击按钮就出图。对团队协作和汇报来说这种交互式图表的感染力远强于静态图片。等你这篇文章里的函数都跑熟了完全可以往这个方向尝试一下。