R语言对某地天气和温度的分析及预测

历史天气数据来源：http://tianqi.2345.com/wea_history/54511.htm，这是北京的历史数据，采样城市北京、上海、苏州、长沙、广州、一共采集了2011-1-1到2015-4-2这四年三个月共1542（354+366+365+365+92）天的天气数据，其中2011-1-17到2011-1-25这十天的数据缺失，查了多个网站都发现这种情况，就没有把数据补齐了。另外，上海、苏州、广州这三个城市2012-1-15这天，长沙2015-2-10这天，五个城市2014-3-8这天的数据缺失，这里根据前后两天的温度取平均值，天气定为阴，风向后风级都根据前后天补齐。下面两个图是在两个不同的网站上查到的历史天气数据都是有缺失的。

将数据做完清洗整理后，存做csv格式，数据输入R中，并查看数据基本结构如下：

[plain] view plaincopy

1. beijing <-read.csv("BeiJing.csv",header=T,stringsAsFactors=FALSE)
2. changsha <-read.csv("Changsha.csv",header=T,stringsAsFactors=FALSE)
3. guangzhou <-read.csv("GuangZhou.csv",header=T,stringsAsFactors=FALSE)
4. shanghai <-read.csv("ShangHai.csv",header=T,stringsAsFactors=FALSE)
5. suzhou <-read.csv("SuZhou.csv",header=T,stringsAsFactors=FALSE)
7. head(suzhou,n=5)

每一列分别是：date(日期)、highestTemp(最高温度)、lowest(最低温度)、weather(天气)、wind(风向)、windForce(风级)。

天气篇

统计的基本天气类型：雪、雨、晴、阴，优先级顺序也是这样，就是说如果是雨雪天气，记为雪，晴转多云记为晴，另外多云和阴都记为阴，这里有个链接说明了气象中阴和多云的差别，主要是云量大小的差异，http://www.guokr.com/question/252793/，这里我把它统一算作阴，还有雾、霾、浮尘天气都记做阴。

[plain] view plaincopy

1. weatherStat <- function(x){
2. count <- numeric(0)
3. count[1:5] <- 0
4. for(i in 1:length(x[[4]])){
5. if(length(grep("雪",x[[4]][i]))>0){
6. count[1] <- count[1]+1
7. }else if(length(grep("雨",x[[4]][i]))>0){
8. count[2] <- count[2]+1
9. }else if(length(grep("晴",x[[4]][i]))>0){
10. count[3] <- count[3]+1
11. }else if((length(grep("阴",x[[4]][i]))>0)||(length(grep("多云",x[[4]][i]))>0)||(length(grep("雾",x[[4]][i]))>0)||(length(grep("霾",x[[4]][i]))>0)||(length(grep("浮尘",x[[4]][i]))>0)){
12. count[4] <- count[4]+1
13. }else{
14. count[5] <- count[5]+1
15. print(x[[4]][i])
16. }
17. }
18. count
19. }
21. statAll<-list(beijing=numeric(0),suzhou=numeric(0),shanghai=numeric(0),changsha=numeric(0),guangzhou=numeric(0))
22. statAll$suzhou<- weatherStat(suzhou)
23. statAll$beijing<- weatherStat(beijing)
24. statAll$shanghai<- weatherStat(shanghai)
25. statAll$changsha<- weatherStat(changsha)
26. statAll$guangzhou<- weatherStat(guangzhou)
28. statAll<- as.data.frame(statAll)
29. statAll<- statAll[-5]#第五行，是统计除这四种天气外还是否有其他情况，五个城市都为0，所以删除这一行，前面雾、霾、浮尘都是在这个类别下发现然后整理到阴下面的
30. rownames(statAll)<- c("雪","雨","晴","阴")
31. colnames(statAll)<- c("北京","苏州","上海","长沙","广州")
32. statAll
33. statAll<- as.matrix(statAll)
34. barplot(statAll,legend=TRUE,col=c("snow3","lavender","khaki1","lemonchiffon"))

这个简单的统计可以看出基本上越往南走，晴天越多，下雨天则是相比而言南方比北方更多，下雪天则是北方比南方多，不过这里北方就只取了北京一个城市。到这里我发现对苏州我似乎还有点误解，相比曾经呆过的长沙而言，这里的雨天还比较少，为什么我总有种这里一到放假就下雨的感觉呢？

下面只看苏州的天气情况，看这几年的一个统计：

[plain] view plaincopy

1. #2011年1:354,2012年355:720,2013年721:1085,2014年1086:1450,2015年1451:1542
2. statSuzhou<-list(one=numeric(0),two=numeric(0),three=numeric(0),four=numeric(0),five=numeric(0))
3. statSuzhou$one<- weatherStat(suzhou[1:354])
4. statSuzhou$two<- weatherStat(suzhou[355:720])
5. statSuzhou$three<- weatherStat(suzhou[721:1085])
6. statSuzhou$four<- weatherStat(suzhou[1086:1450])
7. statSuzhou$five<- weatherStat(suzhou[1451:1542])
8. statSuzhou<- as.data.frame(statSuzhou)
9. statSuzhou<- statSuzhou[-5]#第五行，是统计除这四种天气外还是否有其他情况，五个城市都为0，所以删除这一行
10. rownames(statSuzhou)<- c("雪","雨","晴","阴")
11. colnames(statSuzhou)<-c("2011","2012","2013","2014","2015")
12. statSuzhou
13. statSuzhou<- as.matrix(statSuzhou)
14. barplot(statSuzhou,legend=TRUE,col=c("snow3","lavender","khaki1","lemonchiffon"))

苏州天气，晴天的时间除了2013年都在30%以下，2013年天气都挺好的？后面可以看下气温在2013年有没有什么异常。

那么不同月下雨有什么差异呢？能看出集中的降雨时间吗？接下来根据月份做一个天气的统计。

？查了一下正常的梅雨季节应该是六月到七月，我最前面的梅雨季节的说法不恰当，这里就不改了。

[plain] view plaincopy

1. statMonthly<- data.frame()
2. statMonthly[1:5,1]
3. count<- 1
4. for(i in2011:2015){
5. if(i != 2015){
6. for (j in 1:12) {
7. temp <-nrow(suzhou[grep(paste(i,"-",j,"-",sep=""),suzhou$date)])
8. statMonthly[1:5,((i-2011)*12+j)]<- weatherStat(suzhou[count:(count+temp-1)])
9. print(weatherStat(suzhou[count:(count+temp-1)]))
10. print(c(((i-2011)*12+j)))
11. count <- count + temp
12. }
13. }else{
14. for (j in 1:3) {
15. temp <-nrow(suzhou[grep(paste(i,"-",j,"-",sep=""),suzhou$date)])
16. statMonthly[1:5,((i-2011)*12+j)]<- weatherStat(suzhou[count:(count+temp-1)])
17. print(weatherStat(suzhou[count:(count+temp-1)]))
18. print(c(((i-2011)*12+j)))
19. count <- count + temp
20. }
21. }
22. }
24. statMonthly<- statMonthly[-5]
25. rownames(statMonthly)<- c("雪","雨","晴","阴")
26. colnames(statMonthly)<-c(paste(rep(2011,12),1:12,sep="-"),paste(rep(2012,12),1:12,sep="-"),paste(rep(2013,12),1:12,sep="-"),paste(rep(2014,12),1:12,sep="-"),paste(rep(2015,3),1:3,sep="-"))
27. statMonthly<- as.matrix(statMonthly)
28. barplot(statMonthly,legend=TRUE,col=c("snow3","lavender","khaki1","lemonchiffon"))
29. barplot(statMonthly,col=c("snow3","lavender","khaki1","lemonchiffon"),cex.names=.6)
30. lines(statMonthly[2],type="l",col="red")
32. statMonthly

上面是按月统计的天气结果图，可以看出，雨雪天气在4-8月会比较集中，上图中用红线标出来的部分是对应的每年的4-8月份。并且，2013年，雨雪天相比而言比较少。不过横轴太长，看起来很费力。

下面不考虑2015年的数据，看2011-2015年，每个月的天气统计情况。

[plain] view plaincopy

1. statMonthly2<- data.frame()
2. temp<- 1:51
3. for(i in1:12){
4. if(i==12){
5. for(j in 1:4){
6. statMonthly2[j,12]<- sum(statMonthly[j,temp%%12==0&temp<49])
7. }
8. }else{
9. for(jin 1:4){
10. statMonthly2[j,i]<- sum(statMonthly[j,temp%%12==i&temp<49])
11. }
12. }
13. }
14. rownames(statMonthly2)<- c("雪","雨","晴","阴")
15. colnames(statMonthly2)<- 1:12
16. statMonthly2<- as.matrix(statMonthly2)
17. barplot(statMonthly2,col=c("snow3","lavender","khaki1","lemonchiffon"),cex.names=.8)

这里就能清楚看出来，降雨量比较多的月份从是6-8月，正是梅雨季，但是其他月份的降雨量相差的不是很大。比较明显的是北京，如下图：

这里主要是对天气的一个统计分析，能看出苏州是一个很典型的江南城市，雨雪天气比较多，并且全年都有，6-8月稍多。分析的这四年中，2013年比较特别，雨天比较少天晴的时间比较多。而我比较关注的四月，从统计结果看起来，阴晴雨的时间看起来很平均。又刮风下雨时，请记住这里晴天的几率只有不到百分之三十，所以阴雨天，正常！

温度篇

前面已经讲了苏州的天气特点，还是用相同的数据，做接下来的苏州气温特点的分析预测，是的预测在这里！

首先看下2011年到2015年苏州整体的温度表现是什么样的。

[plain] view plaincopy

1. plot(suzhou$highestTemp,type="l",col="red",main="苏州2011-2015年气温图",xlab="时间轴",ylab="温度℃")
2. lines(suzhou$lowestTemp,type="l",col="blue")
3. legend("topright",c("最高气温","最低气温"),col=c("red","blue"),lty=1)

红色是最高气温，蓝色是最低气温，年度季节性的特征很明显。每年都是先升再降，7、8月份是温度最高的时间，1、2月是温度最低时间。

因为时间太长，横轴没有具体的对应点。

同样，按月取平均值，再来看整体的表现。

[plain] view plaincopy

1. avgTemper <-numeric(0)#月平均气温
2. diffTemper <-numeric(0)#月最高温差
3. length(avgTemper)<- 48
4. length(diffTemper)<- 48
6. for(i in2011:2014){
7. for(j in 1:12){
8. print((i-2011)*12+j)
9. avgTemper[(i-2011)*12+j]<-mean(c(suzhou$highestTemp[grep(paste(i,"-",j,"-",sep=""),suzhou$date)],suzhou$lowestTemp[grep(paste(i,"-",j,"-",sep=""),suzhou$date)]),na.rm=TRUE)
10. diffTemper[(i-2011)*12+j]<- max(suzhou$highestTemp[grep(paste(i,"-",j,"-",sep=""),suzhou$date)]-suzhou$lowestTemp[grep(paste(i,"-",j,"-",sep=""),suzhou$date)],na.rm=TRUE)
11. }
12. }
14. avgTemperTS<- ts(avgTemper,frequency=12,start=c(2011,1))
15. plot.ts(avgTemperTS,main="苏州2011-2014年月平均气温图",xlab="时间",ylab="月平均温度℃")

四年月平均气温有一个很明显的周期性规律，明显能看出来2013年7、8月份气温高于其他三年。

换一种方式，将每年的数据放在同一个时间轴上，可以更明显看出来每年的温度走势都非常接近，并且对比图中，2013年夏季比其他几年都热，查了一下，最高气温能达到41摄氏度，是2013-8-8，当天最高气温41度，最低气温也有31度。

[plain] view plaincopy

1. plot.ts(avgTemperTS[25:36],col="red",main="苏州2011-2014年月平均气温图",xlab="月份",ylab="月平均温度℃")
2. lines(avgTemperTS[1:12],col="black")
3. lines(avgTemperTS[13:24],col="blue")
4. lines(avgTemperTS[37:48],col="green")
5. legend("topright",c("2011","2012","2013","2014"),col=c("black","blue","red","green"),lty=1)

[plain] view plaincopy

1. which.max(suzhou[,2])
2. suzhou[940]

那么苏州的温差又有什么特点呢，什么时候是一年中温差最大的时候？同样看四年的整体表现和对比图。

[plain] view plaincopy

1. diffTemperTS <-ts(diffTemper,frequency=12,start=c(2011,1))
2. plot.ts(diffTemperTS,main="苏州2011-2014年月最大温差图",xlab="时间",ylab="月最大温差℃")
3. plot.ts(diffTemperTS[25:36],col="red",ylim=c(8,23),main="苏州2011-2014年月最大温差图",xlab="月份",ylab="月最大温差℃")
4. lines(diffTemperTS[1:12],col="black")
5. lines(diffTemperTS[13:24],col="blue")
6. lines(diffTemperTS[37:48],col="green")
7. legend("topright",c("2011","2012","2013","2014"),col=c("black","blue","red","green"),lty=1)

上面两个图中，都有一个好大的波峰，没错又是任性的2013年。整体上来看温差比较大的时间是3、4月，正是现在，换季的时节天气很是自在随性。那么那个孤高的波峰是哪一天呢？我们可以看看。

[plain] view plaincopy

1. which.max(suzhou$highestTemp[grep(paste(2013,"-",3,"-",sep=""),suzhou$date)]-suzhou$lowestTemp[grep(paste(2013,"-",3,"-",sep=""),suzhou$date)])
2. suzhou[grep("2013-3-9",suzhou$date)]

最后，我想用时间序列来做温度走势的预测。因为R中周期性的时间序列数据需要每期数据项相同，但是这里2011年-2014年每年的数据项都不同（天气篇中开篇讲到了分别是354、366、365、365），所以我把这个时间序列的周期定为360，按顺序取数据的话，2011年因为1月份少了10天，所以2011、2012年会有一些时间和对应的温度错位。考虑到温度是连续型变量，小范围的变动应该不会对整体的时间序列的预测产生很大的影响。

最高温度和最低温度分别定为一个时间序列，并将其在同一个图上显示出来。

[plain] view plaincopy

1. highestTS <-ts(suzhou$highestTemp[1:(360*4)],frequency=360,start=c(2011,1,1))
2. plot.ts(highestTS,col="red")
3. lowestTS <-ts(suzhou$lowestTemp[1:(360*4)],frequency=360,start=c(2011,1,1))
4. lines(lowestTS,col="blue")

这是周期性特别明显的数据，所以考虑使用HoltWinters指数平滑方法来做时间序列的分析预测，用这个方法分别对最高温时间序列数据和最低温时间序列数据分别做平滑得到平滑模型然后来做气温的预测。

[plain] view plaincopy

1. highestForecasts<- HoltWinters(highestTS)
2. highestForecasts
3. lowestForecasts<- HoltWinters(lowestTS)
4. lowestForecasts

上面模型的平滑参数alpha、beta、gamma，分别是平滑指数、趋势指数、季节指数，参考http://blog.csdn.net/cl1143015961/article/details/42267613，值alpha越大说明对远期也就是以前年份的数据参考权重越大，另外gamma越大说明这个时间序列的季节性特点会越明显。

接下来对这两个时间序列做预测。要用到forecast包，定义预测150天的气温。

[plain] view plaincopy

1. library("forecast")
2. highestForecast2<- forecast.HoltWinters(highestForecasts,h=150)
3. lowestForecast2<- forecast.HoltWinters(lowestForecasts,h=150)
4. plot(highestForecast2$mean,ylim=c(-4,32),col="red",main="苏州2015年最高最低气温预测时序图",xlab="时间",ylab="温度℃")
5. lines(lowestForecast2$mean,col="blue")
6. legend("topright",c("最高气温","最低气温"),col=c("red","blue"),lty=1)

这个是预测的最高最低气温的图，但是到底预测结果可信度高不高呢？

[plain] view plaincopy

1. Box.test(highestForecast2$residuals,lag=20,type="Ljung-Box")
2. Box.test(lowestForecast2$residuals,lag=20,type="Ljung-Box")

参考文章里面讲到了几种检验方式，这里我们用Ljung-Box检验方式，可以看到两个时间序列的检验结果p值都非常小，所以模型的预测是有效的。

突然想起来，2015年不是已经过了三个月了嘛，我们也有前面三个月的数据，那么接下来就把预测结果和实际结果放在一个图中对比看预测效果。本来想把四条线放在一起查看，但是比较混乱所以还是分开看。

[plain] view plaincopy

1. highestFactTS <- ts(suzhou$highestTemp[1451:1542],frequency=360,start=c(2015,1))
2. lowestFactTS <- ts(suzhou$lowestTemp[1451:1542],frequency=360,start=c(2015,1))
3. plot(highestForecast2$mean,col="red",main="苏州2015年最高气温预测观察值时序图",xlab="时间",ylab="温度℃")
4. lines(highestFactTS,col="orange")
5. legend("topleft",c("实际值","预测值"),col=c("orange","red"),lty=1)

[plain] view plaincopy

1. plot(lowestForecast2$mean,col="blue",main="苏州2015年最低气温预测观察值时序图",xlab="时间",ylab="温度℃")
2. lines(lowestFactTS,col="green")
3. legend("topleft",c("实际值","预测值"),col=c("green","blue"),lty=1)

没有重合，怎么可以不重合！开玩笑，重合是不太可能的，从上面两个最高气温和最低气温的预测值与观察值的图都能看出来预测值比实际值都有些偏低。但是预测结果与实际结果相比相差不大，并且整个变化趋势还是接近的。给出这个模型预测的接下来苏州地区50天的温度预测，当然，这个结果只能做一个参考，你是愿意看天气预报呢，还是看它呢？

[plain] view plaincopy

1. predictTemp <-list(date="",highestTemp="",lowestTemp="")
2. predictTemp$date <-c(paste(rep(2015,50),c(rep(4,20),rep(5,30)),c(11:30,1:20),sep="-"))
3. predictTemp$highestTemp <-round(highestForecast2$mean[101:150])
4. predictTemp$lowestTemp <-round(lowestForecast2$mean[101:150])
5. predictTemp <-as.data.frame(predictTemp)

有没有好奇天气预报到底是怎么做出来的？想到这里，顺手百度了一下，大家可以看看，http://www.guokr.com/question/273034/，主要是云层图和历史数据，还要天气预报人员的经验！