Rで、福岡市オープンデータサイトのCSVファイルを扱いやすいように加工する

福岡市のオープンデータのサイトに「福岡市の大気環境測定結果（直近48時間）」というCSVデータが公開されています。


Rでこれを取り込みたいときは、

d <- read.csv("http://ckan.open-governmentdata.org/dataset/ef0a812c-538b-4e54-ae3e-b22022a7c83a/resource/c60ccc29-811d-44f8-9aba-dcb660387927/download/kankyodata48.csv")

とするだけでOK。（URLが長くてごちゃごちゃしてますが、read.csvにダブルクォートで指定しているだけです）

でも、中身を見てみると、いまいち使いにくそうなフォーマットなんですよね。

head(d)とやれば分かりますが、ちょっと列が多いので、列名だけ表示させてみます。

> colnames(d)
 [1] "年月日"       "測定局名称"   "測定項目名称" "緯度"       
 [5] "経度"         "単位"         "測定値.1時."  "測定値.2時."
 [9] "測定値.3時."  "測定値.4時."  "測定値.5時."  "測定値.6時."
[13] "測定値.7時."  "測定値.8時."  "測定値.9時."  "測定値.10時."
[17] "測定値.11時." "測定値.12時." "測定値.13時." "測定値.14時."
[21] "測定値.15時." "測定値.16時." "測定値.17時." "測定値.18時."
[25] "測定値.19時." "測定値.20時." "測定値.21時." "測定値.22時."
[29] "測定値.23時." "測定値.24時."

年月日という列があるように、日付が違えば別の行として表現されています。しかし、同一の日の別の時間帯のデータは、

"測定値.1時.", "測定値.2時.", ..., "測定値.24時."

のように列方向（横方向）に伸びています。

1日分だけ参照するならいいのですが、データは直近48時間ですから、大抵の場合3日の日付にまたがっているわけです。なので、48時間分のデータをまとめて参照するのには、使い勝手が悪そう。

ということで、処理しやすいようにCSVの形式を変更するスクリプトを書いてみました。

url <- "http://ckan.open-governmentdata.org/dataset/ef0a812c-538b-4e54-ae3e-b22022a7c83a/resource/c60ccc29-811d-44f8-9aba-dcb660387927/download/kankyodata48.csv"
 
# 測定値に数値じゃないものが混ざっているので
# 因子型にならないように、as.is=Tで読み込む
d <- read.csv(url, as.is=T)
 
df <- data.frame() # 空のデータフレームを用意しておく
 
for ( i in 1:nrow(d) ) {
  # いったん "20151031 01", "20151031 02",... という形にしてから、
  # 日付時刻型に変換する
  tmp <- paste(d[i,]$年月日, sprintf("%02d", 1:24))
  YmdH <- strptime( tmp, "%Y%m%d %H")
 
  value <- as.vector(t(d[i, 7:30])) # 7～30列が、1～24時の値に対応
  value <- as.numeric(value)        # 未測定部分が非数字なのでwarningが出る
 
  df <- rbind(df, data.frame( 年月日時=YmdH,
                              測定局名称=d[i, ]$測定局名称,
                              測定項目名称=d[i, ]$測定項目名称,
                              緯度=d[i, ]$緯度,
                              経度=d[i, ]$経度,
                              単位=d[i, ]$単位,
                              測定値=value) )
}

未設定のデータ部分に「@」が入っているため、warningが出ます。「NA」という共通語があるのになぜ使わない、と愚痴りたいところですが、as.numericしたときにNAに変換されるので、とりあえずwarningはほっといても大丈夫そうです。

生成されたデータフレームは↓こんな形になっています。

> head(df)
             年月日時 測定局名称 測定項目名称     緯度     経度 単位 測定値
1 2015-11-02 01:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb      0
2 2015-11-02 02:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb     NA
3 2015-11-02 03:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb      0
4 2015-11-02 04:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb      0
5 2015-11-02 05:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb      0
6 2015-11-02 06:00:00       香椎   一酸化窒素 33.67239 130.4378  ppb      0

冗長にはなりましたが、扱いやすい形になったと思います。

試しにプロットしてみましょう。

# プロットしたいものだけになるようにフィルタリング
df2 <- df[df$測定局名称=="香椎" & df$測定項目名称=="微小粒子状物質(PM2.5)" , ]
 
plot(df2$年月日時, df2$測定値, xaxt="n", type="l", xlab="日付時刻", ylab="PM2.5")
axis.POSIXct(1, at=seq(df2$年月日時[1],rev(df2$年月日時)[1], by="24 hour"),
             format="%m/%d %H:00")

香椎のPM2.5だけにフィルタリングしてプロット

ええい、まとめてプロットしてしまえ、という人は↓こんな感じで。

#
# 場所と項目の組み合わせに対して、網羅的にグラフを描画
#
pdf("福岡市大気環境測定結果.pdf", family="Japan1GothicBBB", width=8.27, height=11.69)
par(mfrow=c(4,3))
for ( place in unique(df$測定局名称) ) {
  for ( item in unique(df[df$測定局名称==place, ]$測定項目名称) ) {
    tmp <- df[df$測定局名称==place & df$測定項目名称==item , ]
    plot(tmp$年月日時, tmp$測定値, xaxt="n", type="l",
         main=paste(place, "：", item),  xlab="日付時刻", ylab=tmp$単位[1] )
    axis.POSIXct(1, at=seq(tmp$年月日時[1],rev(tmp$年月日時)[1], by="24 hour"),
    format="%m/%d %H:00")
  }
}
par(mfrow=c(1,1))
dev.off()

上記で生成されたPDFファイルはこちら。

PDFファイルのキャプチャ↓

場所と項目の組み合わせに対して、まとめてプロット

複数の測定局間で比較するなら、測定項目ごとに縦軸のスケールを合わせた方がよさそうですね。あと、「風向」はこの可視化じゃダメですね。なんて、ツッコミもあるとは思いますが、データとしてはかなり扱いやすい形になったんじゃないでしょうか。