De la ayuda sp::over :
x = "SpatialPoints", y = "SpatialPolygons" returns a numeric
vector of length equal to the number of points; the number is
the index (number) of the polygon of ‘y’ in which a point
falls; NA denotes the point does not fall in a polygon; if a
point falls in multiple polygons, the last polygon is
recorded.
Entonces, si convierte su SpatialPolygonsDataFrame en SpatialPolygons , recupera un vector de índices y puede subjugar sus puntos en NA :
> over(pts,as(ply,"SpatialPolygons"))
[1] NA 1 1 NA 1 1 NA NA 1 1 1 NA NA 1 1 1 1 1 NA NA NA 1 NA 1 NA
[26] 1 1 1 NA NA NA NA NA 1 1 NA NA NA 1 1 1 NA 1 1 1 NA NA NA 1 1
[51] 1 NA NA NA 1 NA 1 NA 1 NA NA 1 NA 1 1 NA 1 1 NA 1 NA 1 1 1 1
[76] 1 1 1 1 1 NA NA NA 1 NA 1 NA NA NA NA 1 1 NA 1 NA NA 1 1 1 NA
> nrow(pts)
[1] 100
> pts = pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),]
> nrow(pts)
[1] 54
> head([email protected])
var1 var2
2 0.04001092 v
3 0.58108350 v
5 0.85682609 q
6 0.13683264 y
9 0.13968804 m
10 0.97144627 o
>
Para los que dudan, aquí hay evidencia de que la sobrecarga de conversión no es un problema:
Dos funciones: primero el método de Jeffrey Evans, luego mi original, luego mi conversión pirateada, luego una versión basada en gIntersects basada en la respuesta de Josh O'Brien:
evans <- function(pts,ply){
prid <- over(pts,ply)
ptid <- na.omit(prid)
pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]
return(pt.poly)
}
rowlings <- function(pts,ply){
return(pts[!is.na(over(pts,as(ply,"SpatialPolygons"))),])
}
rowlings2 <- function(pts,ply){
class(ply) <- "SpatialPolygons"
return(pts[!is.na(over(pts,ply)),])
}
obrien <- function(pts,ply){
pts[apply(gIntersects(columbus,pts,byid=TRUE),1,sum)==1,]
}
Ahora, para un ejemplo del mundo real, he dispersado algunos puntos aleatorios sobre el conjunto de datos columbus :
require(spdep)
example(columbus)
pts=data.frame(
x=runif(100,5,12),
y=runif(100,10,15),
z=sample(letters,100,TRUE))
coordinates(pts)=~x+y
Se ve bien
plot(columbus)
points(pts)
Comprueba que las funciones están haciendo lo mismo:
> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] TRUE
Y ejecute 500 veces para la evaluación comparativa:
> system.time({for(i in 1:500){evans(pts,columbus)}})
user system elapsed
7.661 0.600 8.474
> system.time({for(i in 1:500){rowlings(pts,columbus)}})
user system elapsed
6.528 0.284 6.933
> system.time({for(i in 1:500){rowlings2(pts,columbus)}})
user system elapsed
5.952 0.600 7.222
> system.time({for(i in 1:500){obrien(pts,columbus)}})
user system elapsed
4.752 0.004 4.781
Según mi intuición, no es una gran sobrecarga, de hecho podría ser una sobrecarga menor que la conversión de todos los índices de fila a caracteres y viceversa, o ejecutar na.omit para obtener valores perdidos. Lo que incidentalmente conduce a otro modo de falla de la función evans ...
Si una fila del marco de datos de los polígonos es NA (lo que es perfectamente válido), entonces la superposición con SpatialPolygonsDataFrame para los puntos en ese polígono producirá un marco de datos de salida con todos los NA s, que% entonces evans() caerá:
> [email protected][1,]=rep(NA,20)
> [email protected][5,]=rep(NA,20)
> [email protected][17,]=rep(NA,20)
> [email protected][15,]=rep(NA,20)
> set.seed(123)
> pts=data.frame(x=runif(100,5,12),y=runif(100,10,15),z=sample(letters,100,TRUE))
> coordinates(pts)=~x+y
> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] FALSE
> dim(evans(pts,columbus))
[1] 27 1
> dim(rowlings(pts,columbus))
[1] 28 1
>
PERO gIntersects es más rápido, incluso con la necesidad de barrer la matriz para verificar las intersecciones en R en lugar de en el código C. Sospecho que es el prepared geometry skills de GEOS, creando índices espaciales; sí, con prepared=FALSE lleva un poco más de tiempo, unos 5,5 segundos.
Me sorprende que no haya una función para devolver directamente los índices o los puntos. Cuando escribí splancs hace 20 años, las funciones de punto en polígono tenían ambas ...