Analyze IT.

2013-12-04

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 みんなのイチバ  足立江北2丁目店 その他 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サンクス、その他が紫)。


広域図

f:id:EulerDijkstra:20131202180331j:image


詳細図

f:id:EulerDijkstra:20131202180332j:image


2013年も一月を切りました。皆様よいクリスマス、年越しを。

*1:データ数は6182だったがジオコードの結果、同一座標となった店舗からはランダムに選出した。

スパム対策のためのダミーです。もし見えても何も入力しないでください
ゲスト


画像認証

トラックバック - http://d.hatena.ne.jp/EulerDijkstra/20131204/1386118237
リンク元