R语言及ggplot2在环境空气监测数据可视化中的应用探究

　　随着新《环境空气质量标准》(GB 3095-2012)的发布、实施，全国各地都在采用环境空气自动监测以应对其严苛的数据有效性要求。环境空气自动监测系统每天都在产生大量的数据，如何对数据进行有效的统计分析，成为一个新的课题。对数据的可视化，是数据分析的第一步。本文应用R语言及ggplot2图形绘制包以浙江省常山县城区两个空气自动站即环保大楼站及图书馆站2014年全年监测数据为例进行各种可视化分析，以探讨该软件在空气质量数据分析领域应用的巨大潜力。
　　1 软件准备
　　1.1 R语言软件与扩展包
　　R语言是一种区分大小写的解释性语言，其具有的强大统计计算及绘图能力，是从大数据中获取有用信息的绝佳工具，它提供了Windows、OS X、Linux等各大操作系统下的版本，可以直接从网上免费下载、安装、使用。R语言软件在基本安装中提供了大量的数据处理、统计和图形函数，此外各社区还开发了数以千计的扩展包(packages)为R增加了更多令人惊奇的功能，本文用到的ggplot2、plyr、reshape2等即是其中的一部分。
　　ggplot2是目前R语言中的图形绘制扩展包，它为R语言提供了一个基于全面而连贯的语法的绘图系统，它由一系列独立的图形部件组成，并能以许多种不同的方式组合起来，使得数据分析者能用紧凑的语法轻松绘制出复杂的图形，从而使数据分析者更能将注意力集中于数据分析本身。plyr与reshape2是两个功能强大的数据整理扩展包，它们为R语言扩展了对数据变换、数据变形的功能。
　　ggplot2软件包可以在R语言命令提示符后输入install.packages(ggplot2)，选择合适的下载镜像后，就可以自动从网上下载安装。在Windows操作系统下下载的是二进制文件，可以直接使用;在类Unix操作系统下，下载的是包的源代码，经过编译后才能使用，但只要具备gcc等编译工具，安装都可通过简单的命令完成。在使用前，需要输入library(ggplot2)或require(ggplot2)调用该软件包。plyr包等也是如此。
　　1.2 Excel
　　但面对原始数据，微软公司的Excel往往是对数据进行清洗、分列的直观分析利器。各地自动监测站采集并导出的原始数据往往是xml格式的，而Excel处理这种格式极为得心应手。Excel中的数据筛选功能可以轻易地删除大量无效数据，分列功能可以将监测日期与时刻分离以便于下一步的处理，数据格式也可以得到统一，将数据另存为csv格式后即可通过read.csv()函数导入R软件。R软件也可以通过write.csv()函数写出csv格式的文件与Excel交互使用。
　　2 数据处理
　　2.1 数据的收集与整理
　　在数据分析工作中，收集、整理数据的时间往往是占了工作时间的80%以上。在数据收集、整理工作中交互地使用R语言与Excel，可极大地提高工作效率。首先使用Excel将在自动监测过程中产生的无效数据进行筛选、剔除，再按需要，对某些数据进行分列处理，最后导出为csv文件后以便R语言读取。与Matlab等将一切视为矩阵不同，R语言可以灵活地对向量、矩阵、数组、数据框、因子、列表等多种数据结构进行处理。通过read.csv()导入的数据都被视为数据框(data.frame)。数据框是R语言中最常用的数据结构，其每一列都可以包含不同模式的数据，如环境空气监测数据之后，可以增加名为month的一列，用中文一月份表示该数据是一月份的数据，增加site的一列，用中文标注其站点名称，如假设hbdl201401为导入的环保大楼站一月份数据，用R命令行示例如下：site-rep(环保大楼站，dim(hbdl201401))month-rep(一月份，dim(hbdl201401))hbdl201401-cbind(hbdl201401，site，month)同一站点中的数据，可以使用rbind()函数进行按月合并。当多站点数据合并时，由于设备配置的原因，数据框中各列的名称与数量并不完全一致，使用rbind()合并数据时会出错。如常山县环保大楼站除常规六参数外还配有碳黑仪、能见度仪、气象五参数等，而图书馆站中仅有常规六参数。这时的合并可以用到merge()函数。假设hbdl2014为常山县环保大楼站2014年全年小时均值数据，而tsg2014为图书馆站数据，则可用以下命令合并数据框：
　　cs2014-merge(hbdl2014，tsg2014，all=T，sort=F)合并后的数据框cs2014中在图书馆站中所没有的碳黑仪等数据均以NA(缺失)表示。merge()函数可以两两合并大量来源不同的数据框，只需保证列名的统一，数据即可合并无误。在R语言中用命令tail(cs2014)即可列出cs2014的最后6行：
　　2.2 数据的统计变换
　　对数据分析的细化过程往往就是分组的过程。2.1节中为数据增加month一列即是为数据按月份进行分组，site一列即是为站点分组，clock一列即是为数据按时刻进行分组。在实际工作中，分析数据的角度不同，分组也是千奇百怪。如可按当日温度差来分组，也可按平均风速的大小来分组，当然也可以按PM2.5日均值大小来分组来考察其他观测值的情况。数据的变换与变形则是基于统计以及下一步分析的需要，将数据变化成相应的形式。
　　2.2.1 创建分组变量。在数据分析过程中，常常发现这种情况：组别分的太细或者需要根据原数值内容创建分组。
　　第一种情况可以合并组别。例如按气象学上的规定，以每年的三、四、五月为春季，六、七、八月为夏季，九、十、十一月为秋季，十二、一、二月为冬季。需要以季节分组时，可以调用如下命令：
　　c s 2 0 1 4 $ s e a s o n [ c s 2 0 1 4 $ m o n t h = = 三月份|cs2014$month==四月份|cs2014$month==五月份-春
　　c s 2 0 1 4 $ s e a s o n [ c s 2 0 1 4 $ m o n t h = = 六月份|cs2014$month==七月份|cs2014$month==八月份-夏
　　c s 2 0 1 4 $ s e a s o n [ c s 2 0 1 4 $ m o n t h = = 九月份|cs2014$month==十月份|cs2014$month==十一月份-秋
　　c s 2 0 1 4 $ s e a s o n [ c s 2 0 1 4 $ m o n t h = = 十二月份|cs2014$month==一月份|cs2014$month==二月份-冬
　　运行上述命令后，数据框cs2014中就多了一个名为season的列，其标注了每一行数据的季节。第二种情况可以用到R语言中的cut()函数。如原数据中各污染物按月份按时刻进行浓度平均得到数据框cs2014month后，可以以0 0.040mg/m3、0.040 0.080mg/m3、0.080mg/m3为区间将臭氧的时刻平均数据分作3个组别，分别名为弱、中、强，用cut()函数创建分组变量o3level：cs2014month$o3level-cut(cs2014month$o3，breaks=c( 0 ， 0 . 0 4 0 ， 0 . 0 8 0 ， I n f ) ， l a b e l s = c ( 弱 ， 中，强))
　　运行上述命令后，数据框cs2014month中增加了名为o3level的新列，其标注了每一行数据中臭氧的浓度水平。
　　2.2.2 分组变换。plyr包提供了一整套工
　　具集来处理列表(list)、数组(array)和数据框(data.frame)，它可以将复杂的数据分割成几个部分，分别对各个部分进行统计。对于数据框的操作，用的是plyr包中的ddply()函数，如：我们若想得到各个站点的各污染因子日均浓度值，则可调用如下命令：
　　c s d a t e m e a n - d d p l y ( c s 2 0 1 4 ， . ( d a t e ， s i t e ) ，
　　summarise，so2=round(mean(so2，na.rm=TRUE)，
　　3 )，n o 2 = r o u n d(me a n ( n o 2 ， n a . r m = TRU E)，
　　3)，co=round(mean(co，na.rm=TRUE)，2)，
　　pm10=round(mean(pm10，na.rm=TRUE)，3)，
　　pm2.5=round(mean(pm2.5，na.rm=TRUE)，3)，
　　o3=round(mean(sort(o3，decreasing=T)，
　　na.rm=TRUE)，3))
　　其中cs2014指原数据框，.(date，site)括号中指的是分组变量(日期、站点)，summarise命令意为从原数据中总结出一个新数据，round()为小数位数修约函数，mean()为平均函数，sort()为数据排序函数，命令中臭氧数据取其每日最大8小时平均值，其他污染物均为24小时平均值。该命令运行后就得到了一个新的数据框csdatemean。
　　类似的还可以用d d p l y ( ) 函数针对月份以及其他分类变量进行分组统计。如2.2.1节中提到的cs2014month，就可以调用如下命令得到：
　　cs2014month-ddply(cs2014，.(clock，month)，
　　summarise，so2=round(mean(so2，na.rm=TRUE)，
　　3 )，n o 2 = r o u n d(me a n ( n o 2 ， n a . r m = TRU E)，
　　3)，co=round(mean(co，na.rm=TRUE)，2)，
　　pm10=round(mean(pm10，na.rm=TRUE)，3)，
　　pm2.5=round(mean(pm2.5，na.rm=TRUE)，3)，
　　o3=round(mean(o3，na.rm=TRUE)，3))
　　除此之外，plyr包还提供了一些灵活的函数如transform等用于组间数值的统计，详见plyr包的自带说明。
　　2.2.3 数据变形。reshape2包提供的melt()和dcast()函数实现的是类似于Excel中数据透视表的功能，可从繁杂的数据中摘取出想要的信息。其中melt()函数可以将多列数据筛选后融合在一起，即从宽变长，下面举例说明：
　　csclockmean-ddply(cs2014，.(clock，season)，
　　summarise，so2=round(mean(so2，na.rm=TRUE)，
　　3)，no2=round(mean(no2，na.rm=TRUE)，
　　3)，co=round(mean(co，na.rm=TRUE)，2)，
　　pm10=round(mean(pm10，na.rm=TRUE)，3)，
　　pm2.5=round(mean(pm2.5，na.rm=TRUE)，3)，
　　o3=round(mean(o3，na.rm=TRUE)，3))
　　csclockmelt-melt(csclockmean，id.var=c(clock，season)，measure.var=c(so2，no2，pm2.5))
　　执行以上命令后，先得到一个名为csclockmean的数据框，再通过m e l t ( ) 函数得到一个名为csclockmelt的新数据框，用head()函数分别读取两个数据框的前6行：
　　melt()函数中id.var意为标识变量，即在原数据框中需保留的列向量，measure.var意为度量变量，即在原数据框中需要融合的列向量。melt()函数生成的新数据框中保留了id.var的同时，生成了两列：一列名为variable，包含的是measure.var;一列名为value，即原数据框中measure.var的值。由于ggplot2的独特语法，使用melt()函数得到的新数据框可使ggplot2在同一张图中绘制多个参数。
　　3 ggplot2数据可视化应用
　　一般认为R语言有4套图形系统，即graphics、grid、lattice、ggplot2。其中graphics、grid、lattice被认为是基础图形系统，它们的绘图命令精炼、实用、运行速度块，制出的图形美观、实用，但是基础图形系统的制图方式极为生硬，函数各参数繁复，令人往往在绘图上耗费大量的时间，而不是在数据分析上。
　　ggplot2是2005年后才新出现的图形系统，它提供了统一的接口及一些选项，替代了基础图形系统那套繁杂的修补体系，令使用者能更多地关注于数据本身。当然ggplot2中也有大量底层命令可对图形作精细的修改。
　　3.1 数据的简单分布
　　箱线图是一种展示数据分布的方法， 箱内为25% 75%的值，箱中的横线代表中位值，箱子上下竖线表示上下的相邻值，超过相邻值的点称为外部点。使用ggplot2绘制常山县两个站点二氧化硫全年日均值分布的箱线图命令如下所示。
　　ggplot(csdatemean，aes(site，so2))+geom_boxplot()+xlab(空气自动监测站点)+ylab(expression(paste(二氧化硫日均浓度(mg/m^3))))
　　ggplot2有其特殊的语法，其特点是+号，任何图层都可以通过+号来实现。geom_boxplot()是ggplot2中绘制箱线图的函数，类似的还有geom_point()(绘制散点图)、geom_line()(绘制折线图)、geom_histogram()(绘制条形图)等。xlab()添加的是x轴上的文字说明，ylab()中的expression()函数是为了在图中实现mg/m3的表达式。
　　3.2 趋势线的添加
　　在数据可视化工作中，参数与参数之间往往有着或多或少的关系，添加趋势线是探索这种关系的一种桥梁。R语言在数据分析方面最大的优势就是其强大的数据拟合功能，它提供了线性、非线性、神经网络、支持向量机等大量数据拟合工具。这些工具均可以与ggplot2配合使用绘制出趋势图形，以检验其拟合成果。下例使用ggplot2中的stat_smooth()函数自动添加趋势线来探索PM2.5与PM10数据间的关系。
　　ggplot(csdatemean，aes(pm2.5，pm10))+geom_point(colour=grey60)+stat_smooth(method=lm，se=F，colour=black)+xlab(expression(paste(PM2.5日均数据(mg/m^3))))+ylab(expression(paste(PM10日均数据(mg/m^3))))
　　在geom_point()和stat_smooth()两个函数中添加颜色参数，以使得散点图显得颜色淡一点，而让趋势线更明显。stat_smooth()函数中method=lm指按直线方式添加趋势线。从图5中明显可见PM2.5与PM10呈现了一种极强的线性关系。
　　3.3 分组变量的体现
　　在数据可视化的过程中，数据分组可以帮助人们揭示更多的细节。ggplot2可将在数据处理阶段创建的分组变量以颜色、形状、大小等形式映射于图上。
　　如2.2.2节中得到了各污染物按月份的时刻浓度平均数据cs2014month，2.2.1节中用cut()函数将该数据框中臭氧浓度分为弱、中、强三个水平，用ggplot2绘制出二氧化氮的时刻浓度图，并用图中二氧化氮浓度点的大小来代表臭氧的浓度水平，命令如下：
　　ggplot(cs2014month，aes(clock，no2))+geom_point(aes(size=o3level))+xlab(时刻)+ylab(expression(paste(二氧化氮浓度(mg/m^3))))+scale_size_discrete(name=臭氧浓度水平)
　　此例geom_point()函数中的size参数将臭氧浓度水平映射到二氧化氮浓度点的大小上，如此从图中可以明显可见臭氧浓度水平强时，往往是二氧化氮浓度较低时。类似的还可以用colour、shape等将分组变量的各组别映射为颜色、形状等，使绘制出的二维图呈现出更多的内容。
　　3.4 图形分面
　　数据可视化中最实用的技术之一是将分组的数据并列呈现，这样可以使得组间数据的比较变得容易许多，ggplot2中将之称为分面(facet)。ggplot2中可以用facet_grid()和facet_wrap()两个函数都可以绘制出分面图。
　　2.2.3中利用melt()函数得到了csclockmelt数据框，下例利用该数据，绘制四季中二氧化硫、二氧化氮、细颗粒物这三大污染物在24小时中的分布情况图，命令如下。
　　csclockmelt$season-factor(csclockmelt$season，levels=c(春，夏，秋，冬))
　　ggplot(csclockmelt，aes(clock，value))+geom_point(aes(shape=variable))+geom_line(aes(colour=variable))+facet_wrap( season)+xlab
　　(时刻)+ylab(expression(paste(浓度值(mg/m^3))))+scale_colour_hue(name=污染物，labels=c(二氧化硫，二氧化氮，细颗粒物))+scale_shape(name=污染物，labels=c
　　(二氧化硫，二氧化氮，细颗粒物)，solid=F)
　　命令中f a c t o r () 函数为季节排序， s h a p e 、colour将污染物种类映射至散点形状、曲线颜色，facet_wrap()函数自动将分组变量season进行分面处理。如此在相同的坐标系中得到四张并列的图，可令人更直观地比较数据。
　　4 结语
　　R语言最初时只是一个统计分析软件包，它擅长于绘图、分析数据以及利用数据来拟合统计模型，在不擅长的领域如复杂数据存储、大数据处理，它也提供了接口;ggplot2以其自成一派的语法闻名，它不需要考虑某图形需要填充什么颜色，使用什么形状的线段连接散点，不必去纠结基本图形系统中那些繁杂的参数，将更多的时间用于数据分析本身。本文中某些例子由于考虑到各图形的美观及信息的正式性，在绘制时应用了不少ggplot2中的底层命令，显得各例中命令较为复杂。实际应用中均可大大简化，即使如此，与基础图形系统相比，ggplot2的优势也尽显无疑。
　　当然，工具永远只能是工具，无论什么工具都无法替代有效数据的积累和大脑思考问题的角度。各地将气象数据纳入环境空气质量数据一起分析时，思考时角度可以更加多元化;各地的一些污染控制工作内容也可以纳入数据框中作分析，以评价其有效性;省、市以上的监测部门更是可以利用更多的环境空气自动监测数据作一些地域性的探索。