未来はキミドリイロ

サイコメトリックアイドルを目指す心理学徒の勉強部屋です。勉強・趣味などについて書いています。

【R】ボロノイ分割で本当の最寄り駅が知りたい

最近知ったボロノイ分割というのがいろいろと使い道がありそうなので触ってみたくて,ちょちょっとRでやってみました。

ボロノイ分割とは?

私の理解したところでは,ボロノイ図(Voronoi diagram)とは,平面上に設定された「母点」と呼ばれるいくつかの座標をもとに,平面上の座標空間をどの母点に最も近いかによって分割することで作成できる図のことだと言えます。
また,この分割のことをボロノイ分割(英語表記は脚注参照*1)といいます。

せっかくなので備忘もかねて,自分流にかみ砕いてつくり方の概念的なものを図示してみます。
頭を使うのが面倒だったのでペイントで作りました,すいません。

座標平面上に母点を置く

Figure1
f:id:arca821:20160515072042p:plain

黒線で囲まれた内側が(2次元の)座標空間であるとして,この画像のオレンジ色の点のように,母点を置きます。
このとき,母点の数が2つ以下ですと次に支障を来すので,必ず3つ以上の母点を置くようにするのがいいかもしれません*2

ドロネー図の作成

次に,Figure1をもとに,ドロネー図(Delaunay diagram)と呼ばれるものを作成します。
つくり方は簡単で,母点同士を線で結ぶことで得られます。

Figure2
f:id:arca821:20160515073100p:plain

この,母点同士を結ぶ線のことをドロネー辺,ドロネー辺で囲まれた三角の領域のことをドロネー三角形と呼びます。

ドロネー辺の垂直二等分線を考える

Figure2で作成したドロネー図について,今度は全てのドロネー辺における垂直二等分線を考えます。
図では細い線を用いて表現しています。
ペイントを用いて目方で作っているので,厳密に垂直二等分できていないところもあるかもしれませんが目をつぶっていただけると助かります。
全ての辺に対して線を引いている都合でごちゃごちゃしていますが,各線がどのドロネー辺の垂直二等分線になっているのかをちらりと確認しながら見るのがよいかと思います。

Figure3
f:id:arca821:20160515074258p:plain

ボロノイ図を作成する

Figure3はかなりごちゃごちゃした図になっていますが,4つの三角形について個々に見ていくとわかりやすくなります。
各々の三角形は3つの(ドロネー)辺から構成され,それらの辺の垂直二等分線が交わる点は,中学数学でも習う「三角形の重心」となっています。
まず,各三角形の重心を赤丸で書いてみました。
ひとつだけ画像の上端を突き破って外に出てしまっていますが,ちゃんと4つあります。

Figure4
f:id:arca821:20160515080308p:plain

たくさん引いた垂直二等分線の中から,この重心同士を結ぶ線分を考えます。

Figure5
f:id:arca821:20160515080724p:plain

このFigure5に,垂直二等分線のうち重心と座標平面の端とを結ぶ線分になっているものにのみ赤入れをしてあげれば,ボロノイ図の完成です。

Figure6
f:id:arca821:20160515080906p:plain

赤い線(と座標空間の端)で囲まれた部分はボロノイ領域と呼びます。
ボロノイ領域を区分けする線をボロノイ境界,今回は重心となっているボロノイ境界の交点をボロノイ点と呼びます。
見やすいように,ボロノイ領域ごとに色分けしてみました。

Figure7
f:id:arca821:20160515081232p:plain

この図は,例えば水色の部分に含まれる全ての点は,5つある母点のうち水色の領域内にある母点が最も近いものであるということを表します。

本当の最寄り駅を知りたい!

ここまでを踏まえて,ありきたりだとは思いますが「あれ,これ学区とか最寄り駅とか考えるのに使えるんじゃないかなー」と思ったので試しにRで実行してみようというのがこの記事の主旨です。

そこで今回は東京大学本郷キャンパスをピックアップしてみました。
本郷に足を運んでみた方はわかると思うのですが,非常に広大なキャンパスで,しかもアクセスマップを見てみると4つの駅に囲まれています*3
f:id:arca821:20160515082807g:plain
ピックアップした理由はそこでして,最寄り駅っぽいものが4つあるけれど,結局どの駅を使うのが便利なんだろうというのをRを使って見てみようと思っています。

データの準備

緯度経度のデータを使いたいので,駅の緯度経度を調べます。
私は次のサービスを使いました。
http://www.odekakemap.com/station/

一応,今回読み込んだデータはこんな感じになっています,というのを共有しておきます。

station y x
Hongo3chome 35.706671 139.759914
Nezu 35.717373 139.765741
Todaimae 35.717633 139.758025
Yushima 35.708243 139.769916

yが緯度で,xが経度です。

パッケージの準備

今回は5つのパッケージを利用しました。

  • ggplot2:言わずと知れた,図表作成の必須ツール
  • reshape2:名前から想像できるように,データ変形ツール
  • plyr:前処理で便利に使えるパッケージみたいです
  • deldir:今回の目玉,ボロノイ分割を行うためのパッケージ
  • ggmap:Google Mapを利用するためのパッケージ

Rで実行

いろいろ参考にして,こんなかんじでスクリプトを書きました。
文字インデントが下がっているところは,一行が長かったので改行した部分です。

#--------------------------------------------------------#
#下準備
setwd("ワークスペースのディレクトリ")

library(ggplot2)
library(plyr)
library(deldir)
library(reshape2)
library(ggmap)
#--------------------------------------------------------#
#データの読み込み
d<-read.csv("ut.csv",header=TRUE,sep=",")
#--------------------------------------------------------#
#google mapのR上での描画テスト
#読みだデータのうち,緯度・経度の最小~最大の範囲で位置情報を使う
loca<-c(min(d$x),min(d$y),max(d$x),max(d$y))
#zoomは見たい範囲が過不足なく映るように調整する
gmap<-ggmap(get_map(location=loca,zoom=15,source="google"))+xlab("")+ylab("")
#--------------------------------------------------------#
#voronoi分割
res1<-deldir(d$x,d$y)
#linetype=1,実線で分割
gmap+geom_point(aes(x=d$x,y=d$y,colour=factor(d$station)),size=3,data=d)
  +geom_segment(aes(x=x1,y=y1,xend=x2,yend=y2),size=.25,data=res1$dirsgs,linetype=1)
#--------------------------------------------------------#
#ボロノイ領域の着色用にタイルをうんぬんかんぬん
tilelist<-unclass(tile.list(res1))
tilelist2<-lapply(tilelist,function(l){
data.frame(x=l$x,y=l$y)
})

tiledf<-melt(tilelist2,id.vars=c("x"))
tiledf2<-tiledf[,c(4,1,3)]
names(tiledf2)<-c("tile","x","y")
#--------------------------------------------------------#
#polygonを使って着色
gmap+geom_polygon(aes(x=x,y=y,fill=factor(tile)),data=tiledf2,colour="Black",alpha=.3)
  +geom_point(aes(x=x,y=y,colour=factor(station)),size=2,data=d)
  +theme_set(theme_bw(base_size=12,base_family="HiraMaruProN-W4"))

途中でフォントがないですみたいなエラーが出ると思いますが無視で大丈夫です。
こんな結果がアウトプットされます。
f:id:arca821:20160515094556p:plain

これにキャンパスマップを組み合わせれば,例えば東大病院に行きたければ○○駅,三四郎池を見に行きたければ○○駅,という意思決定が簡単にできそうですよね。

使えそうなこと

家探ししているときに最寄りがやけにいっぱいあるけど,ぶっちゃけどれでしょうなんてときにはこういう視点もあるかもしれないなーというかんじでしょうか。
今は母点の数が少ない=ボロノイ領域の数が少ないのであまり恩恵を受けることができませんが,最短ルート探しなんかにも使えそうですよね。
後は,今回は垂直二等分線を使いましたが何も2等分である必要はなく,重み付けもできると思うので池袋なんかの繁華街からの交通費を使って重みづけしてみるなんてこともできそうです。

Wikipedia曰く,画像の圧縮や離散データの集約にも使えるとのことで,なるほどと思いました。

今回はこんなかんじで,それではごきげんよう。

*1:英語版Wikipediaによると,a Voronoi tessellation, a Voronoi decomposition, a Voronoi partition, or a Dirichlet tessellationなどの呼び名があるとのこと

*2:おそらく取りうる点の限界・端を設定すれば2つからでもいけるような気がするのですが,はっきりしたことが言えないので,ここでは説明のため必ず3つ以上の母点をと書いています。

*3:画像は東大アクセスマップ(http://www.u-tokyo.ac.jp/campusmap/map01_02_j.html )よりお借りしました。