オープンデータを解析する- ggplot2を用いたボロノイ分割で厚木市のコンビニ出店を見てみよう！-
オープンデータ
今日もオープンデータを可視化して、何かの役に立ててみようということをやってみたいと思います。昔、同級生が卒業研究で画像認識をやっていて、そのときに「ボロノイ分割」というものがありました。聞いているとですね、いろんな分野に応用されている手法で下のようなものに実際使われているようです。
で、今回はこのボロノイ分割を使って、厚木市の「コンビニ」の所在地を可視化してみようと思います。なぜ厚木市かというと、コンビニのデータが公開されているのって、厚木市しかなくて．．．どこかに落ちてるところがあれば教えてください。
ボロノイ分割とは
ボロノイ分割については、ボロノイ図とはを見ていただくとわかりやすいと思いますが、簡単にいえば、「平面上に、いくつかの点が配置されているとき、その平面上の点を、どの点に最も近いかによって分割してできる図」のことです。イメージがわかりにくいので、下の図を見ていただければ何となくわかってもらえるかと。
今回は、この真ん中の「母点」がコンビニの位置です。で、求めたいのは「ボロノイ境界」ということになります。
データの中身
解析なので、まずはデータから確認します。今回用いたデータはここのページ（厚木市コンビニデータ｜Open data sharing & Download｜LinkData）から落とすことのできる「コンビニデータ」です。このデータを用いて解析を行います。
> head(df) id 住所 緯度 経度 店舗の種類 お店の店名 連絡先 1 神奈川県厚木市三田南２丁目２１番９号 35.47149 139.3483 セブン－イレブン セブンイレブン三田店 046-242-1171 2 神奈川県厚木市妻田南１丁目１７−２５ 35.45299 139.3600 セブン－イレブン セブンイレブン妻田南店 046-223-8455 3 神奈川県厚木市下荻野１０４９−１ 35.48385 139.3415 セブン－イレブン セブンイレブン厚木下荻野店 046-241-5830 4 神奈川県厚木市金田１１９−１ 35.47002 139.3683 セブン－イレブン セブンイレブン厚木金田店 046-225-5572 5 神奈川県厚木市妻田東３丁目２−７ 35.45966 139.3597 セブン－イレブン セブンイレブン厚木妻田店 046-221-0676 6 神奈川県厚木市七沢２３８ 35.44230 139.2991 セブン－イレブン セブンイレブン厚木七沢店 046-248-0531
必要なのは、データの中で「店舗の種類」と「緯度」と「経度」ですから、これらをstore, lat, lonという変数にした次のようなデータフレームを用います。
> head(df)
lat lon store
1 35.47149 139.3483 セブン－イレブン
2 35.45299 139.3600 セブン－イレブン
3 35.48385 139.3415 セブン－イレブン
4 35.47002 139.3683 セブン－イレブン
5 35.45966 139.3597 セブン－イレブン
6 35.44230 139.2991 セブン－イレブン
ボロノイ図を書くために...
さて、ここからは具体的な解析を行っていきますが、いろんなパッケージを使います。
- ggplot2 : 図を書くために使います
- deldir：ボロノイ分割を行うパッケージ
- reshape2：データをggplot2用に変換するのに使います。
- plyr：データの分割・処理・結合 -> 「plyrパッケージで君も前処理スタ☆」改め「plyrパッケージ徹底入門」
- ggmap：地図を書くためのパッケージ
結構、いろんなものに依存しているのでややこしいですがご了承ください。
setwd("~/Desktop/blog/厚木_コンビニ")
library(ggplot2)
library(deldir)
library(reshape2)
library(plyr)
library(ggmap)
#data.frame
df = read.csv("data_conv.csv",sep="\t",header=F)
df = df[,c(3,4,5)]
names(df) = c("lat","lon","store")
#日本語フォントを使えるようにする
quartzFonts(maru=quartzFont(rep("HiraMaruProN-W4",4)),kaku =quartzFont(rep("HiraKakuPro-W6",4)))
#map
#地図の範囲
loc=c(min(df$lon),min(df$lat),max(df$lon),max(df$lat))
#googleから地図データを取得（zoomの大きさを間違えるとこの後エラーがでます。ggplotするところで missing valueとか, 今回は12 or 13）
M=ggmap(get_map(location=loc,zoom=12,source="google"))+xlab("")+ylab("")ここまでで、実は特定の区域を書くことができました。この辺りの流れについては、オープンデータを使ってみよう！-流山市の桜を最適なルートで回ろう！- - Data Science by R and Pythonを参照していただければわかりやすいかと。ここからは、本格的にボロノイ分割を行います。そのためには、まず、先ほど作成したデータフレームdfから、緯度と経度を取ってこないといけません。
############################# #data.frame only lat and lon# ############################# dd <- df[,1:2]
次にhullというものを考えましょう。hullというのは「（果物などの）外皮」のことですが、今回は、データの点集合の一番外側の点ということです。具体的には得られた図を見ていただくのが良いかと思います。
############# #convex hull# ############# hull = chull(dd) hullPointsIndices <- c(hull, hull[1]) #最後に1点加えているのは、点を回って戻ってくる直線を描くため, 下のqplotで意味が分かる。 hullPoly_df <- df[hullPointsIndices,] #上で指定した行だけ抜き出し # plot convex hull qplot(lon, lat, data = dd) + geom_polygon(data = hullPoly_df, colour = 'red', fill = 'red',alpha = I(1/10), size = 1)
外側の点だけ取ってプロットしました。つまり、この範囲にすべての点が入っていることになります。これは何かに使うわけではないですが、図示するときなどには、外枠がある方が見やすいので書いてみました。
次に、ボロノイ分割の計算です。ここで、ボロノイ分割の計算にはdeldirという関数を用いています。その前にrange, expand.rangeという関数を用いていますが、これはコードの中の説明を見てください。簡単にいえば、ボロノイ分割のデータの範囲を絞るために行っています（deldirのrw指定箇所）。
############# #triangulate# ############# #range::データの最大値と最小値を返す #expand_range(x,y):xのうち小さいもののy倍をz=xyとすると、x-z, y+zを計算する #図を書く範囲を限るために計算する xrng <- expand_range(range(dd$lon), .05) yrng <- expand_range(range(dd$lat), .05) #deldir::分割 deldir <- deldir(x=dd$lon,y=dd$lat, rw = c(xrng, yrng))
さて、deldirが返してくる結果には様々なリストが入っています。詳しくは以下のURLから参照していただきたいですが、先ほどの図でいえば、ボロノイ境界のデータが入っています。この他、以下の図（これも先ほどのHPより引用）でいえば、ドロネー図の情報も入っています。
http://cran.r-project.org/web/packages/deldir/deldir.pdf
さて、ドロネー図というものを書いてみましょう。
#################
#Delaunay fugure#
#################
qplot(lon, lat, data = dd) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .5,data = deldir$delsgs, colour = 'red', alpha = .5)
############################
# triangulate with polygons#
############################
trilist <- unclass(triang.list(deldir)) #クラス構造を分解します
tridf <- melt(trilist, id.vars = c('x','y'))[,c(3,1:2)] #ggplot用に加工します, meltは便利！
names(tridf) <- c('tri','x','y')
qplot(lon, lat, data = dd) + geom_polygon(aes(x=x, y=y, fill=factor(tri)), data = tridf, colour = 'black', alpha = .5) + scale_fill_discrete(guide = FALSE)そのドロネー図がこちらになります。この図が、何の役に立つのかは後ほどわかります。
次に、この図に色を付けてみます。これも、何の意味があるかは後ほどです。
さて、ここで使っているgeom_polygonという関数が、空間を色で塗りつぶすための関数です。これを使うとボロノイ図が奇麗に書けます。その、ボロノイ図を次に書いていきましょう。ボロノイ図を書くためのデータも先ほどdeldirで作成したデータの中に入っています。上のものと同じようにして、以下のコマンドを打ってやれば．．．（今回は後ろ側に地図を入れます！）
######### #voronoi# ######### #分割図のみを作成する（色をぬったりしない） M + geom_point(aes(x=lon, y=lat,colour=factor(store)),size=3, data = df) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25,data = deldir$dirsgs, linetype = 2)
図は、こんな感じになります
このままじゃ寂しいので、領域に色をぬりましょう！ということをするためには、geom_polygonです。書いてみると、次のようになります。
#######################
#voronoi with polygons#
#######################
tilelist <- unclass(tile.list(deldir))
tilelist <- lapply(tilelist, function(l){
data.frame(x = l$x, y = l$y)
})
#melt(tilelist,id.var=c("x","y"))がうまくいかないので....
tiledf <- melt(tilelist, id.vars = c("x"))
tiledf <- tiledf[,c(4,1,3)]
names(tiledf) <- c('tile','x','y')
M + geom_polygon(aes(x = x, y = y, fill = factor(tile)),data = tiledf, colour = 'black', alpha = .3) + geom_point(aes(x = lon, y = lat,colour=factor(store)),size=2, data = df)+theme_set(theme_bw(base_size = 12,base_family="HiraMaruProN-W4"))結構奇麗に書けます。ただ問題がありまして．．．解消できてないんですけど、領域の色をコンビニの種類で分けたいんですが、どうもgeom_polygonでは無理かもしれなくて．．．解消法がわかる方教えていただけると嬉しいです。
図からわかりそうなこと。
先ほどお見せしたこの図に戻ってみましょう。すると、すぐわかることとしてコンビニが「どういうところ」に出店しているのかを傾向として見ることができます。図示の力といえますね。具体的には、こんなことが読み取れるんじゃないでしょうか？
まとめ
今回は、ボロノイ分割を用いてコンビニの出店傾向を探ってみました。図示ってやはり強力な手段です。何かの「仮説」を得る場合に、頭の中で想像するだけではなくて、「図示」することは非常に強力な手段になります。
罪滅ぼしのために奇麗に作った図とか...（苦笑）
M + geom_polygon(aes(x = x, y = y, fill = factor(tile)),data = tiledf, alpha = .5) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25,data = deldir$delsgs) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25,data = deldir$dirsgs, linetype = 2) + geom_polygon(data = hullPoly_df, colour = 'red', fill = NA,alpha = I(1/10), size = 2) + geom_point(aes(x=lon, y=lat), data = df,size = 5.75, alpha = .75) + geom_text(aes(x=lon,y=lat,label = 1:nrow(df)),data=df, size = 3, colour = 'green') + scale_fill_discrete(guide = FALSE)
