這裏我使用rgeos拓撲庫中的gDistance函數。我正在使用一個強力雙循環,但它是驚人的快。它需要不到2秒,142點和10個多邊形。我相信有一種更優雅的方式來執行循環。
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 2013-05-08 19:34:39
'spatstat'包可能可以做你想做的。我不是那個工具集的專家,所以無法確定。 – 2013-05-08 19:39:48
在這裏涉及的數字中,10個多邊形和142個點(1420個距離!)的蠻力應該不成問題。如果你不喜歡for循環,'plyr'軟件包會幫助你。 – Gregor 2013-05-08 19:43:18