2016年11月2日水曜日

Rで順列、組み合わせを計算する

最近知ったんですが、Rで順列(Permutation)や、組み合わせ(Combination)を計算するデフォルトの関数ってないんですね。

【2016年11月18日追記 開始】
勘違いしていました。訂正します。
組み合せの数を計算するchooseという関数があります。
例えば、4個から2個選ぶ組み合わせの数は

> choose(4,2)
[1] 6

みたいな感じで求めることができます。
【2016年11月18日追記 終了】

そういうパッケージはあるみたいですけど、わざわざパッケージを使うまでもないような気もしないでもないです。

高校数学を復習する意味でも、自分で実装してみましょう。

まず、順列。

n個の中からk個を取り出して、順番を意識して並べていくような場合の数ですね。

nから初めて、1ずつ減らしながら、r個分かければ計算できます。

    n × (n-1) × (n-2) × ... × {n-(r-2)} × {n-(r-1)}

上記をやりたい場合、受け取った引数を全部かけてくれるprodという関数が使えます。1ずつ減っていく引数は、 ":" (コロン)を使えばよさそうですね。

ということで、

    prod(n:(n-r+1))

とやればいいのですが、これだと r が 0 のときにうまくいきません。0個を選び出す順列の数は 1 なので、関数にしてしまって、↓こんな感じで、どうでしょうか?

# 順列を計算する関数定義
Permu <- function(n, r){
  ifelse( r==0, 1, prod(n:(n-r+1)) )
}


動作確認してみましょう。

> Permu(4, 0)
[1] 1
> Permu(4, 1)
[1] 4
> Permu(4, 2)
[1] 12
> Permu(4, 3)
[1] 24
> Permu(4, 4)

[1] 24

うん、よさそうです。

次は、組み合わせ。

nからr個選ぶ順列を、rの階乗で割れば(r個の順番を無視する)、OKですね。

階乗はfactorialという関数があるので、そのまま使えます。

こちらも関数定義にしてみましょう。

# 組み合わせを計算する関数定義
Combi <- function(n, r){
  Permu(n, r)/factorial(r)
}


こちらも動作確認してみましょう。

> Combi(4, 0)
[1] 1
> Combi(4, 1)
[1] 4
> Combi(4, 2)
[1] 6
> Combi(4, 3)
[1] 4
> Combi(4, 4)

[1] 1

うん、よさそうです。




2016年6月16日木曜日

Rでデータフレームを部分抽出後、ベクトルや行列に変換する

例えば↓こんな感じで、各市町村の人口データがあったとします。

city   <- c("A市", "B町", "C村", "D市")
male   <- c(95, 49, 26, 45)
female <- c(99, 51, 24, 55)
df <- data.frame(city, male, female)
df
  city male female
1  A市   95     99
2  B町   49     51
3  C村   26     24
4  D市   45     55

で、人口や男女比を比較するために、それぞれの市町村に対してグラフを描いてみようとしたとします。

じゃあ、1行ごとにループで回しながら、barplotあたりを使って、グラフを描いてみようとしまして、↓こんな風にスクリプトを書くと(iがループのカウンタだと思ってね)、

# i行目の男女人口を棒グラフにしたい・・・
barplot(df[i, c("male", "female")])

 barplot.default(df[i, c("male", "female")]) でエラー: 
   'height' はベクトルか行列でなければなりません

なんでだ?と思って、切り出される部分をチェックしてみると、

df[i, c("male", "female")]

  male female
1   95     99

class(df[i, c("male", "female")])

[1] "data.frame"

構成要素としては、2つの数値の並びになっているもの、型がデータフレームなもんだから、barplotにそのまま渡すとエラーになってしまうのね。ちなみに、円グラフのpie関数なんかでも同じで、データフレームのままでは受け取ってくれない。

じゃあ、ベクトルに変換してあげよう。

barplot(as.vector(df[i, c("male", "female")], mode="numeric"))

as.vector後にbarplotに渡したもの

ラベル(maleとかfemale)の情報は消えちゃったけど、まあなんとかなったと。

エラーメッセージには「ベクトルか行列でなければなりません」とあったから、今度は行列にしてみよう↓


barplot( as.matrix(df[i, c("male", "female")]) )



as.matrix後にbarplotに渡したもの

列名が残っているおかげで、軸ラベルがそのまま出てくれるので、こっちの方がいいですね。
ちなみに↓こんな風に書いても、同じ動きになります。

barplot( t(t(df[1, c("male", "female")])) )


転置したときに行列に変換されるからですね。文字数的には少なく済みますが、あとで見たときや、他の人が見たときに混乱するかもしれないので、こういう書き方はやめた方がよさそうです。
で、ここからは補足になるんですが、もともとデータフレームを入力とするggplotを使ってみたらどうだろうかと。

前述のデータフレームはwideフォーマットなので、いったんlongフォーマットに変換します。

library(reshape2)

# longフォーマットに変換
df_l <- melt( df, id.vars="city", variable.name="gender",
              value.name="population" )
df_l
  city gender population
1  A市   male         95
2  B町   male         49
3  C村   male         26
4  D市   male         45
5  A市 female         99
6  B町 female         51
7  C村 female         24
8  D市 female         55

「dl_l」という名前のロングフォーマットのデータフレームに変換できたので、これをggplotに渡します。

library(ggplot2)


# A市の行だけ取り出してプロット

ggplot( df_l[city=="A市", ], aes(x=gender, y=population)) +
  geom_bar(stat="identity")



1行だけ取り出してggplot
うん、想定通りに表示できました。

じゃあ、これをループで回す? いやいや、ggplotなら、いっぺんにプロットできますよね。

ggplot( df_l, aes(x=city, y=population, fill=gender)) +
  geom_bar(stat="identity", position="dodge")



各市町村の男女人口や比率を比較したいなら、この可視化がよさげですね。

でも、このデフォルトのmaleとfemaleの色の割り当ては、無用な勘違いを生んでしまいそうです。
ggplotで、塗りつぶし(fill)の色を変更する方法については、

Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集
Rグラフィックスクックブック ―ggplot2によるグラフ作成のレシピ集

↑この書籍の「レシピ3.2 棒をグループ化する」のところなんかに、カラーパレットとしてbrewerパレットを使うレシピが載っていたりします。

その方法を使ってみると↓


はい、それっぽい色になりました。

やり方については、上記の書籍を買って読んでみてください。(最後はアフィリエイトかよ!)




2016年3月30日水曜日

Rを使って福岡市の人口密度ランキング(町丁目単位)

福岡県の市単位(政令指定都市は区単位)の人口密度単位のランキングは、↓こちらにあって、

福岡県の人口密度番付 - 都道府県・市区町村ランキング【日本・地域番付】

  1位 福岡市 中央区(11770人/㎢)
  2位 福岡市 城南区( 8031人/㎢)
  3位 福岡市 南区 ( 7976人/㎢)

となっております。トップは予想通りの福岡市中央区。

ただ、中央区と言ってもそれなりに広いし、人口密度にはムラがあるでしょうし、もしかしたら他の区にもギュギュっと人の密集しているエリアがあるかもしれません。

ということで、Rを使って町丁目単位(博多駅前1丁目とか、天神2丁目とか)でのランキングを出してみました。

とりあえず結果をどうぞ↓



中央区全体の人口密度は 11770人/㎢ でしたが、中央区荒戸2丁目で人口密度を計算すると 49294人/㎢ と、4倍以上の数字を叩き出しました。

また、上位5位の荒戸、平尾、薬院、高砂はいずれも中央区の町名ですが、6位に入っている愛宕浜は西区にある町です。

西区と言えば福岡市の中では辺境ですし(注:著者の個人的なイメージです)、前述の市区単位のランキングでは17位とかなり下位にランキングされてました、

実は、この愛宕浜2丁目は室見川をはさんで早良区のすぐ対岸にあって、比較的都市部に近い場所なんですよね。そして、ほとんどマンションしかありません(他には中学校とマルキョウがある)。大抵のマンションは10階以上あるので、おのずと人口密度が高くなるわけですね。

↓Google Earthで見てみるとこんな感じ



では、最後にRスクリプトを載せておきます。

シェープファイルは↓こちらに書いてある方法で入手して、

e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ

↓このやり方でマージしました。

Rで複数のシェープファイルを結合する - Rプログラミングの小ネタ


library(maptools)
library(ggplot2)
 
# シェープファイルを読み込む
pj <- CRS("+proj=longlat +datum=WGS84")
shp <-readShapePoly("merge.shp", proj4string=pj)
 
# 必要な項目を取り出しデータフレームにする
df <- data.frame(町丁名=shp$MOJI, 面積=shp$AREA, 人口=shp$JINKO)
 
# 人口密度の列を追加
df$人口密度 <- df$人口 / (df$面積 / 1000^2) # 1平方キロあたり
 
# ソートして、上位50だけ取り出す
df_sort <- df[order(df$人口密度, decreasing=T), ]
top50 <- df_sort[1:50, ]
 
# 画像ファイルに出力
png("人口密度ランキング.png", width=640, height=1280)
ggplot(top50, aes(x=reorder(町丁名, 人口密度), y=人口密度)) +
  geom_bar(stat="identity") +
  ggtitle("福岡市の人口密度TOP50") +
  ylab("人口密度 (人/㎢)") +
  theme( plot.title=element_text(size=20),
         axis.title.x=element_blank(),
         axis.title.y=element_blank(),
         axis.text.x=element_text(size=20),
         axis.text.y=element_text(size=20)  ) +
  geom_text(aes(label=paste(round(人口密度), "人/㎢"), y=人口密度-5000),
            size=5, colour="white") +
  coord_flip()
dev.off()




2015年12月18日金曜日

Rで複数のシェープファイルを結合する

地図で可視化する際の行政界のシェープは、統計局のものをよく利用させてもらっています。↓以下参考。

e-Stat(統計局)で公開されているShapeファイルを、Rで表示する - Rプログラミングの小ネタ

統計局で公開されているシェープは比較的小さい単位に分かれています。例えば、福岡市だと、区単位(東区、博多区、中央区、南区、西区、城南区、早良区)の7つのシェープに分かれています。

こういうのをまとめて1つにしたい時ってありますよね。

QGISのようなGUIのあるGISソフトでやる方法もありますが、数が多くなるとしんどい。ということで、Rのスクリプトでシェープの結合をする方法を紹介します。

ざっくり言うと、readShapePolyで読み込んで、spRbindで結合させるだけなんですけどね。

でも、そう簡単にはいかない時があるんですよね。というのも、ポリゴンのIDがかぶっているとspRbindでエラーが出てしまう。例えば、東区のシェープの中のあるポリゴンに「1」というIDが付いていて、博多区のシェープの中のあるポリゴンにも「1」というIDが付いていると、この2つを結合させようとする際にエラーになってしまいます。

なので、東区のポリゴンはのIDは、"1.1", "1.2", "1.3", ...、博多区のポリゴンのIDは"2.1", "2.2", "2.3", ... みたいにリネームしたあとに、結合させましょうというのが、↓このスクリプトのポイントでございます。

library(maptools)
 
# シェープのファイル名を取得
shape_files <- list.files(path="shape_folder", pattern=".*\\.shp", full.names=T)
 
# 緯度経度、世界測地系
pj <- CRS("+proj=longlat +datum=WGS84")
 
# ポリゴンのIDが重複しないように前処理したあと、
# リストに格納しておく
shape_list <- list() # 空のリストを用意
 
for ( i in 1:length(shape_files) ) {
  # シェープの読み込み
  shp <-readShapePoly(shape_files[i], proj4string=pj)
 
  # id文字列のベクトルを作る
  ids <- paste(i, 1:length(shp@polygons), sep=".")
 
  # polygonsのIDを変更する
  for ( j in 1:length(shp@polygons) ) {
    shp@polygons[[j]]@ID <- ids[j]
  }
 
  # dataの行名もそれに対応して変更する
  row.names(shp@data) <- ids
 
  # リストに追加する
  shape_list <- c(shape_list, shp)
}
 
# リストに格納してある、シェープを結合する
merged <- shape_list[[1]]
for ( i in 2:length(shape_list) ) {
  merged <- spRbind(merged, shape_list[[i]])
}
 
# シェープの書き出し
writePolyShape(merged, "merge.shp")


↓あとは、おまけみたいなもんです。

ちゃんとくっついているかplotしてみます。ifelseでごちゃごちゃ書いているのは、区ごとに色分けするためです。

ややこしいのが海の部分で、ここを同じ色で塗りたいのですが、区によってMOJIの属性(通常は町丁目名が入っている)に「博多湾」と入っていたり、「博多港」と入っていたり、「水面」と入っていたり、NULLが入っていたりする。博多湾とか博多港は分かるけど、「水面」ってなんだよ。

あと、地図に文字をプロットするようなときには、文字が大きいと重なっちゃうし、小さいとつぶれてしまうことがよくあるので、PDFに出力するとベクタとして書き出せるので吉です。

library(RColorBrewer)
cols <- brewer.pal(8, "Pastel2")
 
merged$col <- ifelse ( merged$MOJI=="博多湾" |
                       merged$MOJI=="博多港" |
                       merged$MOJI=="水面"   |
                       is.na(merged$MOJI)       , cols[1],
              ifelse ( merged$CSS_NAME=="東区"  , cols[2],
              ifelse ( merged$CSS_NAME=="博多区", cols[3],
              ifelse ( merged$CSS_NAME=="中央区", cols[4],
              ifelse ( merged$CSS_NAME=="南区"  , cols[5],
              ifelse ( merged$CSS_NAME=="西区"  , cols[6],
              ifelse ( merged$CSS_NAME=="城南区", cols[7],
              ifelse ( merged$CSS_NAME=="早良区", cols[8], "NG" ))))))))
 
pdf("out.pdf", family="Japan1GothicBBB", width=11.69, height=16.54)
plot(merged, col=merged$col, lwd=0.1)
text(merged$X_CODE, merged$Y_CODE, merged$MOJI, cex=0.1)
dev.off()


↓こんな感じで表示されます。

福岡市の行政界地図

↓PDFファイル自体も置いておきますね。

福岡市の行政界地図(PDFファイル)




2015年11月18日水曜日

Rでベクトル場を図示する

ベクトル場がxy座標の式で与えられているときに、それがどんな場を表しているかは、具体的にいくつかのベクトルを図示してみると、なんとなく見えてきたりします。

  f(x, y) = (x, y)

みたいな単純なベクトル場でさえ、慣れてないとイメージわかなかったりしますよね。

Rで図示してみましょう。(唐突)

#
# ベクトル場を表す関数
#
f <- function(x, y){
  X <- x
  Y <- y
  return( c(X,Y) )
}
 
#
# 矢印を描画する関数
#
write_arrow <- function(x, y){
  a <- 0.1  # 矢の長さ(比率)
  b <- 0.05 # 矢尻の長さ
  arrow_head <- f(x, y) * a   # aの値で長さを調整
 
  arrows( x,                  # 矢のx座標(from)
          y,                  # 矢のy座標(from)
          x + arrow_head[1],  # 矢のx座標(to)
          y + arrow_head[2],  # 矢のy座標(to)
          length = b )         # 矢尻のサイズ
}
 
#
# ここからメイン
#
n <- 10 # 格子の1辺の数
 
# 描画する領域を準備
plot(0, 0, xlim=c(-n, n), ylim=c(-n, n), type="n", xlab="x", ylab="y")
 
# 各格子点について繰り返す
for(x in -n:n){
  for(y in -n:n){
    write_arrow(x, y)
  }
}

f(x, y) = (x, y)のベクトル場


「ぶゎっ」て感じですね。

「ベクトル場を表す関数」のところを書き換えれば、いろいろ図示できます。


  X <- y
  Y <- x
とすると↓
f(x, y) = (y, x)のベクトル場

  X <- -y
  Y <-  x
とすると↓
f(x, y) = (-y, x)のベクトル場

  X <- 5
  Y <- 0
とすると↓
f(x, y) = (5, 0)のベクトル場

  X <- y^2
  Y <- x^2
とすると↓

f(x, y) = (y^2, x^2)のベクトル場

矢が長すぎる場合は、変数aの値を小さくすれば、矢が短くなります。

  a <- 0.01  # 矢の長さ(比率)

としたもの↓

f(x, y) = (y^2, x^2)のベクトル場(矢短め)

矢の長さだけで表現するのに限界を感じたら、色や太さを使ってもいいかもしれません。(色や太さに指定している18やら5やらの数字は、色や太さがほどよくなるように恣意的に決めたものです)

#
# 矢印を描画する関数
#
write_arrow <- function(x, y){
  a <- 0.01  # 矢の長さ(比率)
  b <- 0.05 # 矢尻の長さ
  arrow_head <- f(x, y) * a   # aの値で長さを調整
 
  l <- sqrt(x^2 + y^2)
  cols <- heat.colors(16)
 
  arrows( x,                  # 矢のx座標(from)
          y,                  # 矢のy座標(from)
          x + arrow_head[1],  # 矢のx座標(to)
          y + arrow_head[2],  # 矢のy座標(to)
          length = b,         # 矢尻のサイズ
          col = cols[18 - round(l)], # 矢の色
          lwd = round(l)/5 )         # 矢の太さ
}


f(x, y) = (y^2, x^2)のベクトル場(ベクトルの大きさを、長さ/色/太さで表現)

こうやって、絵に描いたところで問題が解けるというわけでもないですが、なんかこう、気分がのってくるでしょ。




2015年11月12日木曜日

RでFukuoka City Wi-Fiの利用状況を可視化する

オープンデータとして、Fukuoka City Wi-Fiの利用状況のCSVファイルが公開されています↓

Fukuoka City Wi-Fi 利用状況 - データセット - 福岡市のCKAN

が、例によってなかなか扱いにくそうなフォーマットになっております↓

福岡市のオープンデータサイトで公開されていたCSVの形式

1拠点の1日が1行になっており、その行の中に0時台、1時台、2時台、~、23時台のように横に伸びてデータが格納されています。別の日は、別の行に格納されています。

そのままじゃ扱いにくそうなので、日付時刻型で表現して、形もシンプルにしましょう。

d <- read.csv("201410access.csv")
 
d <- d[d$X != "合計", ] # 合計が入っているレコードは削除する
 
df <- data.frame() # 空のデータフレームを用意しておく
 
for ( i in 1:nrow(d) ) {
  tmp  <- paste(d[i, ]$X, sprintf("%02d", 0:23))
  YmdH <- strptime( tmp, "%Y%m%d %H")
 
  value <- as.vector(t(d[i, 4:27]))
 
  df <- rbind(df, data.frame( 年月日時 = YmdH,
                              拠点名   = d[i, ]$拠点名,
                              認証回数 = value ) )
}


そうすると、冗長だけどシンプルな形式になりました↓

CSVの形式を変えてみたもの

ひと月の合計値で棒グラフを書いてみましょう。認証回数の少ない順にソートしています。なんとなくggplotを使ってみました。(いまいち、まだ慣れていない感が・・・)

ggplotのリファレンスはウェブにも秀逸なものがあるので、細かいところはそちらでチェックしてください。

library(ggplot2)
s      <- sort(tapply(df$認証回数, df$拠点名, sum))
df_sum <- data.frame(拠点名=names(s), 認証回数=s)
 
pdf("認証回数合計グラフ.pdf", family="Japan1GothicBBB", width=11.69, height=8.27)
g <- ggplot(df_sum, aes(x=reorder(拠点名, 認証回数), y=認証回数))
g <- g + geom_bar(width=0.8, stat="identity")
g <- g + geom_text(aes(label=認証回数), angle=90, hjust=-0.1, vjust=0.4, size=3) 
g <- g + theme(axis.text.x  = element_text(angle=90, hjust=1, vjust=0.4, size=8))
g <- g + ggtitle("Fukuoka City Wi-Fi 認証回数の合計値(2014年10月)") + xlab("") + ylab("認証回数の合計")
g <- g + ylim(0,115000) # 数字が切れないギリギリ
plot(g)
dev.off()


Fukuoka City Wi-Fi 認証回数の合計値(2014年10月)(PDFファイル)

↓PDFファイルのスクリーンショット
Fukuoka City Wi-Fi ひと月のアクセス数の合計値

やっぱり、博多駅、天神駅が利用者が多いですね。アゴーラ福岡のゼロ(しかも3ヶ月連続)ってのは何だろう。稼動していなかったのか、本当に誰も使わなかっただけなのか。

せっかく日付時刻型にしたので、時系列の推移も見てみましょう。

pdf("1時間ごとの認証回数.pdf", family="Japan1GothicBBB", width=8.27, height=11.69 * 10)
g <- ggplot(df, aes(x=年月日時, y=認証回数))
g <- g + geom_line(colour="black")
g <- g + facet_wrap( ~ 拠点名, ncol=1)
g <- g + ggtitle("Fukuoka City Wi-Fi 1時間ごと認証回数 (2014年10月)")
plot(g)
dev.off()


ggplotだと、拠点ごとにわけて複数のグラフにするってのが簡単にできますね。

1時間ごとの認証回数.pdf(PDFファイル)

↓PDFファイルのスクリーンショット
Fukuoka City Wi-Fi 1時間ごとの利用者数の推移

なんか心電図みたいですね。これが心電図なら、アゴーラ福岡は、ご臨終ですね。

基本1日ごとに山になっているんですが、朝と晩にピークがあって昼間に下がる感じで2つの山っぽく(M字みたいに)なっているものも多いです。そして、博多駅などは土日祝日は山が小さくなってますね。図書館は、平坦になっていることで休館日が分かります。

公開されているのは拠点ごと(=複数のアクセスポイント)のデータですが、APごとにデータを見られれば、APが足りているのか不足気味なのかとか分かりそうです。

アゴーラは撤退してもいいんじゃないでしょうか。そして、うちの集合住宅の目の前にでも付けてもらえると助かります。

っていうか、

Fukuoka City Wi-Fi各拠点の1時間ごとの利用状況(認証回数)のデータです。月1回中旬頃に前月分データを提供。

なんて、書いてあるくせに、去年(2014年)の10月のデータを最後にデータが更新されていません。

作成者  : 福岡市市長室
メンテナー: 広報課

となっているので、これはメンテナーがさぼっているってことでしょうか。そういうことは、やめんてな。

「エントリを駄洒落で締める四十代」




2015年11月5日木曜日

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ファイルのキャプチャ↓



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

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