1. 内容整体设计与思路拆解1.1 为什么偏偏是openair包气象数据分析这个事在R语言生态里其实有点尴尬。base R自带的基础绘图确实能用ggplot2接手以后可视化质量也有了质的飞跃但真到了气象环境数据分析这个具体场景你会发现一个很现实的问题数据预处理太琐碎了。风向要转成16方位风速要按国标分级污染物浓度要和气象要素做相关性分析还要做风玫瑰图、污染玫瑰图、时间变化趋势、浓度分位数分析……这些功能如果全靠自己一行行写代码代码量会膨胀到让人崩溃而且不同项目的代码风格还不统一换个数据集就得重新调。openair包就是专门来解决这个痛点的。它是英国Kings College London团队开发的开源R包最初是为了支持欧洲空气质量管理研究后来功能不断扩展逐渐覆盖了气象数据分析和环境数据分析的主流场景。它的核心设计思路是“按分析任务封装成函数”比如你想画风玫瑰图直接调windRose()想分析污染物浓度随时间的变化趋势直接调timePlot()或timeVariation()想看看不同风向条件下某个污染物的浓度分布污染玫瑰图一屏搞定。这个设计思路的实际价值在于它把“气象和环境数据探索性分析”这一整条工作流标准化了。你不需要每次从零开始写数据清洗、统计汇总、图表定制的代码只要数据格式符合要求调用对应函数就能拿到标准化的分析结果。对于科研工作者、环保监测站的分析人员、环评工程师来说这套工具能省掉的重复劳动是相当可观的。1.2 它能解决什么问题从我自己的使用体验来看openair包解决的核心问题可以归纳成四类。第一类是风向风速的可视化。这是气象分析里最基础也最频繁的需求风玫瑰图几乎是每个气象分析报告的标配。openair的windRose()函数不仅支持经典的16方位风玫瑰还能按污染物浓度分档、按时间段分面、按季节拆分信息量比传统风玫瑰大得多。第二类是污染物浓度与气象条件的关系分析。比如你要回答“这个测点的PM2.5浓度在什么风向下最高”这类问题用pollutionRose()函数直接就能出图不用自己写风向分箱、浓度统计、极坐标绘制的代码。第三类是时间序列分析。气象和环境监测数据天然是时间序列openair的timePlot()可以做多变量时间序列对比timeVariation()可以拆分出日变化、周变化、月变化甚至季节变化对识别数据规律特别有用。第四类是数据质量检验和基础统计汇总。openair内置了summaryPlot()、timeAverage()等数据概览和聚合工具能快速检查数据完整度、异常值、缺失情况。所以这个包适合谁呢主要三类人做大气环境研究的科研人员、环保系统里做监测数据分析的技术人员、以及环评或咨询公司里需要处理气象和污染物数据的工程师。当然如果你的专业背景是气象学主要做温度、气压、降水这类纯气象要素分析openair也能覆盖一大部分需求只是它的强项更偏向“气象条件与空气质量结合”的场景。1.3 项目分析的整体框架这篇文章我计划按照一个完整的数据分析项目流程来展开从环境准备到数据导入再到核心分析和图形解读最后是常见问题的排查基本上覆盖了你拿到一套气象数据之后从零开始做分析的完整路径。我会结合一个实际案例来演示假设你有某个城市监测站点一年的逐小时气象数据包含风速、风向、温度、湿度、气压和污染物浓度数据PM2.5、PM10、NO2、O3目标是分析这个站点污染物的季节变化特征、日变化规律以及气象条件对污染浓度的影响。这套流程走下来你对openair包的核心功能就有一个完整的认识后续换成自己的数据时核心代码基本可以直接复用。2. 数据准备与预处理openair的前置条件2.1 安装和加载openair包的安装方式和其他CRAN包一样一条命令搞定。不过有一点要提醒openair的依赖包比较多包括lattice、latticeExtra、ggplot2部分版本、cluster、dplyr等第一次安装时如果网络不好或者R版本太旧可能会遇到依赖包安装失败的问题。建议在安装前先把R更新到当前主流版本确保依赖能正常编译或下载。install.packages(openair)如果是国内网络环境CRAN镜像可能会很慢甚至超时失败。我一般会换成国内镜像源安装速度会快很多# 设置清华镜像源速度更稳定 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) install.packages(openair)加载的时候没什么特殊的和普通包一样library(openair)加载之后有一个小习惯建议你培养用sessionInfo()或packageVersion(openair)确认一下版本号。因为这个包更新频率不低不同版本之间API可能会有细微调整你在网上搜索到的教程代码如果运行报错“Error: could not find function”先检查一下版本是不是太旧。2.2 数据格式要求和类型转换openair对数据格式有一套统一的要求这个必须先讲清楚因为据我观察大约有一半的新手问题都出在数据格式上。核心要求所有与时间相关的数据必须有一列名为date的POSIXct类型列。其他变量列名可以自己定义但数据类型要一致数值列必须是numeric不能用character。如果数据里时间是分开的年月日时分几个字段需要先把它们合并成一列再转换为POSIXct类型。来看一个典型的实际操作场景。假设你从监测平台导出一个Excel文件里面有年、月、日、时、风速、风向、温度、PM2.5这些列导入R后会变成这样# 假设已经用 readxl 或 read.csv 读入数据变量名为 raw_data # 检查数据结构 str(raw_data)如果看到date相关字段是character或者年月日分散在多列就需要手动拼接转换# 将年月日时合并为一个时间列 raw_data$date - as.POSIXct( paste(raw_data$year, raw_data$month, raw_data$day, raw_data$hour, sep -), format %Y-%m-%d-%H, tz Asia/Shanghai )如果要处理的是小时数据还可以用更简洁的方式直接从原始时间戳转换# 如果原始数据已经有一个时间字符串列例如 2024-01-01 00:00:00 raw_data$date - as.POSIXct(raw_data$时间, format %Y-%m-%d %H:%M:%S, tz Asia/Shanghai)时区这里要特别注意。国内监测站的数据通常都是北京时间东八区如果你的R会话默认时区是UTC时间轴会整体偏8小时画出来的日变化图会明显错位。建议统一显式指定tz Asia/Shanghai不要依赖系统默认值。如果你的数据本身就是UTC时间比如某些再分析资料那就保持UTC但后续分析时要知道这一点。转换完成之后还要做两件事。第一是确保date列按时间顺序排列openair部分函数对乱序数据虽然不会报错但时间序列图会画得很乱第二是检查重复时间戳尤其是多个监测点数据合并之后可能会出现同一时间多条记录的情况这种数据如果不处理timeAverage()等聚合函数会返回意料之外的结果。2.3 数据概览与质量初检数据准备好之后不要急着画图分析先用openair自带的工具做一轮快速体检。summaryPlot()是openair提供的一个数据概览函数用法很简单summaryPlot(raw_data)这个函数会画出一张综合图上面是各变量的时间序列迷你图下面是数据缺失率的统计和变量分布直方图。通过这张图你能快速掌握几个关键信息哪些变量有缺失值、缺失发生在什么时间段、是否有明显异常值、各变量的取值范围是否合理。举个例子如果你发现某个月的风向数据大量缺失而这个月恰好又有一轮重污染过程那后续分析污染玫瑰图时就要谨慎解释结果。再比如如果PM2.5出现负数那明显是仪器零漂或校准造成的异常值需要根据业务逻辑决定剔除还是置为NA。这个环节的核心思路是先让数据问题暴露出来再决定后续清洗策略。不要一上来就用完整数据跑分析因为异常值对统计结果的影响远比你想象的大。尤其在做浓度分位数、均值计算这类对数值敏感的统计时一个离谱的负值就能让整个月均值失真。数据清洗我一般按这个顺序做# 1. 剔除风速为负值、风向不在0-360度的记录 raw_data$ws[raw_data$ws 0] - NA raw_data$wd[raw_data$wd 0 | raw_data$wd 360] - NA # 2. 剔除超出合理物理范围的极端值按污染物特性设定阈值 raw_data$pm2.5[raw_data$pm2.5 0] - NA raw_data$pm2.5[raw_data$pm2.5 1000] - NA # 3. 剔除重复时间戳 raw_data - raw_data[!duplicated(raw_data$date), ]有两点提醒。第一负浓度值在国标里通常按无效数据处理直接置NA比用插值更稳妥除非你有明确的质控说明可以解释负值的来源。第二如果是做趋势分析、年际比较这类需要连续数据的工作置NA之后可能会导致大量空白时间段这时候再考虑时间插值方法比如线性插值或kriging但插值方法的选择要写进方法学里不能默默处理。3. 核心分析实操从风玫瑰到污染玫瑰3.1 风玫瑰图读懂一个测点的风气候特征风玫瑰图是气象分析的基础图件openair的windRose()函数是我用过所有绘图工具里做风玫瑰最顺手的。基础用法一句话windRose(raw_data, ws ws, wd wd)这里ws和wd是风速和风向列的名称按你自己的数据列名调整就行。默认情况下函数会按标准风速分档比如0-0.5、0.5-2、2-4、4-6、6-8、8-10、10 m/s和16方位N、NNE、NE等汇总中间圆点表示静风比例。但别小看这个函数它的几个高级参数在实际项目中非常关键。第一个是ws.int参数用于调整风速分档间隔。默认值是1但不同城市的风速分布差异很大平原城市可能经常出现大风山区城市则常年静风这时候按默认间隔画出来的图可能全是同一个色阶区分度很差。我一般会先看风速的summary统计再确定合理的分档间隔。第二个是pollutant参数。这是windRose最有价值的功能之一它可以在风玫瑰的每个风向上叠加显示该风向区间内的污染物平均浓度用颜色深浅表示浓度高低。这个图在识别污染来源方向时极其直观。windRose(raw_data, ws ws, wd wd, pollutant pm2.5)第三个是type参数用于按因子分面。比如说你想看看这个站点在双休日和工作日的风玫瑰是否不同或者在不同季节的风玫瑰差异直接用type season或type weekday就能自动拆分windRose(raw_data, ws ws, wd wd, type season)这个功能在实际分析中非常有用。我做过一个案例某个站点的全年风玫瑰显示主导风向是东南风但split by season之后发现冬季实际上是以北风为主东南风只在夏季占主导。如果只看全年风玫瑰很容易得出误导性结论。风玫瑰图的输出格式是lattice图形对象所以可以用openair自带的绘图控制参数来调整标题、图例位置等。也可以直接用print()把它输出为png、pdf文件png(windrose_season.png, width 10, height 8, units in, res 300) windRose(raw_data, ws ws, wd wd, type season) dev.off()这里有个小坑很多人在Windows系统下用png()函数输出中文字体时会出现乱码或方块但windRose默认的文本基本都是英文所以问题不大。如果后续你添加自定义标题时用了中文建议提前用windowsFonts()设置中文字体或者直接输出矢量pdf再用AI软件编辑。3.2 污染玫瑰图污染物浓度的风向响应pollutionRose()是windRose的孪生函数它的侧重点从风的分布变成了污染物浓度的风向分布。通俗地说这张图回答的问题是当风从某一个方向吹来的时候污染物浓度是高还是低这个图在污染溯源分析中几乎是必备工具。比如某城市的监测站点位于工业园区下风向那么当刮特定方向的风时污染物浓度会显著升高。污染玫瑰图能把这种关系可视化得非常清楚。基本用法pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd)默认情况下图上的每个风向扇区会显示该风向下污染物浓度的分位数箱线图或者说类似“风向-浓度”的统计箱线。这种表现方式比单纯看平均值要可靠得多因为平均浓度很容易被少数极端值拉高而分位数箱线图则能体现整个浓度分布状况。如果数据量足够还可以用statistic mean或statistic max来切换显示不同统计量pollutionRose(raw_data, pollutant o3, ws ws, wd wd, statistic mean)对于臭氧这类受光化学反应控制明显的二次污染物用mean统计量通常能反映出清洁风向下浓度低、停滞条件下浓度高的特征。除了单污染物分析pollutionRose还支持按浓度区间叠加多种污染物的联合显示。比如把PM2.5和NO2的浓度同时画在一张污染玫瑰上通过不同颜色标识两种污染物就能快速判断它们是不是同源。这对于做污染源解析的初步研判很有价值。不过这里要提醒一点污染玫瑰图只讨论相关性不能直接证明因果。某个风向高浓度可能只是因为这个风向上游有污染源也可能是因为这个风向对应特定的气象条件比如静稳天气导致污染物累积。所以我在实际项目里一般会配套做风向-风速-浓度的三元关系分析或者进一步用条件概率函数CPF来量化不同风向上浓度超标的条件概率。openair这个方向其实还有更进阶的函数比如polarPlot()它用极坐标的方式展示风速、风向下污染物浓度的插值分布场。polarPlot比污染玫瑰提供的信息密度更高可以同时看出风向、风速和浓度三者之间的关系是识别近距离污染源扩散形态的好工具polarPlot(raw_data, pollutant pm2.5, ws ws, wd wd)这张图画出来的效果是一个彩色极坐标等值面横轴代表风速从中心向外增大角度代表风向颜色代表浓度。如果图像中存在一个清晰的高浓度“斑块”集中在某个风向和风速区间通常暗示那个方向存在局地排放源。3.3 时间变化模式分析日变化、周变化与季节变化气象和环境数据的核心特征之一就是时间周期性。timeVariation()函数是openair里我最常推荐给数据分析师的功能之一因为它用一张图就展示了多个周期尺度的变化规律。用法timeVariation(raw_data, pollutant pm2.5)函数会默认画出四个面板日变化hour of day、周变化day of week、月变化month和整体时间序列。每张面板都显示均值线和基于自助法计算的置信区间方便你判断波动是否统计显著。日变化图的价值在于识别排放源类型和气象过程。比如PM2.5的日变化如果呈现出早晚双峰型一般与交通排放早晚高峰有关如果是夜间单峰型则可能是夜间边界层降低导致污染物在地面累积。O3的日变化则刚好相反午后浓度达到峰值夜间由于滴定消耗降到最低这是光化学反应主导的典型特征。季节变化图能反映污染的季节性规律。华北地区的PM2.5冬高夏低O3则是春夏高、秋冬低这些特征在同一张图上对照观察会非常直观。type参数在这里也同样奏效。比如你想对比工作日和周末的日变化差异可以这样timeVariation(raw_data, pollutant pm2.5, type weekday)输出图会按工作日/周末分面展示日变化曲线。这在研究“周末效应”即污染物浓度的周内波动时非常有用。还有一个小技巧timeVariation不仅支持数值型污染物还可以把气象要素如温度、风速等作为pollutant传入这样就能同时观察气象条件和污染物浓度的同步变化关系。举个例子如果你把pm2.5和ws风速同时传入timeVariation(raw_data, pollutant c(pm2.5, ws))会画出多个变量的日变化对比图。你可能会发现风速在午后增强而PM2.5在午后降低两者呈现负相关关系这为解释污染浓度日变化提供了气象动力学的背景。3.4 时间序列绘图与平滑趋势分析时间序列图是气象分析里最基础也最直接的图形表达。openair的timePlot()函数虽然看起来只是把基础绘图封装了一层但它的几个特性让我在实际工作中离不开了。首先是多变量对比功能。你可以把PM2.5、PM10、NO2等污染物和一两个气象变量放在同一个时间轴里对比用不同颜色区分不需要自己写多面板拼图代码timePlot(raw_data, pollutant c(pm2.5, pm10, no2), group TRUE)group TRUE会把所有变量画在同一个面板里适合量纲一致的浓度对比。如果变量量纲差异太大比如温度和风速一起画建议用group FALSE让每个变量单独一个小面板但共享同一个时间轴这样变季节趋势时非常清楚。其次timePlot支持平滑曲线叠加。通过smooth参数打开loess平滑后原始序列的噪声会被过滤长趋势一目了然timePlot(raw_data, pollutant pm2.5, smooth TRUE)这个功能很适合快速判断一段时间内污染物浓度是上升还是下降趋势。比如分析“某城市过去三年PM2.5年均浓度是否真的在下降”用timePlot加平滑曲线先看一眼再进入正式的统计检验思路会清晰很多。timePlot还有一个常用于成果出图的功能支持布局和样式设置。你可以用lattice图形的维护参数调整x轴时间跨度、y轴范围、配色等。如果你需要出发表级别的图建议输出矢量格式后配合Adobe Illustrator微调。3.5 时间聚合与统计摘要数据聚合是处理气象数据的日常操作。原始监测数据往往是分钟级或小时级但分析报告通常需要日均值、月均值或年均值。openair的timeAverage()函数能把这个过程压缩成一行代码# 将小时数据聚合成日均值 daily_data - timeAverage(raw_data, avg.time day) # 聚合到月均值 monthly_data - timeAverage(raw_data, avg.time month)这个函数的优势在于它对日期时间格式的处理非常稳健它能自动识别时间步长按calendar time聚合同时还支持处理缺失值比例、数据覆盖度等参数。如果某一天的数据缺失率达到50%以上按日均值计算时应该判定为无效日openair里可以通过min.freq参数控制# 只有当一天内有超过75%的小时数据时才计算日均值 daily_data - timeAverage(raw_data, avg.time day, min.freq 0.75)这个细节在中国环境空气质量标准里是有对应要求的有效日均值至少需要一定的有效小时数否则该日数据无效。如果手动写这个逻辑代码会啰嗦且容易出错用timeAverage一行搞定。除了按时间聚合openair还提供timeSelect()函数可以用它对时间子集做筛选分析。比如只看冬季12月-次年2月的数据winter_data - timeSelect(raw_data, start 2024-12-01, end 2025-02-28)不过我更常用的是结合dplyr的filter函数来筛选时间范围因为灵活度更高library(dplyr) winter_data - raw_data %% filter(date 2024-12-01 date 2025-03-01)3.6 相关性分析和散点图矩阵气象要素和污染物浓度之间往往存在复杂的相关关系。比如温度和O3通常正相关风速和PM2.5通常负相关。openair的scatterPlot()函数提供了带回归线的散点图绘制功能并且在分面支持下可以做出非常丰富的对比图。基本用法scatterPlot(raw_data, x ws, y pm2.5, linear TRUE)linear TRUE会在图上叠加线性回归线并显示R²和回归方程。这对快速评估两个变量的相关性强度很有帮助。如果想看不同季节下风速与PM2.5的关系差异加type参数scatterPlot(raw_data, x ws, y pm2.5, type season, linear TRUE)你可能会看到夏季和冬季的回归斜率明显不同冬季风速增大对PM2.5的稀释效果更强而夏季可能因为二次生成机制主导风速与PM2.5的关系没那么紧密。需要特别提醒的是相关性不代表因果气象要素和污染物之间的关系往往受多重因素协同影响散点图只是探索性分析的第一步。更严谨的做法是在看完散点图后引入多元回归、广义加性模型GAM或机器学习模型做进一步的定量分析。3.7 其他几个高频函数openair还有一些使用频率也很高的函数简单提一下方便你在实际项目中按需选用。trendLevel()可以画出污染浓度随时间x轴和另一变量一般是小时或月份y轴变化的二维热力图。比如trendLevel(raw_data, pollutant o3, x day, y hour)可以直观展示O3浓度在一天中的不同时段和全年不同日期的分布情况能快速识别光化学污染的季节高发时段。calendarPlot()把污染物浓度按日历形式展示出来每天一格颜色表示日均浓度等级。这种图在月度报告里非常直观领导或业主一眼就能看出哪段时间污染重calendarPlot(raw_data, pollutant pm2.5, year 2024)smoothTrend()则是用多种方法LOESS、Theil-Sen等估算时间序列的趋势项适合做长期趋势分析smoothTrend(raw_data, pollutant pm2.5)4. 实操过程与核心环节实现一个完整的分析案例4.1 案例数据说明下面我以一个具体的分析案例把前面讲的内容串起来。假设你手头有一份2024年全年逐小时监测数据包含以下变量变量名含义单位date时间POSIXctws风速m/swd风向度0-360temp温度℃rh相对湿度%pm2.5PM2.5浓度μg/m³pm10PM10浓度μg/m³no2NO2浓度μg/m³o3O3浓度μg/m³这份数据是我根据真实监测数据的典型特征模拟的不是实际监测站的原始数据但数据结构和特征规律符合实际情况完全可以用这个流程套自己的数据。4.2 完整流程代码# 加载包 library(openair) library(dplyr) # 从CSV读入数据实际使用中改成你的路径 raw_data - read.csv(airquality_2024.csv, stringsAsFactors FALSE) # 检查数据结构 str(raw_data) # 时间列转换 raw_data$date - as.POSIXct(raw_data$date, format %Y-%m-%d %H:%M:%S, tz Asia/Shanghai) # 数据质量初检 summaryPlot(raw_data) # 过滤异常值 raw_data$ws[raw_data$ws 0] - NA raw_data$wd[raw_data$wd 0 | raw_data$wd 360] - NA raw_data$pm2.5[raw_data$pm2.5 0] - NA raw_data$pm2.5[raw_data$pm2.5 800] - NA # 去除重复时间戳 raw_data - raw_data[!duplicated(raw_data$date), ] # 生成季节和星期因子openair部分函数需要 raw_data$month - format(raw_data$date, %m) # 1. 全年风玫瑰图 png(fig01_windrose_year.png, width 8, height 7, units in, res 300) windRose(raw_data, ws ws, wd wd) dev.off() # 2. 按季节风玫瑰图 png(fig02_windrose_season.png, width 12, height 10, units in, res 300) windRose(raw_data, ws ws, wd wd, type season) dev.off() # 3. PM2.5污染玫瑰图 png(fig03_pollutionrose_pm25.png, width 8, height 7, units in, res 300) pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd) dev.off() # 4. 时间变化模式分析 png(fig04_timevariation_pm25.png, width 12, height 9, units in, res 300) timeVariation(raw_data, pollutant pm2.5) dev.off() # 5. 多污染物时间序列 png(fig05_timeseries_all.png, width 14, height 8, units in, res 300) timePlot(raw_data, pollutant c(pm2.5, pm10, no2, o3)) dev.off() # 6. 风速与PM2.5散点关系 png(fig06_scatter_ws_pm25.png, width 8, height 7, units in, res 300) scatterPlot(raw_data, x ws, y pm2.5, linear TRUE) dev.off() # 7. 按月聚合做趋势分析 monthly_data - timeAverage(raw_data, avg.time month) print(monthly_data) # 8. 日均值数据导出 daily_data - timeAverage(raw_data, avg.time day, min.freq 0.75) write.csv(daily_data, daily_data_2024.csv, row.names FALSE)4.3 关键结果解读与实际分析这套代码跑完之后你会得到一组图和数据表。我自己在实际执行时会重点关注几个分析结论这里分享一个示例性的解读思路。从按季节风玫瑰图来看如果冬季主导风向是北风或西北风且对应风速区间以2-4 m/s为主说明这个站点在冬季主要受偏北气流影响冷空气活动频繁时可能带来上游区域的污染传输。如果夏季主导风向转为东南风且风速相对较小那么本地排放加上静稳条件就可能是夏季污染的主要成因。再看PM2.5污染玫瑰图如果污染浓度高值集中在偏北和偏西方向结合风玫瑰图判断这些方向是否有工业区或交通干道就能初步形成一个污染来源假说。这时候配合polarPlot做进一步验证观察浓度高值区是否集中在特定风速区间可以区分是本地排放通常高浓度出现在低风速区因为扩散条件差还是区域传输高浓度出现在中等风速且风向明确的条件。timeVariation图的解读是最有意思的部分。如果PM2.5日变化曲线呈早高峰7-9时和晚高峰18-21时双峰特征且周变化在工作日偏高、周末偏低大概率交通源占主导。如果夜间浓度比白天还高并且季节变化呈现冬季明显偏高那就要考虑采暖排放源和夜间边界层下降的共同作用。O3的日变化如果呈现下午峰值说明光化学反应条件良好如果O3在夜间不降反升可能是残余层向下混合或者区域传输的影响这种情况在复杂地形区域并不少见。4.4 导出和报告集成分析完成后把核心图表导出成统一风格再嵌入到报告里是实际项目中最常做的事。我通常会这样做把需要用于报告的图统一设置尺寸为16:9或4:3、300 dpi这样放进Word或PPT里清晰度足够。一些重要的统计结果比如月均值、季节均值、超标天数、风向风速与浓度的相关性系数用write.csv导出为表格方便报告写作时直接引用。如果你需要把多个图拼在一张大图上lattice的gridExtra包可以用但我更推荐用openair自身的output lattice接口再结合grid.arrange来做library(gridExtra) p1 - windRose(raw_data, ws ws, wd wd, type season) p2 - pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd) grid.arrange(p1, p2, ncol 2)5. 常见问题与排查技巧实录5.1 “date列不存在”或时间格式错误这个问题的出现频率高得惊人。openair几乎所有函数都要求数据框里有一个名为date的POSIXct列如果列名写成Date、time或日期变量类型是Date函数会直接报错Error in checkPrep(...) : Need a date column or something similar.排查思路很简单先看str(data)确认列名和类型。如果列名不对改列名names(data)[names(data) time] - date如果类型不是POSIXct先用as.POSIXct()转换。还有一点要注意有些用户从Excel导入的数据里日期时间列会被读成numeric尤其是Excel里时间显示成数字的格式这种情况千万不要直接as.POSIXct强制转换先要理解这个数字代表什么——它可能是距离某个基准时间的天数或秒数直接转换会得到完全错误的时间。5.2 时间序列图形错位八小时这个问题几乎每个用中国数据的用户都会遇到。处理办法其实在数据准备部分已经说过显式指定时区tz Asia/Shanghai。如果转换完之后发现图形还是错位的一种可能是你的原始数据里时间是以字符串存储的并且已经包含了时区偏移信息在as.POSIXct()时又重复指定了时区。我的建议是先用attr()查看转换后的date列是否有tz属性确认没有UTC残留。attr(data$date, tzone)如果是空值或UTC但你想用北京时间重新设置attr(data$date, tzone) - Asia/Shanghai另外提醒R的日期时间格式里%z表示时区偏移%Z表示时区缩写容易混淆解析时别用错。5.3 风玫瑰图全部为空或只有静风一个环节有时候画出来的风玫瑰图几乎全是中间的一个圆盘外围没有风向分布看起来像数据没进来。这种情况大概率是风向变量里全是NA值或者风向范围异常。先用summary(data$wd)检查风向变量的统计量。如果min为0、max为360但大量集中在某个值比如0那可能是传感器故障或数据编码问题。另外注意部分气象数据里静风用风向0或999表示0度静风和0度正北风容易被混淆需要结合业务逻辑清洗。如果风向数据本身没问题检查一下是否有缺失值影响绘图sum(is.na(data$wd))如果某个季节或月份的风向缺失率过高可以考虑按邻近时间插值但插值的合理性要在分析方法说明中交代清楚。5.4 污染玫瑰图上浓度颜色分层不明显默认的浓度color scale是等间隔分箱但如果你的数据分布极度偏斜比如绝大多数浓度都在50以下突然有几个500的强污染过程色阶会被少数极高值拉伸导致低浓度区的颜色区分度很低。解决办法调整breaks参数手动设置颜色断点pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd, breaks c(0, 25, 50, 75, 100, 150, 200, 300, 500))也可以把浓度做对数变换后再绘图但这样图例标签的解释成本会增加非学术场合不太推荐。我更常用的是手动断点结合国内空气质量分指数AQI的浓度区间来设置这样图和国标对照起来更直观。5.5 timeAverage聚合结果比预期少了或多了一天timeAverage默认按calendar day聚合但不同时区的日界可能不一样。比如北京时间日界是0点但如果数据里含有UTC时间的观测记录聚合结果就会和你预期的日期差一天。解决办法确保所有数据时间已经统一为同一个时区并且在聚合前用range(data$date)确认时间范围正确。另外如果数据有分钟级或秒级的时间戳timeAverage的avg.time参数还支持1 hour、30 min等细分设置多检查一下参数名别拼错。5.6 图形中文乱码虽然openair默认的图形文本以英文为主但你自己标标题、坐标轴时可能会用中文。在Windows下画图时会显示成方块这是因为R默认的中文字体映射不对。解决办法在输出png前设置字体# Windows下指定中文字体 windowsFonts(yahei windowsFont(Microsoft YaHei))或者在lattice的par.settings里指定字体族。如果你是在Linux服务器上批量出图建议直接用英文标题省心也避免字体依赖。5.7 性能问题数据量太大跑得慢如果数据量是分钟级的多站点多年数据openair绘图函数可能变得很慢。此时建议按站点和时间范围拆分数据分批处理。另外可以先按小时聚合减少数据量后再绘图hourly_data - timeAverage(raw_data, avg.time hour)如果还不够快把不参与分析的列删掉只保留date、ws、wd和要分析的污染物列能明显提升速度。5.8 一个容易踩的潜在坑列名冲突如果你自己构建数据框的时候用了ws或wd以外的列名比如风速列叫wind_speed风向列叫wind_dir那在调用windRose时就必须显式指定windRose(raw_data, ws wind_speed, wd wind_dir)如果列名与函数默认参数不一致又没显式指定函数会按默认的wswd去找列然后报错。这个小问题在一些自己拼数据的场景里挺常见的写代码时留意一下就行。6. 进阶应用思路openair还能干这些事6.1 把openair和其他包组合使用openair的很多图形对象都是lattice对象这意味着你可以用lattice生态的工具做进一步定制。同时它导出的数据比如timeAverage的结果、timeVariation的统计输出可以直接交给dplyr和ggplot2做更个性化的处理和可视化。在实际项目中我经常这样组合使用数据清洗用dplyr完成基础探索性图表用openair快速完成需要合成为汇报材料时把关键统计结果输出成data.frame用ggplot2重新精修成风格统一的图表。openair负责高效探索ggplot2负责最终呈现两者配合效率很高。6.2 多站点对比分析如果你手头有多个站点的监测数据openair同样适用。只要保证所有站点数据的列名一致用dplyr的bind_rows()合并数据然后使用type site参数或者手动构造一个site因子列就能分站点出图。# 假设两个站点数据分别为 site_A, site_B all_sites - bind_rows(site_A %% mutate(site A), site_B %% mutate(site B)) windRose(all_sites, ws ws, wd wd, type site) timeVariation(all_sites, pollutant pm2.5, type site)这种多站点对比图在分析区域污染空间差异时特别直观。比如两个站点一个在城区、一个在郊区对比它们的PM2.5日变化曲线城区可能呈现明显的交通双峰郊区则相对平缓这对理解不同功能区的污染特征很有帮助。6.3 与气象再分析资料结合openair不仅能分析站点观测数据也能处理网格化的再分析数据如ERA5前提是先将网格数据提取到目标站点位置整理成标准的时间序列数据框。提取之后类似的分析流程完全一致。我就做过一个案例把ERA5的10米风速和2米温度提取到某城市的监测站点位置然后和地面观测的PM2.5浓度一起分析气象条件对污染过程的影响。再分析数据的一个优势是变量齐全、时间连续对于观测缺测时期的数据补充很有价值。6.4 污染过程案例分析如果你关注的是某一次重污染过程的演变openair的timePlot配合timeVariation足够用了。更精细的做法是锁定一个污染过程的时间窗口利用timeSelect()截取数据然后对比污染过程前中后三个阶段的污染玫瑰图变化process_before - timeSelect(raw_data, start 2024-01-05, end 2024-01-09) process_during - timeSelect(raw_data, start 2024-01-10, end 2024-01-14) process_after - timeSelect(raw_data, start 2024-01-15, end 2024-01-19) pollutionRose(process_before, pollutant pm2.5, ws ws, wd wd) pollutionRose(process_during, pollutant pm2.5, ws ws, wd wd) pollutionRose(process_after, pollutant pm2.5, ws ws, wd wd)通过对比三个阶段的污染玫瑰图和风场变化往往能判断这次重污染以本地积累为主还是区域传输为主这在写污染过程分析报告时是很关键的内容。7. 写在最后的经验总结用openair这套包分析气象数据最大的感受是它真正做到了“把时间留给思考而不是写代码”。气象数据的探索性分析其实有大量重复模式——画风玫瑰、看时间变化、分析浓度和气象的关系、做统计摘要这些事如果每次都用base R或ggplot2从零实现不仅代码量大、容易出差错而且由于各个项目中绘图参数的差异最后出图的风格也很不统一。openair把这些模式固化成了标准化的函数反而给了分析人员更多时间去思考数据本身的科学含义。但我也想强调openair不是万能的。它适合做探索性分析和基础图件产出如果你的分析需要深入到复杂的统计建模、机器学习预测或更精细的空间插值还是需要配合其他R包来完成。把它定位成“气象环境数据分析的瑞士军刀”比较合适——万能工具不现实但日常八成以上的图形探索需求它都能解决。最后再分享一个实用习惯每次拿到一套新数据我都会先跑一遍summaryPlot()和windRose()哪怕最终报告里根本不用这两张图。因为这两张图能在三分钟内让我对数据质量、风场特征和浓度水平有个整体印象后面的分析思路会清晰很多。数据分析和写作一样动笔之前先搭好框架后面每个环节都会顺很多。 R语言openair包：气象与空气质量数据可视化分析实战
1. 内容整体设计与思路拆解1.1 为什么偏偏是openair包气象数据分析这个事在R语言生态里其实有点尴尬。base R自带的基础绘图确实能用ggplot2接手以后可视化质量也有了质的飞跃但真到了气象环境数据分析这个具体场景你会发现一个很现实的问题数据预处理太琐碎了。风向要转成16方位风速要按国标分级污染物浓度要和气象要素做相关性分析还要做风玫瑰图、污染玫瑰图、时间变化趋势、浓度分位数分析……这些功能如果全靠自己一行行写代码代码量会膨胀到让人崩溃而且不同项目的代码风格还不统一换个数据集就得重新调。openair包就是专门来解决这个痛点的。它是英国Kings College London团队开发的开源R包最初是为了支持欧洲空气质量管理研究后来功能不断扩展逐渐覆盖了气象数据分析和环境数据分析的主流场景。它的核心设计思路是“按分析任务封装成函数”比如你想画风玫瑰图直接调windRose()想分析污染物浓度随时间的变化趋势直接调timePlot()或timeVariation()想看看不同风向条件下某个污染物的浓度分布污染玫瑰图一屏搞定。这个设计思路的实际价值在于它把“气象和环境数据探索性分析”这一整条工作流标准化了。你不需要每次从零开始写数据清洗、统计汇总、图表定制的代码只要数据格式符合要求调用对应函数就能拿到标准化的分析结果。对于科研工作者、环保监测站的分析人员、环评工程师来说这套工具能省掉的重复劳动是相当可观的。1.2 它能解决什么问题从我自己的使用体验来看openair包解决的核心问题可以归纳成四类。第一类是风向风速的可视化。这是气象分析里最基础也最频繁的需求风玫瑰图几乎是每个气象分析报告的标配。openair的windRose()函数不仅支持经典的16方位风玫瑰还能按污染物浓度分档、按时间段分面、按季节拆分信息量比传统风玫瑰大得多。第二类是污染物浓度与气象条件的关系分析。比如你要回答“这个测点的PM2.5浓度在什么风向下最高”这类问题用pollutionRose()函数直接就能出图不用自己写风向分箱、浓度统计、极坐标绘制的代码。第三类是时间序列分析。气象和环境监测数据天然是时间序列openair的timePlot()可以做多变量时间序列对比timeVariation()可以拆分出日变化、周变化、月变化甚至季节变化对识别数据规律特别有用。第四类是数据质量检验和基础统计汇总。openair内置了summaryPlot()、timeAverage()等数据概览和聚合工具能快速检查数据完整度、异常值、缺失情况。所以这个包适合谁呢主要三类人做大气环境研究的科研人员、环保系统里做监测数据分析的技术人员、以及环评或咨询公司里需要处理气象和污染物数据的工程师。当然如果你的专业背景是气象学主要做温度、气压、降水这类纯气象要素分析openair也能覆盖一大部分需求只是它的强项更偏向“气象条件与空气质量结合”的场景。1.3 项目分析的整体框架这篇文章我计划按照一个完整的数据分析项目流程来展开从环境准备到数据导入再到核心分析和图形解读最后是常见问题的排查基本上覆盖了你拿到一套气象数据之后从零开始做分析的完整路径。我会结合一个实际案例来演示假设你有某个城市监测站点一年的逐小时气象数据包含风速、风向、温度、湿度、气压和污染物浓度数据PM2.5、PM10、NO2、O3目标是分析这个站点污染物的季节变化特征、日变化规律以及气象条件对污染浓度的影响。这套流程走下来你对openair包的核心功能就有一个完整的认识后续换成自己的数据时核心代码基本可以直接复用。2. 数据准备与预处理openair的前置条件2.1 安装和加载openair包的安装方式和其他CRAN包一样一条命令搞定。不过有一点要提醒openair的依赖包比较多包括lattice、latticeExtra、ggplot2部分版本、cluster、dplyr等第一次安装时如果网络不好或者R版本太旧可能会遇到依赖包安装失败的问题。建议在安装前先把R更新到当前主流版本确保依赖能正常编译或下载。install.packages(openair)如果是国内网络环境CRAN镜像可能会很慢甚至超时失败。我一般会换成国内镜像源安装速度会快很多# 设置清华镜像源速度更稳定 options(repos c(CRAN https://mirrors.tuna.tsinghua.edu.cn/CRAN/)) install.packages(openair)加载的时候没什么特殊的和普通包一样library(openair)加载之后有一个小习惯建议你培养用sessionInfo()或packageVersion(openair)确认一下版本号。因为这个包更新频率不低不同版本之间API可能会有细微调整你在网上搜索到的教程代码如果运行报错“Error: could not find function”先检查一下版本是不是太旧。2.2 数据格式要求和类型转换openair对数据格式有一套统一的要求这个必须先讲清楚因为据我观察大约有一半的新手问题都出在数据格式上。核心要求所有与时间相关的数据必须有一列名为date的POSIXct类型列。其他变量列名可以自己定义但数据类型要一致数值列必须是numeric不能用character。如果数据里时间是分开的年月日时分几个字段需要先把它们合并成一列再转换为POSIXct类型。来看一个典型的实际操作场景。假设你从监测平台导出一个Excel文件里面有年、月、日、时、风速、风向、温度、PM2.5这些列导入R后会变成这样# 假设已经用 readxl 或 read.csv 读入数据变量名为 raw_data # 检查数据结构 str(raw_data)如果看到date相关字段是character或者年月日分散在多列就需要手动拼接转换# 将年月日时合并为一个时间列 raw_data$date - as.POSIXct( paste(raw_data$year, raw_data$month, raw_data$day, raw_data$hour, sep -), format %Y-%m-%d-%H, tz Asia/Shanghai )如果要处理的是小时数据还可以用更简洁的方式直接从原始时间戳转换# 如果原始数据已经有一个时间字符串列例如 2024-01-01 00:00:00 raw_data$date - as.POSIXct(raw_data$时间, format %Y-%m-%d %H:%M:%S, tz Asia/Shanghai)时区这里要特别注意。国内监测站的数据通常都是北京时间东八区如果你的R会话默认时区是UTC时间轴会整体偏8小时画出来的日变化图会明显错位。建议统一显式指定tz Asia/Shanghai不要依赖系统默认值。如果你的数据本身就是UTC时间比如某些再分析资料那就保持UTC但后续分析时要知道这一点。转换完成之后还要做两件事。第一是确保date列按时间顺序排列openair部分函数对乱序数据虽然不会报错但时间序列图会画得很乱第二是检查重复时间戳尤其是多个监测点数据合并之后可能会出现同一时间多条记录的情况这种数据如果不处理timeAverage()等聚合函数会返回意料之外的结果。2.3 数据概览与质量初检数据准备好之后不要急着画图分析先用openair自带的工具做一轮快速体检。summaryPlot()是openair提供的一个数据概览函数用法很简单summaryPlot(raw_data)这个函数会画出一张综合图上面是各变量的时间序列迷你图下面是数据缺失率的统计和变量分布直方图。通过这张图你能快速掌握几个关键信息哪些变量有缺失值、缺失发生在什么时间段、是否有明显异常值、各变量的取值范围是否合理。举个例子如果你发现某个月的风向数据大量缺失而这个月恰好又有一轮重污染过程那后续分析污染玫瑰图时就要谨慎解释结果。再比如如果PM2.5出现负数那明显是仪器零漂或校准造成的异常值需要根据业务逻辑决定剔除还是置为NA。这个环节的核心思路是先让数据问题暴露出来再决定后续清洗策略。不要一上来就用完整数据跑分析因为异常值对统计结果的影响远比你想象的大。尤其在做浓度分位数、均值计算这类对数值敏感的统计时一个离谱的负值就能让整个月均值失真。数据清洗我一般按这个顺序做# 1. 剔除风速为负值、风向不在0-360度的记录 raw_data$ws[raw_data$ws 0] - NA raw_data$wd[raw_data$wd 0 | raw_data$wd 360] - NA # 2. 剔除超出合理物理范围的极端值按污染物特性设定阈值 raw_data$pm2.5[raw_data$pm2.5 0] - NA raw_data$pm2.5[raw_data$pm2.5 1000] - NA # 3. 剔除重复时间戳 raw_data - raw_data[!duplicated(raw_data$date), ]有两点提醒。第一负浓度值在国标里通常按无效数据处理直接置NA比用插值更稳妥除非你有明确的质控说明可以解释负值的来源。第二如果是做趋势分析、年际比较这类需要连续数据的工作置NA之后可能会导致大量空白时间段这时候再考虑时间插值方法比如线性插值或kriging但插值方法的选择要写进方法学里不能默默处理。3. 核心分析实操从风玫瑰到污染玫瑰3.1 风玫瑰图读懂一个测点的风气候特征风玫瑰图是气象分析的基础图件openair的windRose()函数是我用过所有绘图工具里做风玫瑰最顺手的。基础用法一句话windRose(raw_data, ws ws, wd wd)这里ws和wd是风速和风向列的名称按你自己的数据列名调整就行。默认情况下函数会按标准风速分档比如0-0.5、0.5-2、2-4、4-6、6-8、8-10、10 m/s和16方位N、NNE、NE等汇总中间圆点表示静风比例。但别小看这个函数它的几个高级参数在实际项目中非常关键。第一个是ws.int参数用于调整风速分档间隔。默认值是1但不同城市的风速分布差异很大平原城市可能经常出现大风山区城市则常年静风这时候按默认间隔画出来的图可能全是同一个色阶区分度很差。我一般会先看风速的summary统计再确定合理的分档间隔。第二个是pollutant参数。这是windRose最有价值的功能之一它可以在风玫瑰的每个风向上叠加显示该风向区间内的污染物平均浓度用颜色深浅表示浓度高低。这个图在识别污染来源方向时极其直观。windRose(raw_data, ws ws, wd wd, pollutant pm2.5)第三个是type参数用于按因子分面。比如说你想看看这个站点在双休日和工作日的风玫瑰是否不同或者在不同季节的风玫瑰差异直接用type season或type weekday就能自动拆分windRose(raw_data, ws ws, wd wd, type season)这个功能在实际分析中非常有用。我做过一个案例某个站点的全年风玫瑰显示主导风向是东南风但split by season之后发现冬季实际上是以北风为主东南风只在夏季占主导。如果只看全年风玫瑰很容易得出误导性结论。风玫瑰图的输出格式是lattice图形对象所以可以用openair自带的绘图控制参数来调整标题、图例位置等。也可以直接用print()把它输出为png、pdf文件png(windrose_season.png, width 10, height 8, units in, res 300) windRose(raw_data, ws ws, wd wd, type season) dev.off()这里有个小坑很多人在Windows系统下用png()函数输出中文字体时会出现乱码或方块但windRose默认的文本基本都是英文所以问题不大。如果后续你添加自定义标题时用了中文建议提前用windowsFonts()设置中文字体或者直接输出矢量pdf再用AI软件编辑。3.2 污染玫瑰图污染物浓度的风向响应pollutionRose()是windRose的孪生函数它的侧重点从风的分布变成了污染物浓度的风向分布。通俗地说这张图回答的问题是当风从某一个方向吹来的时候污染物浓度是高还是低这个图在污染溯源分析中几乎是必备工具。比如某城市的监测站点位于工业园区下风向那么当刮特定方向的风时污染物浓度会显著升高。污染玫瑰图能把这种关系可视化得非常清楚。基本用法pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd)默认情况下图上的每个风向扇区会显示该风向下污染物浓度的分位数箱线图或者说类似“风向-浓度”的统计箱线。这种表现方式比单纯看平均值要可靠得多因为平均浓度很容易被少数极端值拉高而分位数箱线图则能体现整个浓度分布状况。如果数据量足够还可以用statistic mean或statistic max来切换显示不同统计量pollutionRose(raw_data, pollutant o3, ws ws, wd wd, statistic mean)对于臭氧这类受光化学反应控制明显的二次污染物用mean统计量通常能反映出清洁风向下浓度低、停滞条件下浓度高的特征。除了单污染物分析pollutionRose还支持按浓度区间叠加多种污染物的联合显示。比如把PM2.5和NO2的浓度同时画在一张污染玫瑰上通过不同颜色标识两种污染物就能快速判断它们是不是同源。这对于做污染源解析的初步研判很有价值。不过这里要提醒一点污染玫瑰图只讨论相关性不能直接证明因果。某个风向高浓度可能只是因为这个风向上游有污染源也可能是因为这个风向对应特定的气象条件比如静稳天气导致污染物累积。所以我在实际项目里一般会配套做风向-风速-浓度的三元关系分析或者进一步用条件概率函数CPF来量化不同风向上浓度超标的条件概率。openair这个方向其实还有更进阶的函数比如polarPlot()它用极坐标的方式展示风速、风向下污染物浓度的插值分布场。polarPlot比污染玫瑰提供的信息密度更高可以同时看出风向、风速和浓度三者之间的关系是识别近距离污染源扩散形态的好工具polarPlot(raw_data, pollutant pm2.5, ws ws, wd wd)这张图画出来的效果是一个彩色极坐标等值面横轴代表风速从中心向外增大角度代表风向颜色代表浓度。如果图像中存在一个清晰的高浓度“斑块”集中在某个风向和风速区间通常暗示那个方向存在局地排放源。3.3 时间变化模式分析日变化、周变化与季节变化气象和环境数据的核心特征之一就是时间周期性。timeVariation()函数是openair里我最常推荐给数据分析师的功能之一因为它用一张图就展示了多个周期尺度的变化规律。用法timeVariation(raw_data, pollutant pm2.5)函数会默认画出四个面板日变化hour of day、周变化day of week、月变化month和整体时间序列。每张面板都显示均值线和基于自助法计算的置信区间方便你判断波动是否统计显著。日变化图的价值在于识别排放源类型和气象过程。比如PM2.5的日变化如果呈现出早晚双峰型一般与交通排放早晚高峰有关如果是夜间单峰型则可能是夜间边界层降低导致污染物在地面累积。O3的日变化则刚好相反午后浓度达到峰值夜间由于滴定消耗降到最低这是光化学反应主导的典型特征。季节变化图能反映污染的季节性规律。华北地区的PM2.5冬高夏低O3则是春夏高、秋冬低这些特征在同一张图上对照观察会非常直观。type参数在这里也同样奏效。比如你想对比工作日和周末的日变化差异可以这样timeVariation(raw_data, pollutant pm2.5, type weekday)输出图会按工作日/周末分面展示日变化曲线。这在研究“周末效应”即污染物浓度的周内波动时非常有用。还有一个小技巧timeVariation不仅支持数值型污染物还可以把气象要素如温度、风速等作为pollutant传入这样就能同时观察气象条件和污染物浓度的同步变化关系。举个例子如果你把pm2.5和ws风速同时传入timeVariation(raw_data, pollutant c(pm2.5, ws))会画出多个变量的日变化对比图。你可能会发现风速在午后增强而PM2.5在午后降低两者呈现负相关关系这为解释污染浓度日变化提供了气象动力学的背景。3.4 时间序列绘图与平滑趋势分析时间序列图是气象分析里最基础也最直接的图形表达。openair的timePlot()函数虽然看起来只是把基础绘图封装了一层但它的几个特性让我在实际工作中离不开了。首先是多变量对比功能。你可以把PM2.5、PM10、NO2等污染物和一两个气象变量放在同一个时间轴里对比用不同颜色区分不需要自己写多面板拼图代码timePlot(raw_data, pollutant c(pm2.5, pm10, no2), group TRUE)group TRUE会把所有变量画在同一个面板里适合量纲一致的浓度对比。如果变量量纲差异太大比如温度和风速一起画建议用group FALSE让每个变量单独一个小面板但共享同一个时间轴这样变季节趋势时非常清楚。其次timePlot支持平滑曲线叠加。通过smooth参数打开loess平滑后原始序列的噪声会被过滤长趋势一目了然timePlot(raw_data, pollutant pm2.5, smooth TRUE)这个功能很适合快速判断一段时间内污染物浓度是上升还是下降趋势。比如分析“某城市过去三年PM2.5年均浓度是否真的在下降”用timePlot加平滑曲线先看一眼再进入正式的统计检验思路会清晰很多。timePlot还有一个常用于成果出图的功能支持布局和样式设置。你可以用lattice图形的维护参数调整x轴时间跨度、y轴范围、配色等。如果你需要出发表级别的图建议输出矢量格式后配合Adobe Illustrator微调。3.5 时间聚合与统计摘要数据聚合是处理气象数据的日常操作。原始监测数据往往是分钟级或小时级但分析报告通常需要日均值、月均值或年均值。openair的timeAverage()函数能把这个过程压缩成一行代码# 将小时数据聚合成日均值 daily_data - timeAverage(raw_data, avg.time day) # 聚合到月均值 monthly_data - timeAverage(raw_data, avg.time month)这个函数的优势在于它对日期时间格式的处理非常稳健它能自动识别时间步长按calendar time聚合同时还支持处理缺失值比例、数据覆盖度等参数。如果某一天的数据缺失率达到50%以上按日均值计算时应该判定为无效日openair里可以通过min.freq参数控制# 只有当一天内有超过75%的小时数据时才计算日均值 daily_data - timeAverage(raw_data, avg.time day, min.freq 0.75)这个细节在中国环境空气质量标准里是有对应要求的有效日均值至少需要一定的有效小时数否则该日数据无效。如果手动写这个逻辑代码会啰嗦且容易出错用timeAverage一行搞定。除了按时间聚合openair还提供timeSelect()函数可以用它对时间子集做筛选分析。比如只看冬季12月-次年2月的数据winter_data - timeSelect(raw_data, start 2024-12-01, end 2025-02-28)不过我更常用的是结合dplyr的filter函数来筛选时间范围因为灵活度更高library(dplyr) winter_data - raw_data %% filter(date 2024-12-01 date 2025-03-01)3.6 相关性分析和散点图矩阵气象要素和污染物浓度之间往往存在复杂的相关关系。比如温度和O3通常正相关风速和PM2.5通常负相关。openair的scatterPlot()函数提供了带回归线的散点图绘制功能并且在分面支持下可以做出非常丰富的对比图。基本用法scatterPlot(raw_data, x ws, y pm2.5, linear TRUE)linear TRUE会在图上叠加线性回归线并显示R²和回归方程。这对快速评估两个变量的相关性强度很有帮助。如果想看不同季节下风速与PM2.5的关系差异加type参数scatterPlot(raw_data, x ws, y pm2.5, type season, linear TRUE)你可能会看到夏季和冬季的回归斜率明显不同冬季风速增大对PM2.5的稀释效果更强而夏季可能因为二次生成机制主导风速与PM2.5的关系没那么紧密。需要特别提醒的是相关性不代表因果气象要素和污染物之间的关系往往受多重因素协同影响散点图只是探索性分析的第一步。更严谨的做法是在看完散点图后引入多元回归、广义加性模型GAM或机器学习模型做进一步的定量分析。3.7 其他几个高频函数openair还有一些使用频率也很高的函数简单提一下方便你在实际项目中按需选用。trendLevel()可以画出污染浓度随时间x轴和另一变量一般是小时或月份y轴变化的二维热力图。比如trendLevel(raw_data, pollutant o3, x day, y hour)可以直观展示O3浓度在一天中的不同时段和全年不同日期的分布情况能快速识别光化学污染的季节高发时段。calendarPlot()把污染物浓度按日历形式展示出来每天一格颜色表示日均浓度等级。这种图在月度报告里非常直观领导或业主一眼就能看出哪段时间污染重calendarPlot(raw_data, pollutant pm2.5, year 2024)smoothTrend()则是用多种方法LOESS、Theil-Sen等估算时间序列的趋势项适合做长期趋势分析smoothTrend(raw_data, pollutant pm2.5)4. 实操过程与核心环节实现一个完整的分析案例4.1 案例数据说明下面我以一个具体的分析案例把前面讲的内容串起来。假设你手头有一份2024年全年逐小时监测数据包含以下变量变量名含义单位date时间POSIXctws风速m/swd风向度0-360temp温度℃rh相对湿度%pm2.5PM2.5浓度μg/m³pm10PM10浓度μg/m³no2NO2浓度μg/m³o3O3浓度μg/m³这份数据是我根据真实监测数据的典型特征模拟的不是实际监测站的原始数据但数据结构和特征规律符合实际情况完全可以用这个流程套自己的数据。4.2 完整流程代码# 加载包 library(openair) library(dplyr) # 从CSV读入数据实际使用中改成你的路径 raw_data - read.csv(airquality_2024.csv, stringsAsFactors FALSE) # 检查数据结构 str(raw_data) # 时间列转换 raw_data$date - as.POSIXct(raw_data$date, format %Y-%m-%d %H:%M:%S, tz Asia/Shanghai) # 数据质量初检 summaryPlot(raw_data) # 过滤异常值 raw_data$ws[raw_data$ws 0] - NA raw_data$wd[raw_data$wd 0 | raw_data$wd 360] - NA raw_data$pm2.5[raw_data$pm2.5 0] - NA raw_data$pm2.5[raw_data$pm2.5 800] - NA # 去除重复时间戳 raw_data - raw_data[!duplicated(raw_data$date), ] # 生成季节和星期因子openair部分函数需要 raw_data$month - format(raw_data$date, %m) # 1. 全年风玫瑰图 png(fig01_windrose_year.png, width 8, height 7, units in, res 300) windRose(raw_data, ws ws, wd wd) dev.off() # 2. 按季节风玫瑰图 png(fig02_windrose_season.png, width 12, height 10, units in, res 300) windRose(raw_data, ws ws, wd wd, type season) dev.off() # 3. PM2.5污染玫瑰图 png(fig03_pollutionrose_pm25.png, width 8, height 7, units in, res 300) pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd) dev.off() # 4. 时间变化模式分析 png(fig04_timevariation_pm25.png, width 12, height 9, units in, res 300) timeVariation(raw_data, pollutant pm2.5) dev.off() # 5. 多污染物时间序列 png(fig05_timeseries_all.png, width 14, height 8, units in, res 300) timePlot(raw_data, pollutant c(pm2.5, pm10, no2, o3)) dev.off() # 6. 风速与PM2.5散点关系 png(fig06_scatter_ws_pm25.png, width 8, height 7, units in, res 300) scatterPlot(raw_data, x ws, y pm2.5, linear TRUE) dev.off() # 7. 按月聚合做趋势分析 monthly_data - timeAverage(raw_data, avg.time month) print(monthly_data) # 8. 日均值数据导出 daily_data - timeAverage(raw_data, avg.time day, min.freq 0.75) write.csv(daily_data, daily_data_2024.csv, row.names FALSE)4.3 关键结果解读与实际分析这套代码跑完之后你会得到一组图和数据表。我自己在实际执行时会重点关注几个分析结论这里分享一个示例性的解读思路。从按季节风玫瑰图来看如果冬季主导风向是北风或西北风且对应风速区间以2-4 m/s为主说明这个站点在冬季主要受偏北气流影响冷空气活动频繁时可能带来上游区域的污染传输。如果夏季主导风向转为东南风且风速相对较小那么本地排放加上静稳条件就可能是夏季污染的主要成因。再看PM2.5污染玫瑰图如果污染浓度高值集中在偏北和偏西方向结合风玫瑰图判断这些方向是否有工业区或交通干道就能初步形成一个污染来源假说。这时候配合polarPlot做进一步验证观察浓度高值区是否集中在特定风速区间可以区分是本地排放通常高浓度出现在低风速区因为扩散条件差还是区域传输高浓度出现在中等风速且风向明确的条件。timeVariation图的解读是最有意思的部分。如果PM2.5日变化曲线呈早高峰7-9时和晚高峰18-21时双峰特征且周变化在工作日偏高、周末偏低大概率交通源占主导。如果夜间浓度比白天还高并且季节变化呈现冬季明显偏高那就要考虑采暖排放源和夜间边界层下降的共同作用。O3的日变化如果呈现下午峰值说明光化学反应条件良好如果O3在夜间不降反升可能是残余层向下混合或者区域传输的影响这种情况在复杂地形区域并不少见。4.4 导出和报告集成分析完成后把核心图表导出成统一风格再嵌入到报告里是实际项目中最常做的事。我通常会这样做把需要用于报告的图统一设置尺寸为16:9或4:3、300 dpi这样放进Word或PPT里清晰度足够。一些重要的统计结果比如月均值、季节均值、超标天数、风向风速与浓度的相关性系数用write.csv导出为表格方便报告写作时直接引用。如果你需要把多个图拼在一张大图上lattice的gridExtra包可以用但我更推荐用openair自身的output lattice接口再结合grid.arrange来做library(gridExtra) p1 - windRose(raw_data, ws ws, wd wd, type season) p2 - pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd) grid.arrange(p1, p2, ncol 2)5. 常见问题与排查技巧实录5.1 “date列不存在”或时间格式错误这个问题的出现频率高得惊人。openair几乎所有函数都要求数据框里有一个名为date的POSIXct列如果列名写成Date、time或日期变量类型是Date函数会直接报错Error in checkPrep(...) : Need a date column or something similar.排查思路很简单先看str(data)确认列名和类型。如果列名不对改列名names(data)[names(data) time] - date如果类型不是POSIXct先用as.POSIXct()转换。还有一点要注意有些用户从Excel导入的数据里日期时间列会被读成numeric尤其是Excel里时间显示成数字的格式这种情况千万不要直接as.POSIXct强制转换先要理解这个数字代表什么——它可能是距离某个基准时间的天数或秒数直接转换会得到完全错误的时间。5.2 时间序列图形错位八小时这个问题几乎每个用中国数据的用户都会遇到。处理办法其实在数据准备部分已经说过显式指定时区tz Asia/Shanghai。如果转换完之后发现图形还是错位的一种可能是你的原始数据里时间是以字符串存储的并且已经包含了时区偏移信息在as.POSIXct()时又重复指定了时区。我的建议是先用attr()查看转换后的date列是否有tz属性确认没有UTC残留。attr(data$date, tzone)如果是空值或UTC但你想用北京时间重新设置attr(data$date, tzone) - Asia/Shanghai另外提醒R的日期时间格式里%z表示时区偏移%Z表示时区缩写容易混淆解析时别用错。5.3 风玫瑰图全部为空或只有静风一个环节有时候画出来的风玫瑰图几乎全是中间的一个圆盘外围没有风向分布看起来像数据没进来。这种情况大概率是风向变量里全是NA值或者风向范围异常。先用summary(data$wd)检查风向变量的统计量。如果min为0、max为360但大量集中在某个值比如0那可能是传感器故障或数据编码问题。另外注意部分气象数据里静风用风向0或999表示0度静风和0度正北风容易被混淆需要结合业务逻辑清洗。如果风向数据本身没问题检查一下是否有缺失值影响绘图sum(is.na(data$wd))如果某个季节或月份的风向缺失率过高可以考虑按邻近时间插值但插值的合理性要在分析方法说明中交代清楚。5.4 污染玫瑰图上浓度颜色分层不明显默认的浓度color scale是等间隔分箱但如果你的数据分布极度偏斜比如绝大多数浓度都在50以下突然有几个500的强污染过程色阶会被少数极高值拉伸导致低浓度区的颜色区分度很低。解决办法调整breaks参数手动设置颜色断点pollutionRose(raw_data, pollutant pm2.5, ws ws, wd wd, breaks c(0, 25, 50, 75, 100, 150, 200, 300, 500))也可以把浓度做对数变换后再绘图但这样图例标签的解释成本会增加非学术场合不太推荐。我更常用的是手动断点结合国内空气质量分指数AQI的浓度区间来设置这样图和国标对照起来更直观。5.5 timeAverage聚合结果比预期少了或多了一天timeAverage默认按calendar day聚合但不同时区的日界可能不一样。比如北京时间日界是0点但如果数据里含有UTC时间的观测记录聚合结果就会和你预期的日期差一天。解决办法确保所有数据时间已经统一为同一个时区并且在聚合前用range(data$date)确认时间范围正确。另外如果数据有分钟级或秒级的时间戳timeAverage的avg.time参数还支持1 hour、30 min等细分设置多检查一下参数名别拼错。5.6 图形中文乱码虽然openair默认的图形文本以英文为主但你自己标标题、坐标轴时可能会用中文。在Windows下画图时会显示成方块这是因为R默认的中文字体映射不对。解决办法在输出png前设置字体# Windows下指定中文字体 windowsFonts(yahei windowsFont(Microsoft YaHei))或者在lattice的par.settings里指定字体族。如果你是在Linux服务器上批量出图建议直接用英文标题省心也避免字体依赖。5.7 性能问题数据量太大跑得慢如果数据量是分钟级的多站点多年数据openair绘图函数可能变得很慢。此时建议按站点和时间范围拆分数据分批处理。另外可以先按小时聚合减少数据量后再绘图hourly_data - timeAverage(raw_data, avg.time hour)如果还不够快把不参与分析的列删掉只保留date、ws、wd和要分析的污染物列能明显提升速度。5.8 一个容易踩的潜在坑列名冲突如果你自己构建数据框的时候用了ws或wd以外的列名比如风速列叫wind_speed风向列叫wind_dir那在调用windRose时就必须显式指定windRose(raw_data, ws wind_speed, wd wind_dir)如果列名与函数默认参数不一致又没显式指定函数会按默认的wswd去找列然后报错。这个小问题在一些自己拼数据的场景里挺常见的写代码时留意一下就行。6. 进阶应用思路openair还能干这些事6.1 把openair和其他包组合使用openair的很多图形对象都是lattice对象这意味着你可以用lattice生态的工具做进一步定制。同时它导出的数据比如timeAverage的结果、timeVariation的统计输出可以直接交给dplyr和ggplot2做更个性化的处理和可视化。在实际项目中我经常这样组合使用数据清洗用dplyr完成基础探索性图表用openair快速完成需要合成为汇报材料时把关键统计结果输出成data.frame用ggplot2重新精修成风格统一的图表。openair负责高效探索ggplot2负责最终呈现两者配合效率很高。6.2 多站点对比分析如果你手头有多个站点的监测数据openair同样适用。只要保证所有站点数据的列名一致用dplyr的bind_rows()合并数据然后使用type site参数或者手动构造一个site因子列就能分站点出图。# 假设两个站点数据分别为 site_A, site_B all_sites - bind_rows(site_A %% mutate(site A), site_B %% mutate(site B)) windRose(all_sites, ws ws, wd wd, type site) timeVariation(all_sites, pollutant pm2.5, type site)这种多站点对比图在分析区域污染空间差异时特别直观。比如两个站点一个在城区、一个在郊区对比它们的PM2.5日变化曲线城区可能呈现明显的交通双峰郊区则相对平缓这对理解不同功能区的污染特征很有帮助。6.3 与气象再分析资料结合openair不仅能分析站点观测数据也能处理网格化的再分析数据如ERA5前提是先将网格数据提取到目标站点位置整理成标准的时间序列数据框。提取之后类似的分析流程完全一致。我就做过一个案例把ERA5的10米风速和2米温度提取到某城市的监测站点位置然后和地面观测的PM2.5浓度一起分析气象条件对污染过程的影响。再分析数据的一个优势是变量齐全、时间连续对于观测缺测时期的数据补充很有价值。6.4 污染过程案例分析如果你关注的是某一次重污染过程的演变openair的timePlot配合timeVariation足够用了。更精细的做法是锁定一个污染过程的时间窗口利用timeSelect()截取数据然后对比污染过程前中后三个阶段的污染玫瑰图变化process_before - timeSelect(raw_data, start 2024-01-05, end 2024-01-09) process_during - timeSelect(raw_data, start 2024-01-10, end 2024-01-14) process_after - timeSelect(raw_data, start 2024-01-15, end 2024-01-19) pollutionRose(process_before, pollutant pm2.5, ws ws, wd wd) pollutionRose(process_during, pollutant pm2.5, ws ws, wd wd) pollutionRose(process_after, pollutant pm2.5, ws ws, wd wd)通过对比三个阶段的污染玫瑰图和风场变化往往能判断这次重污染以本地积累为主还是区域传输为主这在写污染过程分析报告时是很关键的内容。7. 写在最后的经验总结用openair这套包分析气象数据最大的感受是它真正做到了“把时间留给思考而不是写代码”。气象数据的探索性分析其实有大量重复模式——画风玫瑰、看时间变化、分析浓度和气象的关系、做统计摘要这些事如果每次都用base R或ggplot2从零实现不仅代码量大、容易出差错而且由于各个项目中绘图参数的差异最后出图的风格也很不统一。openair把这些模式固化成了标准化的函数反而给了分析人员更多时间去思考数据本身的科学含义。但我也想强调openair不是万能的。它适合做探索性分析和基础图件产出如果你的分析需要深入到复杂的统计建模、机器学习预测或更精细的空间插值还是需要配合其他R包来完成。把它定位成“气象环境数据分析的瑞士军刀”比较合适——万能工具不现实但日常八成以上的图形探索需求它都能解决。最后再分享一个实用习惯每次拿到一套新数据我都会先跑一遍summaryPlot()和windRose()哪怕最终报告里根本不用这两张图。因为这两张图能在三分钟内让我对数据质量、风场特征和浓度水平有个整体印象后面的分析思路会清晰很多。数据分析和写作一样动笔之前先搭好框架后面每个环节都会顺很多。