[R]8月の年代別最高気温の分布をヒストグラムに描いて、閾値を超える日数をカウントする

「日本の夏は暑くなっているのか!?」という話題に釣られて、色々なヒストグラムをRで描いたので、自分へのメモのためにスクリプトやら何やらを貼り付けておきます。

描いたヒストグラムはこんなの。

神戸の8月の最高気温のヒストグラム


期間15年のヒストグラムを縦3×横2=6個並べて表示し、ヒストグラムに猛暑日の閾値(35℃)を赤い線で書き込むとともに、閾値を超えた日数をタイトルに加えるというもの。

まず、元となるデータを気象庁サイトからダウンロードしてきます。ダウンロードしてきたテキストファイルは最初の方が

ダウンロードした時刻:2015/08/11 10:40:47,,,,,
,,,,,
,,,東京,東京,東京
年,月,日,最高気温(℃),最高気温(℃),最高気温(℃)
,,,,,
,,,,品質情報,均質番号
1875,8,1,31.9,8,1
1875,8,2,25,8,1
1875,8,3,26.8,8,1

こんな風に色々な情報が書かれていますが、後々の処理のためにテキストエディタでヘッダを書き換えております(Rでヘッダの行数とか列名を指定してあげればいいのですが、めんどうだったのでついw)。

y,m,d,max,a,b
1875,8,1,31.9,8,1
1875,8,2,25,8,1
1875,8,3,26.8,8,1

と、こんな感じにヘッダを整理しておきます。
また、ファイル名もそのままだと「data.csv」ってファイル名になっていますが、地点名.csv(東京.csv)の形で保存しておきます(複数地点のグラフを簡単にかき分けられるようにするため)。
そして、ヒストグラムを描いたRのスクリプトです。

setwd("C:/ほげほげ/東京最高気温/dat")

#初期設定
city <- "神戸" #観測所名
step <- 15 #1つのヒストグラムの期間
start <- seq(1925,2014,step) #1925~2014年の間でstep間隔の値を発生

#グラフ設定
par(mfrow=c(3,2)) #縦3×横2に並べる
par(oma = c(2, 0, 0, 0)) #余白設定
options(digits=3)
#length(start) …数は6個?

#気温データファイルを開く
Temp <- read.csv(paste(city, ".csv", sep = ""), header=TRUE)

#ヒストグラムの期間繰り返す
for(i in start){
    #気温のy列(年)がヒストグラム対象年のレコードのみを抽出する
    T <- Temp[i <= Temp$y & Temp$y <= i+step-1,]
     
    #猛暑日日数をカウントする(T35["TRUE"]に猛暑日(T$max >=35)日数が入る
    T35 <- table(T$max >= 35)
     
    #それぞれのヒストグラムのタイトル(年と猛暑日日数)
    title <- paste(i, "~", i+step-1, "年 (猛暑日:", T35["TRUE"], "日)", sep="") 
     
    #ヒストグラム作図
    hist(
        T$max,               #最高気温の配列
        col = "#aaaaFF40",   #ヒストグラムの背景色
        border = "#0000AA",  #ヒストグラムの境界色
        ylim=c(0,120),       #y軸の範囲(0~120)
        breaks=seq(17,40,1), #x軸の範囲と間隔(17(℃)~40(℃)、1℃刻み)
        main=title,          #ヒストグラムのタイトル
        xlab="最高気温[℃]", #X軸のラベル
        ylab="日数"          #Y軸のラベル
    )
    #猛暑日の閾値:35℃に赤い線を引く
    lines(c(35,35), c(0,120), col="red", xlim=c(17, 40))
    
    #日数、平均値、中央値を凡例で表示
    leg <- c(paste("n =", length(T$max)))
    leg <- c(leg, paste("平均値 =", round(mean(T$max), digits = 2)))
    leg <- c(leg, paste("中央値 =", round(median(T$max), digits = 2)))
    leg <- c(leg, paste("SD=", round(sd(T$max), digits = 2)))
    legend("topleft", legend = leg, cex = 0.9, bty="n")
}
 
#全体のタイトル
title <- paste(city,"の8月の最高気温")
mtext(title, outer=TRUE, side=1, cex=1.5, line = 0.2)

<謝辞>
最後に、データを提供してくださった気象庁、「R」の関係者の皆様、Twitterで色々とアドバイスを下さった皆様にお礼申し上げます。


関係記事・参考資料等
R-Tips「55. 複数の図」
グラフの作り方講座にげんめっ!(引き続き東京は暑くなっているのか問題を例に)-Togetterまとめ
過去の気象データ・ダウンロード-気象庁