前几天发现了一个很有趣的包——openair，可以将年度时间序列刻画成周年日历热图，感觉这种形式非常适合用于呈现年度空气质量可视化，所以抓空爬了一些大连市2016年年度空气质量数据拿来玩玩，目标网站网页结构比较简单，爬取过程很轻松，界面部分很规律，感觉这个代码可以作为模板用，感兴趣的小伙伴儿可以试着玩一玩！

library(RCurl)

library(XML)

library(dplyr)

library(ggplot2)

library(stringr)

library(rvest)

library(lubridate)

library("DT")

library(openair)

library(ggplot2)

数据爬取过程：

构造月度url地址（网站是按照月度数据存储的，需要按月爬取）

urlbase<-"https://www.aqistudy.cn/historydata/"

url<-"https://www.aqistudy.cn/historydata/monthdata.php?city=大连"

rd <- getURL(url,.encoding="UTF-8")

rdhtml <- htmlParse(rd,encoding="UTF-8")

otherpage<-getNodeSet(rdhtml,"//a")

allurl<-laply(otherpage,xmlGetAttr,name='href')%>%grep("2016",.,value=T)%>%sub("麓贸脕卢","大连",.)

#以上还是编码出了问题，不知道那个乱码是什么鬼！只能强行替换了！

allurl<-paste0(urlbase,allurl)

以上过程也可先通过观察大连市的月度空气质量url地址规律，然后通过paste函数直接生成。

allurl<-paste0("https://www.aqistudy.cn/historydata/daydata.php?city=大连&month=",201601:201612)

这是简单粗暴的方式，但是不保证任何网址都可以使用

先写完一个看下具体情况

tbls<-read_html(allurl[1],encoding="utf-8")%>%html_table(.,header=TRUE,trim=TRUE);tbls<-tbls[[1]]

编写单次爬取函数，使用for循环遍历网址进行数据获取（原谅我又用了for循环）

mytable<-data.frame()

for (i in allurl){

Sys.sleep(sample(1:5,1))

fun<-function(m){

table<-read_html(m,encoding="utf-8")%>%html_table(.,header=TRUE,trim=TRUE)

table<-table[[1]]

}

mytable<-rbind(mytable,fun(i))

}

使用动态表格查看数据

datatable(mytable)

https://ask.hellobi.com/uploads/article/20170329/9e8ab75a0cee5deaba897cb3bffc1fca.png

#备份一份数据，以防原数据损坏

mytableb<-mytable

调整时间变量

mytable$日期<-as.Date(mytable$日期);names(mytable)[1]<-"date"

AQI指数年度分布日力图

calendarPlot(mytable,pollutant="AQI",year=2016)

https://ask.hellobi.com/uploads/article/20170329/4297db9350ec3dcfaf55c495e4a64fd0.png

PM2.5指数年度分布日力图

calendarPlot(mytable,pollutant="PM2.5",year=2016)

https://ask.hellobi.com/uploads/article/20170329/2ba8975667c50fa924aa2637a4462a3d.png

--------------------------------------------------------------------------------

接下来使用ggplot函数制作同样的日力图

dat <- mytable

这次使用lubridate包来处理时间日期变量（超级好用）

dat$month<-as.numeric(as.POSIXlt(dat$date)$mon+1)

dat$monthf<-factor(dat$month,levels=as.character(1:12),labels=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"),ordered=TRUE)

dat$weekday<-as.POSIXlt(dat$date)$wday

dat$weekdayf<-factor(dat$weekday,levels=rev(0:6),labels=rev(c("Sun","Mon","Tue","Wed","Thu","Fri","Sat")),ordered=TRUE)

dat$week <- as.numeric(format(dat$date,"%W"))

dat<-ddply(dat,.(monthf),transform,monthweek=1+week-min(week))

AQI指数为污染级别以上的天数分布

windowsFonts(myFont = windowsFont("微软雅黑"))

ggplot(dat,aes(monthweek,weekdayf,fill=AQI))+

  geom_tile(colour='white') +

  facet_wrap(~monthf ,nrow=3) +

  scale_fill_gradient(space="Lab",limits=c(100, max(dat$AQI)),low="yellow", high="red") +

  labs(title="大连市2016年空气日历热图",subtitle="AQI指数为污染级别以上的天数分布（AQI>=100）",x="Week of Month",y="")+

  theme(title=element_text(family="myFont"))

https://ask.hellobi.com/uploads/article/20170329/e694a737399ecfa5b00f05672cb48f41.png

PM2.5指数为污染级别以上的天数分布

ggplot(dat,aes(monthweek,weekdayf,fill=PM2.5))+

  geom_tile(colour='white') +

  facet_wrap(~monthf ,nrow=3) +

  scale_fill_gradient(space="Lab",limits=c(75, max(dat$PM2.5)),low="yellow", high="red") +

  labs(title="大连市2016年气温日历热图",subtitle="PM2.5指数为污染级别以上的天数分布（PM2.5>=75）",x="Week of Month",y="")+

  theme(title=element_text(family="myFont"))

https://ask.hellobi.com/uploads/article/20170329/6d4844807e33832c4ad0af58a706b584.png

大体来看，大连整个2016年度污染天气相对来讲，还算是挺良心的，跟帝都比起来要好很多。AQI和PM2.5在污染级别以上的天数不超过两个月。

-------------------------------------------------------------------------------------------------------

接下来呢，我们做一些详细的统计工作，具体就是从时间细分维度来查看季度、月度、周度等平均AQI、PM2.5指数分布情况

data3<-mytable[,c(1,2,4,5)]

data3<-transform(data3,Quarter=quarter(date),Month=month(date),Week=week(date))

data3$质量等级<-factor(data3$质量等级,levels=c("重度污染","中度污染","轻度污染","良","优"),labels=c("重度污染","中度污染","轻度污染","良","优"),order=T)

首选查看五个污染级别在2016年度出现的频率：

countd<-count(data3$质量等级)%>%arrange(.,-freq)

ggplot(countd,aes(reorder(x,freq),freq))+

geom_bar(fill="#0C8DC4",stat="identity")+

coord_flip()+

labs(title="大连市2016年度空气质量分布",subtitle="污染级别频率分布图",caption="https://www.aqistudy.cn/")+

geom_text(aes(label=freq),hjust=1,colour="white",size=8)+

theme_bw()+

theme(

      text=element_text(family="myFont"),

      panel.border=element_blank(),

      panel.grid.major=element_line(linetype="dashed"),

      panel.grid.minor=element_blank(),

      plot.caption=element_text(hjust=0,size=10),

      axis.title=element_blank()

      )

https://ask.hellobi.com/uploads/article/20170329/d7be764b6f1aa70955f415d3142a4918.png

基于季度空气质量平均水平分布图：

Quarter<-aggregate(AQI~Quarter,data=data3,FUN=mean)

ggplot(Quarter,aes(reorder(Quarter,-AQI),AQI))+

geom_bar(fill="#0C8DC4",stat="identity")+

labs(title="大连市2016年度空气质量分布",subtitle="AQI污染指数季度指标平均分布图",caption="https://www.aqistudy.cn/")+

geom_text(aes(label=round(AQI)),vjust=1.5,colour="white",size=8)+

theme_bw()+

theme(

      text=element_text(family="myFont"),

      panel.border=element_blank(),

      panel.grid.major=element_line(linetype="dashed"),

      panel.grid.minor=element_blank(),

      plot.caption=element_text(hjust=0,size=10),

      axis.title=element_blank()

      )

https://ask.hellobi.com/uploads/article/20170329/f0bbcb8fc866a86626dcc51d2135ae3c.png

基于月度空气质量平均水平分布图：

Month<-aggregate(AQI~Month,data=data3,FUN=mean)

ggplot(Month,aes(reorder(Month,-AQI),AQI))+

geom_bar(fill="#0C8DC4",stat="identity")+

labs(title="大连市2016年度空气质量分布",subtitle="AQI污染指数月度指标平均分布图",caption="https://www.aqistudy.cn/")+

geom_text(aes(label=round(AQI)),vjust=1.5,colour="white",size=6)+

theme_bw()+

theme(

      text=element_text(family="myFont"),

      panel.border=element_blank(),

      panel.grid.major=element_line(linetype="dashed"),

      panel.grid.minor=element_blank(),

      plot.caption=element_text(hjust=0,size=10),

      axis.title=element_blank()

      )

基于月度空气质量平均水平分布图：

https://ask.hellobi.com/uploads/article/20170329/f226d81f6f84e7e5b1bfa7bdb5848a6a.png

Week<-aggregate(AQI~Week,data=data3,FUN=mean)

ggplot(Week,aes(factor(Week,order=T),AQI,group=1))+

geom_line(col="#0C8DC4")+

labs(title="大连市2016年度空气质量分布",subtitle="AQI污染指数周度指标平均分布图",caption="https://www.aqistudy.cn/")+

geom_text(aes(label=ifelse(Week>100,Week,"")),vjust=1.5,colour="white",size=6)+

theme_bw()+

theme(

      text=element_text(family="myFont"),

      panel.border=element_blank(),

      panel.grid.major=element_line(linetype="dashed"),

      panel.grid.minor=element_blank(),

      plot.caption=element_text(hjust=0,size=10),

      axis.title=element_blank()

      )

https://ask.hellobi.com/uploads/article/20170329/46e89ef69ff8ed7f76a95b8e1a4c061e.png

从以上周度AQI平均指标上来看，大连市2016年度的周度平均AQI指数大部分周都在100以下，看到这个感觉生活在大连还是蛮幸福的，看着北上的小伙伴儿隔三差五的在朋友圈晒人间仙境也是一件很有趣的事哈哈！