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



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

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






Rで、変数の値を変数名や列名にする

タイトルの意味分かりました?

例えば、x00, x01, x02, ..., x99 という名前の100個の変数を作りたいとか、それらに値を設定したいというときに使える方法です。(そんなシチュエーションがあるかどうかは別として)

上記の変数に、0, 1, 2, ..., 99 の値を設定したいとき、

x00 <- 0
x01 <- 1
x02 <- 2
...
x99 <- 99


を全部書けばいいわけですが、それは大変なので、下記のようにできます。

for(i in 0:99){
  # 命令の文字列を作る
  s <- paste("x", sprintf("%02d",i), " <- ", i, sep="")

  # 文字列を命令として実行する
  eval(parse(text=s))
}

ls() # x00~x99の100個の変数ができていることを確認



このevalとparseのコンビを、列名に使いたい場合は、↓こんな感じでどうでしょうか。

ちなみにRには"letters"なる定数が用意されていて、その中身にはaからzまでの文字が格納されています。

> letters
 [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j"
[11] "k" "l" "m" "n" "o" "p" "q" "r" "s" "t"
[21] "u" "v" "w" "x" "y" "z"

これらを列名にもつ、データフレームを作ってみましょう。

> s1 <- paste(letters, rep("=1:10", 26), sep="")
> s2 <- paste(s1, collapse=",")
> s3 <- paste("df <- data.frame(", s2, ")" )
> eval(parse(text=s3))
> df

    a  b  c ...  x  y  z
1   1  1  1 ...  1  1  1
2   2  2  2 ...  2  2  2
3   3  3  3 ...  3  3  3
4   4  4  4 ...  4  4  4
5   5  5  5 ...  5  5  5
6   6  6  6 ...  6  6  6
7   7  7  7 ...  7  7  7
8   8  8  8 ...  8  8  8
9   9  9  9 ...  9  9  9
10 10 10 10 ... 10 10 10

ダブルクォートの何に予約語が入ると、何がなんだか分からなくなりますね。IDEの強調表示とかで乗り切りましょう。

正直あまりエレガントな方法ではない気もしますが、文字列操作さえできればなんでもできるという、根拠のない自信を持つことはできました。





Rで10進数→16進数、16進数→10進数の変換

10進数から16進数に変換する方法の1つとして、sprintfを使う方法があります。

> sprintf("%X", 0:31) 
 [1] "0"  "1"  "2"  "3"  "4"  "5"  "6"  "7"
 [9] "8"  "9"  "A"  "B"  "C"  "D"  "E"  "F"
[17] "10" "11" "12" "13" "14" "15" "16" "17"
[25] "18" "19" "1A" "1B" "1C" "1D" "1E" "1F"


16進といえば、色のRGB値を指定するときに使ったりしますが、その場合は "7" とか "F" ではなく "07" とか "0F" とかの形で欲しいですよね。

そういうときは、

> sprintf("%02X", 0:31)
 [1] "00" "01" "02" "03" "04" "05" "06" "07"
 [9] "08" "09" "0A" "0B" "0C" "0D" "0E" "0F"
[17] "10" "11" "12" "13" "14" "15" "16" "17"
[25] "18" "19" "1A" "1B" "1C" "1D" "1E" "1F"


とすればOK。

sprintfの書式を知っている人には蛇足ですが、「%02X」の意味は、
  "X": 16進で
  "2": 2桁で
  "0": 桁が足りないときは0詰めで
という指示になります。

16進数でのRGB指定を使ってちょっと遊んでみましょう。

# まず枠だけプロット
plot( 0, 0, xlim=c(0,255), ylim=c(0,255),
      type="n", xlab="R", ylab="G" )

# 16進の値を取り出すためのテーブル
h <- sprintf("%02X", 0:255)

for(r in 1:256) {
  for(g in 1:256) {
    col <- paste("#", h[r], h[g], "00", sep="")
    points(r, g,  col=col)
  }
}

colオプションでRとGを変化させながらプロット

青(B)を固定して、赤(R)と緑(G)を変化させたグラデーションができました。

また、10進から16進への変換を行う方法として、as.hexmode関数を使うというのもあります。

> as.hexmode(31)
[1] "1f"

表示結果は文字列(character)っぽいですが、実は違っていて、hexmodeという16進の型になっているみたいです。

> class(as.hexmode(31))
[1] "hexmode"

なので、

> a <- as.hexmode(7)
> a
[1] "7"

> b <- as.hexmode(8)
> b
[1] "8"

> a + b
[1] "f"

> a * b
[1] "38"

のように、16進数として演算も行われます。

この この "38"(さん はち)が、10進での56(7×8の結果)になっているかどうかは、下記のようにして確認できます。

> 0x38
[1] 56

16進から10進への変換は、「0x」を付けるだけで行えるんですね。




2015年11月1日日曜日

Rを使って郵便番号界の地図を作る(統計局のShape+日本郵便のCSV)

郵便番号がどんな感じで割り当てられているか、地図上で確認したかったんです。

GPSで取った地点情報がプライバシーに配慮云々で郵便番号になっていた、なんてことがありまして、その粒度がどんなものか確認したかったってのが発端です。

検索してみても無償のものは意外に見つからないもので、国際航業さんが売っていたりします↓

PAREA-Zip郵便番号界 | 国際航業株式会社

ノースリーブな(=無い袖は振れない)わたくしとしては、自分で作るしかないわけで。

郵便番号のリストは↓ここにありました。

郵便番号データダウンロード - 日本郵便

住所文字列と郵便番号が入っていればどれでもいいんですが、とりあえず、「住所の郵便番号(CSV形式)」の「読み仮名データの促音・拗音を小書きで表記するもの」というのをダウンロードしました。

次は地図です。こういうときに使う境界データとしては、統計局のe-StatでShapeファイルが公開されているので、これが便利。入手の方法は下記を参照。

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

今回は、我らがサラリーマンの聖地 新橋を有する港区を対象にしてみました。

とりあえず地図をプロットするところまでやってみましょう。

library(maptools)     # GIS系のパッケージ
# 座標系の情報を持つオブジェクトを作る
# longlatは緯度経度指定、WGS84は世界測地系であることを表す
pj <- CRS("+proj=longlat +datum=WGS84")

# シェープファイルの読み込み
map <- readShapePoly("h22ka13103.shp", proj4string=pj)

# とりあえずプロット
plot(map)

shapeをプロットしただけのもの

この地図上に追加で郵便番号をプロットすればよさそうですね。

読み込んだシェープファイルにはポリゴンの境界情報だけでなく、それぞれのポリゴンの属性の情報も含まれています。

こうやれば確認できます

head(map@data)

       AREA  (略)  KEN_NAME SITYO_NAME GST_NAME  (略)         MOJI
0 729237.30          東京都       <NA>     港区       元赤坂2丁目
1 148246.50          東京都       <NA>     港区       北青山1丁目
2  89514.52          東京都       <NA>     港区       元赤坂1丁目
3  93129.48          東京都       <NA>     港区         赤坂3丁目
4 151487.30          東京都       <NA>     港区         赤坂4丁目
5 199895.30          東京都       <NA>     港区       北青山2丁目

「MOJI」っていう列に住所が入っているようです。これだけにアクセスしたい場合は、「map$MOJI」と打てばOK。厳密には「map@data$MOJI」なんでしょうけど、省略できるようです。

これを郵便局のCSVファイルとマッチさせればよさそう。

# 郵便番号データの読み込み
post_data <- read.csv("13TOKYO.CSV", head=F)
post_data <- post_data[post_data$V8=="港区", ] # 港区だけ取り出し

# ○丁目の部分をカットして、MOJI2の列に格納
map$MOJI2 <- sub("[0-9].*", "", map$MOJI)

# MOJI2と郵便番号データの文字列がマッチしたら、
# 郵便番号を取り出して設定する
for ( i in 1:nrow(map@data) ) {
  # 複数マッチしたら、安直に1つ目を採用
  matched_index    <- grep(map$MOJI2[i], post_data$V9)[1]
  map$POST[i] <- post_data$V3[matched_index]
}

plot(map)
text( coordinates(map), as.character(map$POST), cex=0.5)

郵便番号をtext関数で追加プロットしたもの

一応、最低限のところまではできました。でも、文字が見にくいとか、地名も欲しいとか、色塗ったほうがキャッチーだよね(?)とか、ということで付け足してみたのが以下。

# ハイフンを入れる
post <- paste ( substr(map$POST, 1, 3), substr(map$POST, 4, 7), sep="-" )

# 地名&郵便番号の文字列を作る
labs <- paste(map$MOJI, post, sep="\n")

# 塗りわけのパレット
library(RColorBrewer) # カラーパレット
cols <- brewer.pal(12, "Set3")
cols <- rep(cols, 10) # 12色じゃ足りないので繰り返し

pdf("郵便番号界マップ(港区).pdf", family="Japan1GothicBBB", width=16.54, height=11.69)

# 郵便番号を因子型にすると、番号ごとに塗りわけできる
plot(map, col=cols[factor(map$POST)], border="gray", lwd=0.1)
text( coordinates(map), labs, cex=0.3)

dev.off()

上記のサンプルコードで作ったのが↓これです。

郵便番号界マップ(港区)(PDFファイル)


↓PDFファイルのスクリーンショット
郵便番号に、色と地名を追加したもの

なかなかいい感じになったんじゃないでしょうか。郵便番号ごとに塗り分けるために、因子型にするなんてのは、ちょっとしたtipsです。

が、ツッコミが入る前に白状しておかねばなりますまい。この郵便番号地図、正しくないところがあります。

住所のマッチングが結構甘いです。grepでマッチしたものが複数ある場合、単純に1つ目のものを採用しています。

例えば、↓こういうパターンに対してはうまく機能しています。

"1070052", "赤坂(次のビルを除く)"
"1076090", "赤坂赤坂アークヒルズ・アーク森ビル(地階・階層不明)"
"1076001", "赤坂赤坂アークヒルズ・アーク森ビル(1階)"
・・・

が、↓こんなやつだとうまくいきませんね。

"1050014", "芝(1~3丁目)"
"1080014", "芝(4、5丁目)"

芝1丁目、芝2丁目、・・・、芝5丁目を全部別の行にしてくれていたとしたら、まだプログラム側でなんとかしようなんて気も起きますが、「1~3」とか「4、5」とか、機械的に処理することを全く考慮していないデータなんですもの。今後の日本郵便のがんばりに期待したいところですね。(上から目線)

ということで、こんな間違いだらけのマップじゃ使えないよってかたは、国際航業さんから買ってください。




2015年5月26日火曜日

Rのデータフレームを降順でソートする(decreasing=Tを使わずに)

ほんと、ちょっとした小ネタです。

↓こちらの本を読んでいて知った(気付いた)方法。

RとRubyによるデータ解析入門
RとRubyによるデータ解析入門

まず、サンプル用のデータフレームを作ります。

name <- c("Ana", "Bill", "Chris", "Dan", "Emily")
age  <- c(18, 59, 57, 6, 31)
tbl <- data.frame(name, age)
tbl
 
   name age
1   Ana  18
2  Bill  59
3 Chris  57
4   Dan   6
5 Emily  31


このデータフレームを年齢順にソートしたい場合は、ソート結果のインデックスを返してくれるorder関数を使って、↓こんな風にしますよね。

tbl[order(tbl$age), ]

   name age
4   Dan   6
1   Ana  18
5 Emily  31
3 Chris  57
2  Bill  59


order関数のデフォルトは昇順なので、降順にしたい場合は、↓こうですよね。

tbl[order(tbl$age, decreasing=T), ]

   name age
2  Bill  59
3 Chris  57
5 Emily  31
1   Ana  18
4   Dan   6


でも、"decreasing=T"って打ち込むのって、けっこう面倒じゃないですか? あと、"decreasing"の綴りが思い付かなかったり(英語力・・・)。


【クイズ】マッチ棒を1本足すだけで、逆にしてください

tbl[order(-tbl$age), ]

   name age
2  Bill  59
3 Chris  57
5 Emily  31
1   Ana  18
4   Dan   6


おお、マイナスをかければ、大小関係が逆になりますからね。

でも、可読性は低いと思うので、長い期間使うつもりのスクリプトとか、誰かと共有するスクリプトを書くときは、やめた方が良さそう。「なんで、ここマイナスかけてるの?」とか、聞かれてもめんどくさいし。

なので、個人的でテンポラリな確認作業とかやっているときの、タイピングの省力化などにご活用ください。




2015年4月28日火曜日

Rのhistとggplot(geom_histogram)を比較する

ggplot2パッケージって使ったことありませんでした。記法がとっつきにくいというのもあったのですが、あの見た目を素直に「きれい」とは認めたくなかったというか。「ggplotの方がかっこいい」という人への反発みたいな部分も多分にあるのですが。

「白地に線画」みたいなビジュアルが好きなんです。下手に色をつけたらかえって趣味が悪くなったりするじゃないですか。自分自身センスがないことを自覚しているので、色やら塗りつぶしやら飾りっぽいものやらは極力使わないようにしていたんですよね。

でも今回、ggplot2パッケージを使ってみようと思ったのは、ベースとなるレイヤーを用意して、次に重ねたいレイヤーを(場合によっては複数)用意して、最後にバンっとplotする、っていうコンセプトがなんだか使いやすいのではないかという気がしたから。なので、徐々に試していこうかなと思っています。

前置きが長くなりました。

今回はヒストグラムですが、多くの人に馴染みのあるデフォルトのhistでの描画と、ggplot2パッケージを使った場合とで、どんな風に違うのか試してみました。

正規分布に従う乱数を発生させて、普通にヒストグラムを描いてみます。


data1 <- rnorm(1000)
 
hist(data1)

通常のhist関数で描画

ggplot2パッケージを使うと↓こんな描き方になります。

install.packages("ggplot2")
library(ggplot2)
 
df1 <- data.frame(data1) # データフレームにしておく必要がある
 
g <- ggplot(df1, aes(x=data1)) # 対象がdata1列であることを指定
g <- g + geom_histogram()      # ヒストグラムを描画することを指定
plot(g)                        # 描画

ggplot2パッケージを使ったもの

先ほどの、デフォルトのhistと階級幅を合わせる(=0.5)には、binwidthを指定します↓

g <- ggplot(df1, aes(x=data1))
g <- g + geom_histogram(binwidth=0.5)
plot(g)

ggplot2(階級幅を指定)

ちなみにまとめて書いてしまう方法もあって、これくらいレイヤー数だったら、この方がシンプルかもしれません。

ggplot(df1, aes(x=data1)) + geom_histogram(binwidth=0.5)


で、見た目どうでしょうか、黒のベタ塗りはけっこう好みが分かれそうです。

でも、もちろん、指定すればどうにでもできます。

g <- ggplot(df1, aes(x=data1))
g <- g + geom_histogram(colour="black", fill="white", binwidth=0.5)
plot(g)

ggplot2(ボーダーと塗りつぶしの色を指定)

colourで枠線の色を、fillで塗りつぶしの色を指定しました。私はこっちの方が好きですね。




2015年4月17日金曜日

Rでブートストラップ法

「ブートストラップ」という言葉は聞いたことはあるけど、具体的にどんなことなのかはよく知らない、というのが私のレベルです。

が、「統計学入門」(東京大学教養学部統計学教室 編)という教科書↓

統計学入門 (基礎統計学)
統計学入門 (基礎統計学)

の中の演習問題に出てきた(でも、教科書の本編には解説はない)ので、Rでちょっとやってみました。

第3章 練習問題
3.4 <ブートストラップ> 以下の手法をパソコンで行え。

i)  [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

ii)  i)を11回繰返して実行し、11個のランダムな番号を得たのち、表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、相関係数rを計算せよ。

iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップという。

「表3.1、図3.1」というのは、教科書の本編のところに出てきた表と図で、相関関係が認められる例として紹介されていたものです。

統計学入門 (基礎統計学) より

統計学入門 (基礎統計学) より

このデータを使う必要があるので、まずは入力しましょう。

> bro <- c(71, 68, 66, 67, 70, 71, 70, 73, 72, 65, 66)
> sis <- c(69, 64, 65, 63, 65, 62, 65, 64, 66, 59, 62)

書籍に載っていた図と同じような感じでプロットしてみました↓

> plot(bro, sis, xlim=c(68,70), ylim=c(59,70), asp=1, pch=16)

書籍に載っていたデータをRでプロットしたもの

書籍に出てくる図とほぼ同じものができました。ただ、点が10個しかなく1個少ないです。これは(70, 65)という全く同じデータが2つあって同じ座標にプロットされているためです。書籍の方では少しずらして、分かるようになっていました。

あと、書籍だと(64,62)あたりに点がありますが(一番左側のもの)、これは間違いでしょうね。点を数えてみると12個あるし・・・

使うデータは入力したので、演習問題をやっていきましょう。

i)  [1, 11]に属する、整数の乱数を発生させる手続きを考えよ。

乱数の従う分布について特に言及していないので、一様分布でしょうね。

まず思いつくのは、

> ceiling(runif(1, 0, 11))   # 1~11の整数の乱数

runif関数が返す値をceiling関数で整数化してやる感じですね。

でも、Rで整数の乱数を発生させるんだったら、sample関数を使った方がシンプルな気がします↓

> sample(1:11, 1, replace=T) # 1~11の整数の乱数

1, 2, ..., 11 の中から1個選ぶ、復元抽出(replace=T)で、という感じですね。

ii)  i)を11回繰返して実行し、11個のランダムな番号を得たのち、

これはsample関数の引数を変えるだけでOKですね。

> rnd <- sample(1:11, 11, replace=T)
> rnd
 [1]  8  1  8  2  4 11  5  3  9  8 10

表3.1、図3.1のデータからその番号のデータを取り出し(ただし、同一番号があれば重複して取り出し、そのデータは複数個あると考える)、

番号に応じたデータを取り出すには、先ほど発生させたrndをインデックスに指定ればいいですね。

> bro[rnd]
 [1] 73 71 73 68 67 66 70 66 72 73 65

> sis[rnd]
 [1] 64 69 64 64 63 62 65 65 66 64 59

相関係数rを計算せよ。

> cor(bro[rnd], sis[rnd])
[1] 0.5357521

iii) ii)を200通り繰返し、相関係数の値 r1, r2, ..., r200 を得たとき、そのヒストグラムを作れ。以上の統計手法をブートストラップ(太字)という。

> r <- numeric(200) # 200回分の格納用
>
> for ( i in 1:200 ) {
+   rnd <- sample(1:11, 11, replace=T)
+   r[i] <- cor(bro[rnd], sis[rnd])
+ }
>
> hist(r)

ブートストラップ法で作成したヒストグラム

これで宿題は終わりました。単位は取れそうです。

元のデータの相関係数をそのまま出すと、

> cor(bro, sis)
[1] 0.5580547

になるのですが、ヒストグラムのモードは違う場所(0.65あたり)に出てますね。分布も左に歪んでいるし。

以上の統計手法をブートストラップという。

とのことですが、ヒストグラムの意味(モードのずれや歪み)や使い方については、もうちょっと勉強しないといけなさそうです。