自分のための「Rで重心ボロノイ基準？」

先週、谷岡氏と二人で相談して毎週「意味不明」な分析手法を開発することにした。
今週の課題は、「重心ボロノイ」を用いた方法の検討。
この手法は、ボロノイセルのサンプルを、そのセルの重心へと移動し続けるというもの。

逐次計算を行うので、どこで止めるかが問題になるが、
現段階では、重心の移動量がある閾値よりも小さくなったときに終了。
距離の計算方法も色々と設定できそうではあるが、今は考慮しない。

この方法、何に役立つのか？

ふむ。最終的に得られる重心ボロノイは、同じ条件下での領域だから、
つまりは、地理空間に適用した場合、最終的に得られた重心ボロノイは、
潜在的に存在する地理的要因を取り除いたものになるのか？？？

例えば、土地起伏、植生、河川などの水場、「不吉」な空間...etc
そういった地理的要因を取り除いた空間におけるボロノイってことか？
いや、よく解らない。まぁ、これで論文書くわけじゃないし...まぁ、いいや。

実用化するには、もう少し、色々と検討する余地があるが、
重心ボロノイの面積（の逆数）や、初期ボロノイのサンプル点からの移動量、
そういった値の分散を取るのも面白い。潜在的な要因を探ることができそう。

今度、何かの研究で使ってみるかな？

今回は、Rでコロプレスを描く方法も利用する。
長いソースコードをコピペすると、時々、エラーが出るのでスクリプトファイルに。
とりあえず、スクリプトファイルは以下の場所からダウンロード可能。

https://docs.google.com/open?id=0B9OtIs5BjRVhQXRaUk5ZOUEyUXM

とにかく、実行コード

---

# 最初にR スクリプトのソースを読み込む。
source("path/to/source/cis.geom.r")

# シェープファイルを読み込む
shp <- readShapePoints("path/to/shapefile/shapes.shp")

# 重心ボロノイを描く。
cv <- cis.centroidalVoronoi(shp)

# 結果をプロットする。
par(mfrow=c(1,2))

plot(deldir(shp@coords[,1], shp@coords[,2]),wlines="tess", pch=16, main="Original")

cis.choropleth(cv, "a", main="Centroidal Voronoi")

っで、Rスクリプに含まれているソースコード。

---

cis.centroidalVoronoi <- function(shp, eps=0.05, verbose=0){
# Load required packages.
require(voronoi)
require(maptools)
require(utils)

# Build original voronoi diagram.
obj <- deldir(shp@coords[,1], shp@coords[,2], verbose=0)

# Build centroidal voronoi.
c.obj <- centroidal(obj, eps=eps)

n <- length(tile.list(c.obj$T))
l <- list()
a <- vector(length=n)

# Create a text progressbar.
pb <- txtProgressBar(min=0, max=n)
for(i in 1:n){
 # Get coordinates of the tesseract.
 x <- c(tile.list(c.obj$T)[[i]]$x, tile.list(c.obj$T)[[i]]$x[1])
 y <- c(tile.list(c.obj$T)[[i]]$y, tile.list(c.obj$T)[[i]]$y[1])

 # Create a polygon from the set of coordinates.
 p <- Polygon(matrix(c(x,y),ncol=2), hole=FALSE)

 # Convert a polygon to the polygons.
 l[i] <- Polygons(list(p), i)
 a[i] <- l[[i]]@area^(-1)

 # Update the text progressbar.
 setTxtProgressBar(pb, i)
}

# Create a spatial polygon from the list of polygons.
sp <- SpatialPolygons(l)
sp <- addAttrToGeom(sp, data.frame(a), FALSE)

# Return centroidal voronoi with spatial polygons dataframe.
return(sp)
}