如何找到最接近R中某点的多边形？

Question

我有一个空间点数据框和一个空间多边形数据框。例如，我的多边形将是曼哈顿中每个块的多边形。分数是人，分散在各处，有时落在街道中间，不是多边形的一部分。

我知道如何检查多边形中是否包含一个点，但是如何将点分配给它们最接近的多边形？

## Make some example  data
set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
p <- shapefile(system.file("external/lux.shp", package="raster"))
p2 <- as(1.5*extent(p), "SpatialPolygons")
proj4string(p2) <- proj4string(p)
pts <- spsample(p2-p, n=10, type="random")

## Plot to visualize
plot(pts, pch=16, cex=.5,col="red")
plot(p, col=colorRampPalette(blues9)(12), add=TRUE)

Answer 1

这是一个答案，使用的方法基于几年前this excellent answer中mdsumner所描述的方法。

一个重要的注释（在2015年2月8日作为编辑添加）：rgeos，这里用于计算距离，期望它操作的几何将以平面坐标投影。对于这些示例数据，这意味着它们应首先转换为UTM坐标（或一些其他平面投影）。如果您错误地将数据保留在原始的lat-long坐标中，则计算的距离将是不正确的，因为它们将纬度和经度的处理程度视为具有相等的长度。

library(rgeos)

##  First project data into a planar coordinate system (here UTM zone 32)
utmStr <- "+proj=utm +zone=%d +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
crs <- CRS(sprintf(utmStr, 32))
pUTM <- spTransform(p, crs)
ptsUTM <- spTransform(pts, crs)

## Set up container for results
n <- length(ptsUTM)
nearestCantons <- character(n)

## For each point, find name of nearest polygon (in this case, Belgian cantons)
for (i in seq_along(nearestCantons)) {
    nearestCantons[i] <- pUTM$NAME_2[which.min(gDistance(ptsUTM[i,], pUTM, byid=TRUE))]
}

## Check that it worked
nearestCantons
# [1] "Wiltz"            "Echternach"       "Remich"           "Clervaux"        
# [5] "Echternach"       "Clervaux"         "Redange"          "Remich"          
# [9] "Esch-sur-Alzette" "Remich"   

plot(pts, pch=16, col="red")
text(pts, 1:10, pos=3)
plot(p, add=TRUE)
text(p, p$NAME_2, cex=0.7)

Answer 2

我在这里参加派对的时间已经很晚了，但我刚刚找到了这个帖子，并且为了它的价值，提供了这个建议。 RANN软件包的nn2功能允许您仅在一些有限半径内搜索（最近点），这可以节省大量时间。我的建议是在多边形上添加点，将点与多边形相关联，然后搜索最近的点。看起来gDistance方法在点数不多时更快，但是nn2方法可以更好地扩展到更大的问题，因为它搜索有限的半径（当然，如果半径内没有点，它将无法找到匹配，所以必须正确选择半径）。我是新手，所以这可能不是最佳的。如果gDistance允许限制搜索会很好。

## Make some example  data
set.seed(1)
library(raster)
library(rgdal)
library(rgeos)
library(RANN)
library(spatialEco)

p <- shapefile(system.file("external/lux.shp", package="raster"))

##  Project data into a planar coordinate system (here UTM zone 32)
utmStr <- "+proj=utm +zone=%d +datum=NAD83 +units=m +no_defs +ellps=GRS80 +towgs84=0,0,0"
crs <- CRS(sprintf(utmStr, 32))
pUTM <- spTransform(p, crs)

# the points of interest (all within some threshold distance of the polygons)
ptsUTM <- spsample(gBuffer(pUTM,width=2000)-pUTM, n=10000, type="random")
## Plot to visualize
plot(ptsUTM, pch=16, cex=.5,col="red")
plot(pUTM, col=colorRampPalette(blues9)(12), add=TRUE)

# the gDistance method
starttime <- Sys.time()
## Set up container for results
n <- length(ptsUTM)
nearestCantons <- character(n)
## For each point, find name of nearest polygon (in this case, Belgian cantons)
for (i in seq_along(nearestCantons)) {
    nearestCantons[i] <- pUTM$NAME_2[which.min(gDistance(ptsUTM[i,], pUTM, byid=TRUE))]
}
Sys.time()-starttime

# the nn2 method
starttime <- Sys.time()
## create search points and associate with polygon attributes
rp <- spsample(pUTM,n=10000,type="regular")
rp2 <- point.in.poly(rp,pUTM)
# search for nearest point (with radius)
nn <- nn2(coordinates(rp2),coordinates(ptsUTM),k=1,searchtype="radius",radius=5000)$nn.idx
nearestCantons2 <- rp2$NAME_2[nn]
Sys.time()-starttime