[R]POSIXct(POSIXlt)型を使って特定の月(年・日etc)の値を抽出する

R

Rを使って、特定の月のデータを抽出するコードを書いてみました。特定の月を抽出するためには、日時を表す値から月を抽出する必要があります。このコードを書こうと思い至ったのは、気象庁Webサイトからダウンロードしたデータを使って、最高気温の月別最大値とかを計算してみたり、月ごとの最高気温のヒストグラムを描いてみようと思ったのが発端です。単に抽出ならC#とか他の言語で書いたコードで月ごとにデータをまとめるという方法もあるのですが、できればRだけで済ませたいと思って、今回はRのPOSIXct(POSIXlt)型のmonths関数を使って特定の月のデータを抽出してみたのです。

まず、気象庁のデータは下記のような形式になっています。最初の6行(平均気温は5行)がヘッダになっていて、それ以下の行には日付、値(最高気温)、品質情報、均質情報が記録されています。品質情報とは統計を行うのに十分な品質か否かを表しています。均質情報とは観測環境の変化を表すもので、観測時期の古いものから1から順番に付けられ、観測環境が変化するたびに番号が増加するようです(詳細は気象庁サイトへ)。

ダウンロードした時刻：2020/08/23 16:04:58

,東京,東京,東京
年月日,最高気温(℃),最高気温(℃),最高気温(℃)
,,,
,,品質情報,均質番号
1981/1/1,8.9,8,1
1981/1/2,7.5,8,1
1981/1/3,9.0,8,1
1981/1/4,9.6,8,1
1981/1/5,8.4,8,1

さて、とりあえずコードを載せてみます。7行目でデータを読み込んでいます。その際にヘッダは読み飛ばしています。10行目で品質情報が8のデータのみを抽出しています。

さてさて、ここからが今回のミソとなる部分です。
13行目でas.POSIXctを使って、データの一列目をPOSIXct型に変換しています。as.POSIXctのformatパラメータで「年/月/日」と指定しています。ここの詳細は三重大学・奥村先生のサイト「時系列データ」を参考にさせていただきました。
20行目で一番目のデータ(1981/1/1)の月を整数で取得しています。当初は18行目のmonths関数を使ったのですが、返り値が “1月” で、言語に依存しそうな文字列型でした。これでは使い勝手が非常に悪いので、format関数で文字列にした(19行目)後に、as.integer関数で整数に変換しています。
23行目は指定した月の最高気温を返す関数で、この関数を25～28行目で呼び出して指定した月の最高気温などを計算しています。

setwd("D:/hogehoge/work")

#ヘッダ行数(最低気温・最高気温:6、平均気温:5)
iSkip <- 6

#とりあえず読み込む
dfAll <- read.csv("data.csv",  skip=iSkip, header=F, stringsAsFactors=F)

#正常値の行(品質情報=8)のみを抽出
dfObj <- dfAll[dfAll$V3==8,]

#日付を文字列からPOSIXltに変換
dats <- as.POSIXct(dfObj$V1, format="%Y/%m/%d", tz="Japan")
temp <- dfObj$V2
dfTemp <- data.frame(Dat=dats, Temp=temp)

#1行目の月を取得する
months(dats[1])
format(dats[1],"%m")
as.integer(format(dats[1],"%m"), "%d")

#指定した月の気温を返す関数
fMTmp <- function(x){dfTemp[as.integer(format(dfTemp$Dat,"%m"), "%d")==x,]$Temp}

#1-12月の平均気温の月別最高、最低、平均、標準偏差
t_max <- sapply(1:12, function(x){max(fMTmp(x))})
t_min <- sapply(1:12, function(x){min(fMTmp(x))})
t_mean <- sapply(1:12, function(x){mean(fMTmp(x))})
t_sd <- sapply(1:12, function(x){sd(fMTmp(x))})

参考: 東京の気温の最高・最低気温のヒストグラムを描いてみました。

setwd("D:/hogehoge/work")
#ヘッダ行数(最低気温・最高気温:6、平均気温:5)
iSkip <- 6

#凡例の場所
sPos <- c("topright","topright","topright","topright","topleft","topleft","topleft","topleft","topleft","topleft","topright","topright")

#グラフ設定
par(mfcol=c(3,2))      #複数グラフ(縦×横に並べる)
par(mar=c(3, 3, 1.5, 1)) #余白設定
par(mgp=c(2, 0.7, 0))  #グラフ枠と軸との間隔設定
par(xaxs = "i")        #データ範囲とグラフの範囲(X軸)
par(yaxs = "i")        #データ範囲とグラフの範囲(Y軸)

dfAll <- read.csv("data.csv",  skip=iSkip, header=F, stringsAsFactors=F)

#正常値の行(品質情報=8)のみを抽出、日付を文字列からPOSIXltに変換
dats <- as.POSIXct(dfAll[dfAll$V3==8,]$V1, format="%Y/%m/%d", tz = "Japan")
temp <- dfAll[dfAll$V3==8,]$V2

#データフレームに
dfTemp <- data.frame(Dat=dats, Temp=temp)

for(iM in 1:6){
 temp_m <- dfTemp[as.integer(format(dfTemp$Dat,"%m"), "%d")==iM,]$Temp
 hist(temp_m,
  breaks=seq(-10,40,1), xlim=c(0,40),
  main=sprintf("%d月", iM), xlab="最高気温[℃]", ylab="n(日数)"
 )
 
 legend(sPos[iM], bty="n",
  legend=c(
   sprintf("n=%d", length(temp_m)),
   sprintf("Max=%0.2f", max(temp_m)),
   sprintf("Min=%0.2f", min(temp_m)),
   sprintf("Mean=%0.2f", mean(temp_m)),
   sprintf("SD=%0.2f", sd(temp_m))
  )
 )
}
最高気温(1～6月@東京)

最高気温(1～6月@東京)

最高気温(7～12月@東京)

最高気温(7～12月@東京)

最低気温(1～6月@東京)

最低気温(1～6月@東京)

最低気温(7～12月@東京)

最低気温(7～12月@東京)