ggplot2

wind rose(風配図) 3

ggplot2,XML パッケージ

(過去の記事)
wind rose(風配図)
wind rose(風配図) 2

今回は、ggplot2 パッケージを使ってwind rose(風配図)を作成します。

(データ)
気象庁 過去の気象データ検索

鳥取県：鳥取　2015年1月1日（１時間ごとの値）

福井県：小浜　2015年1月1日（１時間ごとの値）

鳥取県：鳥取

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
library(ggplot2)
library(XML)
#library(RCurl)
#
#空のデータフレームを作成
weather<-read.table(text = "",col.names = c("date","time","pressure","precipitation","temperature","wind_velocity","wind_direction"))
#2015年1月の風配図
year<-2015;month<-1
for (i in 1:31){
#データを入手するURLを指定
html.name <-paste("http://www.data.jma.go.jp/obd/stats/etrn/view/hourly_s1.php?prec_no=69&block_no=47746&year=",year,"&month=",month,"&day=",i ,sep="",collapse="")
t1<-readHTMLTable(html.name)
#時間、現地気圧、降水量、気温、風速、風向データ
t2<-t1$tablefix1[,c(1,2,4,5,9,10)]
t2<-t2[-1,]
date<-paste(year,"/",month,"/",i,sep="",collapse="")
d1<- data.frame(matrix(rep(date,24),ncol=1))
t2<-cbind(d1,t2)
names(t2)<-c("date","time","pressure","precipitation","temperature","wind_velocity","wind_direction")
#ポイント：as.numeric(as.vector())　と記述
#as.numeric()ではうまく型変更できない。
t2[,2]<-as.numeric(as.vector(t2[,2]))
t2[,3]<-as.numeric(as.vector(t2[,3]))
t2[,4]<-as.numeric(as.vector(t2[,4]))
t2[,5]<-as.numeric(as.vector(t2[,5]))
t2[,6]<-as.numeric(as.vector(t2[,6]))
weather<-rbind(weather,t2)
}
#数カ月分のデータをダウンロードしあとで月ごとに作図する場合
#zoo クラス　xts クラスにするとよい
#library(xts)
#t<-as.xts(read.zoo(weather))
#wind<-data.frame(month=substring(index(t),1,7),direction=as.vector(chartr("東西南北", "EWSN",t$wind_direction)),velocity=as.numeric(t$wind_velocity))
#
wind<-data.frame(direction=as.vector(chartr("東西南北", "EWSN",weather$wind_direction)),velocity=as.numeric(as.vector(weather$wind_velocity)))
#directionにはいっているのは方向とは限らない。"静穏"や"///"や"x"
wind<- subset(wind,direction=="N"|direction=="NNE"|direction=="NE"|direction=="ENE"|direction=="E"|direction=="ESE"|direction=="SE"|direction=="SSE"|direction=="S"|direction=="SSW"|direction=="SW"|direction=="WSW"|direction=="W"|direction=="WNW"|direction=="NW"|direction=="NNW")
#
#save(wind,file="wind_tottori201501.RData")
#load("wind_tottori201501.RData")
#
ord<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
wind$direction<- factor(wind$direction,levels=ord)
#
wind$cut<-cut(wind$velocity,breaks=c(0,3,6,10,20,100),labels=c("[0,3]","(3,6]","(6,10]","(10,20]" ,"(20,  ") ,right =TRUE, include.lowest = TRUE) # 
#
table(wind$cut,wind$direction)
#
#円が16個に分割されない（その方向のデータがない場合）
#ggplot(wind, aes(x =factor(direction),fill = cut)) + geom_bar(width =0.8) + coord_polar(theta = "x",start = -0.5) + 
#     scale_fill_brewer(palette = "Accent")
#
#ダミーのデータ（すべての方向、風速は欠損値）を結合する。
dummy<-data.frame(direction=levels(wind$direction),velocity=rep(NA,16))
newwind<-rbind(wind[,c("direction","velocity")],dummy)
ord<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
newwind$direction<- factor(newwind$direction,levels=ord)
newwind$cut<-cut(newwind$velocity,breaks=c(0,3,6,10,20,100),labels=c("[0,3]","(3,6]","(6,10]","(10,20]" ,"(20,  ") ,right =TRUE, include.lowest = TRUE) 
#png("tottoroWR201501.png")
ggplot(newwind, aes(x =factor(direction),fill = cut)) + geom_bar(width =0.8) + coord_polar(theta = "x",start = -0.2) + 
     scale_fill_brewer(palette = "Set2") +
     ylim(0,170) +
     theme(legend.position="bottom")+
     labs(fill = "Velocity",x="Direction",y="",title=paste0("風配図(鳥取:",year,"年",month,"月)"))
#dev.off()

福井県：小浜

グラフ作成部分のコード

2015年
1
2
3
4
5
6
7
8
9
10
11
12
library(ggplot2)
#保存したデータを読み込んだ場合は以下が必要
#days<-seq(as.Date("2015-01-01"),as.Date("2015-12-31"), by="days")
#year<-"2015"
#
#1年分のデータならdummyデータは不要のはず
#png("obama2015.png")
ggplot(wind, aes(x =factor(direction),fill = cut)) + geom_bar(width =0.8) + coord_polar(theta = "x",start = -0.2) + 
     scale_fill_brewer(palette = "Set2") +
     theme(legend.position="bottom")+
     labs(fill = "Velocity",x="Direction",y="",title=paste0("風配図(小浜:",year,"年)"))
#dev.off()
2015年 月ごと
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#ダミーのデータ（すべての方向、風速は欠損値）を結合する。
dummy<-data.frame(month=rep(unique(substring(days,1,7))[1],16),direction=levels(wind$direction),velocity=rep(NA,16))
for (i in 2:12){
	dummy<-rbind(dummy,data.frame(month=rep(unique(substring(days,1,7))[i],16),direction=levels(wind$direction),velocity=rep(NA,16)))
	}
#
newwind<-rbind(wind[,c("month","direction","velocity")],dummy)
#
ord<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
newwind$direction<- factor(newwind$direction,levels=ord)
newwind$cut<-cut(newwind$velocity,breaks=c(0,3,6,10,20,100),labels=c("[0,3]","(3,6]","(6,10]","(10,20]" ,"(20, ") ,right =TRUE, include.lowest = TRUE) 
#
#png("obama2015_2.png")
ggplot(newwind,aes(x =factor(direction),fill=cut)) + geom_bar(width = 0.8) + coord_polar(theta = "x",start = -0.2) +
         scale_fill_brewer(palette = "Set2") +
         theme(legend.position="bottom")+
         ylim(0,170) +
         labs(fill = "Velocity",x="Direction",y="",title=paste0("風配図(小浜:",year,"年)"))+
	 facet_wrap(~month) 
#dev.off()
アニメーション

gganimate パッケージではうまくいかなかったので

1
2
3
4
5
6
7
8
9
10
for (i in 1:12){
 p<-ggplot(subset(newwind,month==unique(substring(days,1,7))[i]), aes(x =factor(direction),fill=factor(cut))) + geom_bar(width = 0.8) + coord_polar(theta = "x",start = -0.2) +
         scale_fill_brewer(palette = "Set2") +
         theme(legend.position="bottom")+
         labs(fill = "Velocity",x="Direction",y="",title=paste0("風配図(小浜:",year,"年",i,"月)")) +
         ylim(0,170)
  ggsave(file = paste0("/home/user/Rwork/anime/obamaWR2015_",unique(substring(days,6,7))[i],".png"), plot = p, dpi = 100, width = 4.8, height = 4.8)
		}
#
system("convert -delay 200 -loop 0 /home/user/Rwork/anime/obama*.png /home/user/Rwork/anime/obama2015movie.gif")

複数月に渡るデータを取り込む
days<-seq(as.Date(“2015-01-01”),as.Date(“2015-12-31”), by=”days”)
for 文
year<-substring(days[i],1,4)
・・・

とした。