deldirでコンビニ位置データのボロノイ分割を行う。

この記事はRAdventCalendar2013に参加しています。
この記事は4日目になります。

これがはじめてというわけではないですがRのdeldirパッケージを使ってみました。
パッケージdeldirはドロネー（Delaunay）分割とディリクレ(Dirichlet)分割（ボロノイ分割、ティーセン分割）に関連する処理を含むパッケージです。

ボロノイ分割とドロネー分割の違いについてはこちら

パッケージdeldirの胆になるのが関数delirです。分割に当たってはLee-Schaterアルゴリズムを使用しているようです。

deldir関数の解説

deldir(x, y, dpl=NULL, rw=NULL, eps=1e-09, sort=TRUE, plotit=FALSE,
digits=6, z=NULL, zdum=NULL, suppressMsge=FALSE, ...)

引数

x,y:分割したい座標の入力点
dpl:リスト
rw:分割を行う境界領域を指定する。デフォルトだと入力点の端に対してプラス10％になる。
sort:TRUEがデフォルト。TRUEにセットするとデータが分割前にソートされ計算が効率的になる。
FALSEだと追加点が最後に付加されるので、デバッグ用に設定するための変数である。
plotit:plotを行うか否かのフラグである。
digists:数値の丸め桁数、デフォルトは6桁である。
z:各入力点に対応する重みの値である。
zdum：各ダミー入力点に対応する重みの値である。

返り値

delsgs:ドロネー三角形エッジリストを表すデータフレームである。エッジの座標点（x1,y1,x2,y2）が入力され、各ドロネー点のIDがそれに続く。
dirsgs:ディリクレ領域（ボロノイ領域）エッジリストを表すデータフレームである。エッジの座標点（x1,y1,x2,y2）が入力され、各ドロネー点のIDがそれに続く。
7番目の点はエッジの両端が境界領域と接しているかどうかを表す。
summary:分割されるデータ点とダミー点の情報を含むデータフレーム
-x,y(,z):点とその重みの位置情報
-n.tri:当該点におけるドロネー三角形の数
-del.area:当該点におけるドロネー三角形の面積の3分の1
-del.wts:del.areaを比率に変化させたもの。
-n.tside:ディリクレタイルの他のディリクレタイルとの接触数
-nbpt:ディリクレタイルの境界領域との接触数
-dir.area:ディリクレタイルの面積
-dir.wts:dir.areaの面積を比率に変化させたもの。
n.data:分割されたデータ点、重複は除く。
n.dum:分割されたデータ点、重複点は除く。
del.area:三角に分割された点の凸包の面積、del.areaの合計と一致。
dir.area:境界領域によって囲まれている面積、dir.areaの合計と一致。
rw:正方形領域の枠。

パッケージdeldirその他の関数

tile.list():deldirを受け取って、ボロノイ分割されたタイルを返す。
plot.deldir():deldirを受け取って、プロットする。

以下では、前回の記事でも用いたコンビニデータから東京のコンビニ6009店*1をボロノイ分割してポリゴンをシェープファイルに出力するコードを記載します。

コード

library(deldir)
library(maptools)


#データは事前に準備している。
#> head(df)
#storename        storetype        x        y
#1 セブンイレブン  世田谷中央病院店 セブンイレブン 139.6511 35.64166
#2 ファミリーマート  池袋二丁目店 ファミリーマート 139.7090 35.73495
#3 みんなのイチバ  足立江北２丁目店 その他 139.7637 35.76921


#ボロノイ分割を行う。rwオプションを指定し外枠を拡張させない。
dd<-deldir(df$x, df$y, rw=c(min(df$x),max(df$x),min(df$y),max(df$y)))
#タイルを取得する。
w<-tile.list(dd)
#ポリゴンベクターの初期化
polys<-vector(mode='list', length=length(w))
for (i in seq(along=polys)) {
  #ボロノイ点を取得している
  pcrds<-cbind(w[[i]]$x, w[[i]]$y)
  #閉路化してポリゴン座標に変更する。
  pcrds<-rbind(pcrds, pcrds[1,])
  #ポリゴンを作成する。
  polys[[i]]<-Polygons(list(Polygon(pcrds)), ID=as.character(i))
}

#スパシャルポリゴンオブジェクトの作成
SP<-SpatialPolygons(polys)
#ボロノイ分割したデータをSpatialDataFrameクラスに変更する。
voronoi<-SpatialPolygonsDataFrame(SP,data=data.frame(x=df$x,y=df$y,storename=df$storename, storetype=df$storetype,
  row.names=sapply(slot(SP, 'polygons'),function(x) slot(x, 'ID'))))
#shpファイルにして出力する。
writeSpatialShape(voronoi, fn="vornoi-combini")

結果
以下は結果のシェープファイルをQGISで表示させたものです（赤がセブンイレブン、青がローソン、緑がファミリーマート、オレンジがサークルK・サンクス、その他が紫）。

広域図