ここでは、rgeosトポロジライブラリでgDistance関数を使用しています。私はブルートフォースダブルループを使用していますが、それは驚くほど速いです。 142ポイントと10ポリゴンでは2秒以下です。私はループを実行するためにもっと優雅な方法があると確信しています。
require(rgeos)
# CREATE SOME DATA USING meuse DATASET
data(meuse)
coordinates(meuse) <- ~x+y
pts <- meuse[sample(1:dim(meuse)[1],142),]
data(meuse.grid)
coordinates(meuse.grid) = c("x", "y")
gridded(meuse.grid) <- TRUE
meuse.grid[["idist"]] = 1 - meuse.grid[["dist"]]
polys <- as(meuse.grid, "SpatialPolygonsDataFrame")
polys <- polys[sample(1:dim(polys)[1],10),]
plot(polys)
plot(pts,pch=19,cex=1.25,add=TRUE)
# LOOP USING gDistance, DISTANCES STORED IN LIST OBJECT
Fdist <- list()
for(i in 1:dim(pts)[1]) {
pDist <- vector()
for(j in 1:dim(polys)[1]) {
pDist <- append(pDist, gDistance(pts[i,],polys[j,]))
}
Fdist[[i]] <- pDist
}
# RETURN POLYGON (NUMBER) WITH THE SMALLEST DISTANCE FOR EACH POINT
(min.dist <- unlist(lapply(Fdist, FUN=function(x) which(x == min(x))[1])))
# RETURN DISTANCE TO NEAREST POLYGON
(PolyDist <- unlist(lapply(Fdist, FUN=function(x) min(x)[1])))
# CREATE POLYGON-ID AND MINIMUM DISTANCE COLUMNS IN POINT FEATURE CLASS
[email protected] <- data.frame([email protected], PolyID=min.dist, PDist=PolyDist)
# PLOT RESULTS
require(classInt)
(cuts <- classIntervals([email protected]$PDist, 10, style="quantile"))
plotclr <- colorRampPalette(c("cyan", "yellow", "red"))(20)
colcode <- findColours(cuts, plotclr)
plot(polys,col="black")
plot(pts, col=colcode, pch=19, add=TRUE)
min.distベクトルは、ポリゴンの行番号を表します。たとえば、このベクトルをそのまま使用して、最も近いポリゴンをサブセット化することができます。
near.polys <- polys[unique(min.dist),]
PolyDistベクトルには、フィーチャの投影単位の実際のデカルト距離が含まれています。
最後の段落で判断すると、数学的な問題があるように見えます。フォローの比較セットを作るよりも優れたアルゴリズムを見つけるのは正しいでしょうか?それは数学SEの方が適しているかもしれません。 – Frank
'spatstat'パッケージはあなたが望むことをすることができます。私はそのツールセットの専門家ではないので、確かに確認することはできません。 –
ここに含まれる数字では、10ポリゴンと142ポイント(1420距離!)のブルートフォースが問題になるはずはありません。 'plyr'パッケージは、あなたがループを気に入らない場合に役立ちます。 – Gregor