Come trovo il poligono più vicino a un punto in R?

Question

Ho un frame di dati punti spaziali e un frame di dati poligoni spaziali. Ad esempio, i miei poligoni sarebbero un poligono per ogni blocco di Manhattan. E i punti sono persone, che sono sparse dappertutto, a volte cadendo nel mezzo di una strada, che non fa parte di un poligono.

So come controllare se un punto è contenuto in un poligono, ma come posso assegnare i punti al poligono più vicino?

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

Ecco una risposta che utilizza un approccio basato su quello descritto da mdsumner in this excellent answer da qualche anno fa.

Una nota importante (aggiunto come EDIT 2015/02/08): rgeos, che viene qui usato per calcolare distanze, prevede che le geometrie su cui opera saranno proiettate in coordinate planari. Per questi dati di esempio, ciò significa che dovrebbero essere prima trasformati in coordinate UTM (o qualche altra proiezione planare). Se si commette l'errore di lasciare i dati nelle loro coordinate lat-long originali, le distanze calcolate saranno errate, in quanto avranno trattato gradi di latitudine e longitudine come aventi uguale lunghezza.

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)