Как найти многоугольник, ближайший к точке в R?

У меня есть кадр данных пространственных точек и кадр данных пространственных полигонов. Например, мои полигоны будут многоугольником для каждого блока в Манхэттене. И точки - это люди, которые разбросаны повсюду, иногда падая посреди улицы, которая не является частью многоугольника.

Я знаю, как проверить, содержится ли точка внутри многоугольника, но как я могу назначить точки их ближайшему многоугольнику?

## 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)

1 ответ

Вот ответ, который использует подход, основанный на том, что описывается mdsumner в этом отличном ответе через несколько лет назад.

Одна важная нота (добавлена ​​как EDIT от 2/8/2015): rgeos, которая здесь используется для вычисления расстояний, ожидает, что геометрии, на которых она будет проецироваться в планарных координатах. Для этих примерных данных это означает, что они должны быть сначала преобразованы в координаты UTM (или какой-либо другой планарной проекции). Если вы допустили ошибку, оставив данные в своих исходных координатах с длинной длиной, рассчитанные расстояния будут неправильными, так как они будут обрабатывать градусы широты и долготы как имеющие равную длину.

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)

licensed under cc by-sa 3.0 with attribution.