I have a spatial points data frame and a spatial polygons data frame. For example, my polygons would be a polygon for each block in Manhattan. And the points are people, which are scattered all over, sometimes falling in the middle of a street, which is not part of a polygon.
I know how to check if a point is contained inside a polygon, but how could I assign points to their closest polygon?
## 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)
Here's an answer that uses an approach based on that described by mdsumner in this excellent answer from a few years back.
One important note (added as an EDIT on 2/8/2015): rgeos, which is here used to compute distances, expects that the geometries on which it operates will be projected in planar coordinates. For these example data, that means that they should be first transformed into UTM coordinates (or some other planar projection). If you make the mistake of leaving the data in their original lat-long coordinates, the computed distances will be incorrect, as they will have treated degrees of latitude and longitude as having equal lengths.
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"
crs <- CRS(sprintf(utmStr, 32))
pUTM <- spTransform(p, crs)
ptsUTM <- spTransform(pts, crs)
## Set up containers for results
n <- length(ptsUTM)
nearestCanton <- character(n)
distToNearestCanton <- numeric(n)
## For each point, find name of nearest polygon (in this case, Belgian cantons)
for (i in seq_along(nearestCanton)) {
gDists <- gDistance(ptsUTM[i,], pUTM, byid=TRUE)
nearestCanton[i] <- pUTM$NAME_2[which.min(gDists)]
distToNearestCanton[i] <- min(gDists)
}
## Check that it worked
data.frame(nearestCanton, distToNearestCanton)
# nearestCanton distToNearestCanton
# 1 Wiltz 15342.222
# 2 Echternach 7470.728
# 3 Remich 20520.800
# 4 Clervaux 6658.167
# 5 Echternach 22177.771
# 6 Clervaux 26388.388
# 7 Redange 8135.764
# 8 Remich 2199.394
# 9 Esch-sur-Alzette 11776.534
# 10 Remich 14998.204
plot(pts, pch=16, col="red")
text(pts, 1:10, pos=3)
plot(p, add=TRUE)
text(p, p$NAME_2, cex=0.7)