グラフまとめ

wind rose(風配図)

climatol パッケージの rosavent 関数

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

(データ使用に関する注意事項)
気象庁ホームページについて

(パッケージもしくは関数の入手は)
CLIMATOL
「Warning! The package at CRAN repositories is unusable with recent versions of R!!! 」という警告あり

  • R2.13.0ではCRANを通してインストールできる
  • 新しいバージョンのRの場合、パッケージをダウンロードしてインストールするか。関数をダウンロードしてsource関数で取り込む。

この記事では以下のようなコードでHPから直接取り込みます。

source(“http://www.meteobal.com/climatol/rosavent.R“)

climatol パッケージのサンプルデータ windfr の書式を見てみると、
data(windfr)
library(xtable)
print(xtable(windfr),”html”)








N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
0-3 59 48 75 90 71 15 10 11 14 20 22 22 24 15 19 33
3-6 3 6 29 42 11 3 4 3 9 50 67 28 14 13 15 5
6-9 1 3 16 17 2 0 0 0 2 16 33 17 6 5 9 2
&gt 9 0 1 2 3 0 0 0 0 0 1 4 3 1 1 2 0

やること

  1. 気象庁のHPから鳥取市の気象データを入手し、 XML パッケージを利用してデータフレーム化。
  2. グラフを描く

風向きと風速はrosavent関数が扱える書式に合うように整形する必要あり。

気象庁のHPの2015年1月1日の鳥取市のデータをXMLパッケージを使って加工

この記事では１日分のデータですが、例えば for文 ＋ rbind 関数を使うと１月分とか１年分（複数年分）のデータの取り込みも可です。
もちろん観測地点のコードを変えると他の観測地のデータも取り込めます。

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
library("XML")
source("http://www.meteobal.com/climatol/rosavent.R")
year<-2015;month<-1;i<-1
#データを入手する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]))
tottori<-t2
#データセットの要約をみる
summary(tottori)
#library(xtable)
#print(xtable(summary(tottori)),"html",include.rownames=F)
levels(tottori$wind_direction)











date time pressure precipitation temperature wind_velocity wind_direction
2015/1/1:24 Min. : 1.00 Min. :1008 Min. :0.000 Min. :-1.3000 Min. : 2.700 西南西 :7
1st Qu.: 6.75 1st Qu.:1009 1st Qu.:0.000 1st Qu.:-0.0250 1st Qu.: 3.950 西 :5
Median :12.50 Median :1009 Median :0.500 Median : 0.6000 Median : 5.300 北西 :4
Mean :12.50 Mean :1011 Mean :1.478 Mean : 0.7333 Mean : 5.417 北北西 :3
3rd Qu.:18.25 3rd Qu.:1013 3rd Qu.:1.750 3rd Qu.: 1.3750 3rd Qu.: 7.075 西北西 :2
Max. :24.00 Max. :1016 Max. :8.500 Max. : 2.7000 Max. :10.100 東南東 :1
NA’s :1 (Other):2

[1] “西” “西南西” “西北西” “東南東” “南西” “南南西” “北西” “北北西”

何度も気象庁のHPにアクセスするのを控えるために作成したデータセットを保存する。

1
2
#ここではRDataとして保存します
save("tottori", file="/tottori.RData")

グラフを描きます。
まず、気温、気圧、降水量

1
2
3
4
5
6
7
8
9
10
11
12
13
#期間を指定しての取り出しを容易にするためxts classにします
library(xts)
t<-strptime(paste(tottori$date, tottori$time), "%Y/%m/%d %H", tz="")
r <- as.POSIXct(round(range(t), "days"))
#気温、気圧、降水量
x.xts<- xts(tottori[,3:5],t)
#png("tottori20150101_1.png",width=1000,height=800)
par(mfrow=c(3,1))
plot.zoo(x.xts[,3],main="気温（℃）",xlab="",ylab="気温（℃）")
plot.zoo(x.xts[,1],main="気圧（hPa）",xlab="",ylab="気圧（hPa）")
plot.zoo(x.xts[,2],main="降水量（mm）",type="h",lend=1,xlab="",ylab="降水量（mm）")
par(mfrow=c(1,1))
#dev.off()

風向き、風速

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
wind<-data.frame(direction=as.vector(chartr("東西南北", "EWSN",t2$wind_direction)),velocity=as.numeric(as.vector(t2$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")
windrose<-data.frame(matrix(rep(0,16*5),ncol=16))
colnames(windrose)<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
rownames(windrose)<-c("[0,3]","(3,6]","(6,10]","(10,20] ","(20, ")
for ( i in 1:nrow(wind)){
	if (wind$velocity[i]<=3) {
	windrose[1,as.vector(wind$direction)[i]]<-windrose[1,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>3 & wind$velocity[i]<=6) {
	windrose[2,as.vector(wind$direction)[i]]<-windrose[2,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>6 & wind$velocity[i]<=10) {
	windrose[3,as.vector(wind$direction)[i]]<-windrose[3,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>10 & wind$velocity[i]<=20) {
	windrose[4,as.vector(wind$direction)[i]]<-windrose[4,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>20) {
	windrose[5,as.vector(wind$direction)[i]]<-windrose[5,as.vector(wind$direction)[i]]+1
	}
	}
#風配図
library(knitr)
kable(windrose)
#library(climatol)
#png("tottori20150101_2.png",width=800,height=800)
rosavent(windrose,5,5,ang=-3*pi/16,main="windrose(Tottori 2015/1/1)")
#dev.off()
N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
[0,3] 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0
(3,6] 0 0 0 0 0 1 0 0 0 0 1 5 1 1 3 1
(6,10] 0 0 0 0 0 0 0 0 0 0 0 2 3 0 1 2
(10,20] 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0
(20, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

色を変更

1
2
3
#png("tottori20150101_3.png",width=800,height=800)
rosavent(windrose,5,5,ang=-3*pi/16,main="windrose(Tottori 2015/1/1)",col=c("azure","cyan","yellow","magenta","red"))
#dev.off()

R+XMLパッケージを使って鳥取市の2000年から2015年4月までの気象データを入手し、加工。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
tottori <- read.table("tottori2000_2014.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE,fileEncoding="cp932")
tottori2015 <- read.table("tottori201501_04.csv", header=TRUE, sep=",", na.strings="NA", dec=".", strip.white=TRUE,fileEncoding="cp932")
tottori<-rbind(tottori,tottori2015)
library(xts)
source("http://www.meteobal.com/climatol/rosavent.R")
t<-strptime(paste(tottori$date, tottori$time), "%Y/%m/%d %H", tz="")
r <- as.POSIXct(round(range(t), "days"))
#気温、気圧、降水量
x.xts<- xts(tottori[,3:5],t)
#plot.zoo
par(mfrow=c(3,1))
plot.zoo(x.xts[,3],main="気温（℃）")
plot.zoo(x.xts[,1],main="気圧（hPa）")
plot.zoo(x.xts[,2],main="降水量（mm）",type="h",lend=1)
par(mfrow=c(1,1))
#各月ごとの平均値
par(mfrow=c(3,1))
plot.zoo(apply.monthly(x.xts[,3],mean,na.rm=T),main="気温（℃）")
plot.zoo(apply.monthly(x.xts[,1],mean,na.rm=T),main="気圧（hPa）")
plot.zoo(apply.monthly(x.xts[,2],sum,na.rm=T),main="降水量（mm）",type="h",lend=1)
par(mfrow=c(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
#風向き、風速
wind.xts<- xts(tottori[,7:6],t)
#2014/3～2014/5,2014/6～2014/8,2014/9～2014/11,2014/12～2015/2
#(例) 2014/12～2015/2
period<-"2014-12::2015-02"
wind<-data.frame(direction=as.vector(chartr("東西南北", "EWSN",wind.xts[period,1])),velocity=as.numeric(as.vector(wind.xts[period,2])))
table(wind$velocity,wind$direction)
#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")
windrose<-data.frame(matrix(rep(0,16*5),ncol=16))
colnames(windrose)<-c("N","NNE","NE","ENE","E","ESE","SE","SSE","S","SSW","SW","WSW","W","WNW","NW","NNW")
rownames(windrose)<-c("[0,3]","(3,6]","(6,10]","(10,20] ","(20, ")
for ( i in 1:nrow(wind)){
	if (wind$velocity[i]<=3) {
	windrose[1,as.vector(wind$direction)[i]]<-windrose[1,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>3 & wind$velocity[i]<=6) {
	windrose[2,as.vector(wind$direction)[i]]<-windrose[2,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>6 & wind$velocity[i]<=10) {
	windrose[3,as.vector(wind$direction)[i]]<-windrose[3,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>10 & wind$velocity[i]<=20) {
	windrose[4,as.vector(wind$direction)[i]]<-windrose[4,as.vector(wind$direction)[i]]+1
	} else 	if (wind$velocity[i]>20) {
	windrose[5,as.vector(wind$direction)[i]]<-windrose[5,as.vector(wind$direction)[i]]+1
	}
	}
#風配図
library(xtable)
print(xtable(windrose,digits=0),"html")
#library(climatol)
#source("http://www.meteobal.com/climatol/rosavent.R")
rosavent(windrose,5,5,ang=-3*pi/16,main=paste("風配図(",period,")"),col=c("azure","cyan","yellow","magenta","red") )








N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
[0,3] 55 33 27 31 84 416 134 31 24 19 14 31 23 35 48 49
(3,6] 72 26 11 4 2 149 69 55 60 28 6 6 17 40 109 193
(6,10] 19 4 3 0 0 1 2 9 43 6 0 0 2 26 16 38
(10,20] 1 0 0 0 0 0 0 3 10 0 0 0 0 0 0 2
(20, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0








N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
[0,3] 102 66 41 49 129 417 177 61 50 35 32 31 45 42 68 108
(3,6] 52 11 9 8 4 94 69 67 132 49 6 5 16 47 70 186
(6,10] 0 1 2 1 0 0 3 6 10 5 0 0 0 7 3 2
(10,20] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0
(20, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0








N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
[0,3] 62 69 48 40 128 606 185 41 30 28 26 23 21 21 34 43
(3,6] 99 61 45 4 16 154 42 42 33 18 9 14 19 13 22 95
(6,10] 13 10 8 0 0 0 0 5 4 1 0 7 9 4 1 21
(10,20] 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3
(20, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0








N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW
[0,3] 27 22 14 14 57 352 159 57 36 41 27 37 40 14 23 19
(3,6] 61 19 5 4 9 166 107 45 51 32 33 87 73 54 61 94
(6,10] 17 3 0 0 0 0 2 10 18 0 0 46 65 27 32 83
(10,20] 1 0 0 0 0 0 0 0 1 0 0 4 7 1 0 0
(20, 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0