¿Cómo superponer un polígono sobre SpatialPointsDataFrame y preservar los datos SPDF?

Tengo un SpatialPointsDataFramecon algunos datos adicionales. Me gustaría extraer esos puntos dentro de un polígono y, al mismo tiempo, preservar el SPDFobjeto y sus datos correspondientes.

Hasta ahora he tenido poca suerte y he recurrido a la coincidencia y la fusión a través de una ID común, pero esto solo funciona porque tengo datos cuadriculados con IDS individuales.

Aquí hay un ejemplo rápido, estoy buscando puntos dentro del cuadrado rojo.

library(sp)
set.seed(357)
pts <- data.frame(x = rnorm(100), y = rnorm(100), var1 = runif(100), var2 = sample(letters, 100, replace = TRUE))
coordinates(pts) <- ~ x + y
class(pts)
plot(pts)
axis(1); axis(2)

ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol = 2, byrow = TRUE)
ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID = 1)))
ply <- SpatialPolygonsDataFrame(Sr = ply, data = data.frame(polyvar = 357))
plot(ply, add = TRUE, border = "red")

El enfoque más obvio sería usar over, pero esto devuelve los datos del polígono.

> over(pts, ply)
    polyvar
1        NA
2       357
3       357
4        NA
5       357
6       357

De la sp::overayuda:

 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 SpatialPolygonsDataFramea SpatialPolygons, obtiene un vector de índices y puede subconjugar 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(pts@data)
         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 gIntersectsla 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 columbusconjunto de datos:

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)

Verifique que las funciones estén haciendo lo mismo:

> identical(evans(pts,columbus),rowlings(pts,columbus))
[1] TRUE

Y corra 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 menos una sobrecarga que convertir todos los índices de fila en caracteres y viceversa, o ejecutar na.omit para obtener valores perdidos. Lo que incidentalmente conduce a otro modo de falla de la evansfunción ...

Si una fila del marco de datos de polígonos es todo NA(lo cual es perfectamente válido), entonces la superposición con SpatialPolygonsDataFramepuntos en ese polígono producirá un marco de datos de salida con todos los NAs, que evans()luego caerán:

> columbus@data[1,]=rep(NA,20)
> columbus@data[5,]=rep(NA,20)
> columbus@data[17,]=rep(NA,20)
> columbus@data[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 gIntersectses más rápido, incluso teniendo que barrer la matriz para verificar las intersecciones en R en lugar de en el código C. Sospecho que son las prepared geometryhabilidades de GEOS, la creación de índices espaciales, sí, con prepared=FALSEesto 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í hace splancs20 años, las funciones de punto en el polígono tenían ambas ...

sp proporciona una forma más corta para seleccionar entidades basadas en la intersección espacial, siguiendo el ejemplo de OP:

pts[ply,]

a partir de:

points(pts[ply,], col = 'red')

Detrás de escena esto es la abreviatura de

pts[!is.na(over(pts, geometry(ply))),]

Lo que hay que tener en cuenta es que hay un geometrymétodo que elimina atributos: overcambia el comportamiento si su segundo argumento tiene atributos o no (esto fue la confusión de OP). Esto funciona en todas las clases de Spatial * sp, aunque algunos overmétodos lo requieren rgeos, vea esta viñeta para obtener detalles, por ejemplo, el caso de coincidencias múltiples para polígonos superpuestos.

Estabas en el camino correcto con más. Los nombres de fila del objeto devuelto corresponden al índice de fila de los puntos. Puede implementar su enfoque exacto con solo unas pocas líneas de código adicionales.

library(sp)
set.seed(357)

pts <- data.frame(x=rnorm(100), y=rnorm(100), var1=runif(100), 
                  var2=sample(letters, 100, replace=TRUE))
  coordinates(pts) <- ~ x + y

ply <- matrix(c(-1,-1, 1,-1, 1,1, -1,1, -1,-1), ncol=2, byrow=TRUE)
  ply <- SpatialPolygons(list(Polygons(list(Polygon(ply)), ID=1)))
    ply <- SpatialPolygonsDataFrame(Sr=ply, data=data.frame(polyvar=357))

# Subset points intersecting polygon
prid <- over(pts,ply)
  ptid <- na.omit(prid) 
    pt.poly <- pts[as.numeric(as.character(row.names(ptid))),]  

plot(pts)
  axis(1); axis(2)
    plot(ply, add=TRUE, border="red")
      plot(pt.poly,pch=19,add=TRUE) 

¿Es esto lo que buscas?

Una nota, en edición: la llamada de ajuste a apply()es necesaria para que esto funcione con SpatialPolygonsobjetos arbitrarios , que posiblemente contengan más de una función de polígono. Gracias a @Spacedman por animarme a demostrar cómo aplicar esto al caso más general.

library(rgeos)
pp <- pts[apply(gIntersects(pts, ply, byid=TRUE), 2, any),]


## Confirm that it works
pp[1:5,]
#              coordinates       var1 var2
# 2 (-0.583205, -0.877737) 0.04001092    v
# 3   (0.394747, 0.702048) 0.58108350    v
# 5    (0.7668, -0.946504) 0.85682609    q
# 6    (0.31746, 0.641628) 0.13683264    y
# 9   (-0.469015, 0.44135) 0.13968804    m

plot(pts)
plot(ply, border="red", add=TRUE)
plot(pp, col="red", add=TRUE)

Aquí hay un posible enfoque usando el rgeospaquete. Básicamente, hace uso de la gIntersectionfunción que le permite intersecar dos spobjetos. Al extraer los ID de esos puntos que se encuentran dentro del polígono, puede subconjugar su original SpatialPointsDataFramey conservar todos los datos correspondientes. El código casi se explica por sí mismo, pero si tiene alguna pregunta, ¡no dude en preguntar!

# Required package
library(rgeos)

# Intersect polygons and points, keeping point IDs
pts.intersect <- gIntersection(ply, pts, byid = TRUE)

# Extract point IDs from intersected data
pts.intersect.strsplit <- strsplit(dimnames(pts.intersect@coords)[[1]], " ")
pts.intersect.id <- as.numeric(sapply(pts.intersect.strsplit, "[[", 2))

# Subset original SpatialPointsDataFrame by extracted point IDs
pts.extract <- pts[pts.intersect.id, ]

head(coordinates(pts.extract))
              x          y
[1,] -0.5832050 -0.8777367
[2,]  0.3947471  0.7020481
[3,]  0.7667997 -0.9465043
[4,]  0.3174604  0.6416281
[5,] -0.4690151  0.4413502
[6,]  0.4765213  0.6068021

head(pts.extract)
         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

Hay una solución extremadamente simple usando la spatialEcobiblioteca.

library(spatialEco)

# intersect points in polygon
  pts <- point.in.poly(pts, ply)

# check plot
  plot(ply)
  plot(a, add=T)

# convert to data frame, keeping your data
  pts<- as.data.frame(pts)

Comprueba el resultado:

pts

>             x          y       var1 var2 polyvar
> 2  -0.5832050 -0.8777367 0.04001092    v     357
> 3   0.3947471  0.7020481 0.58108350    v     357
> 5   0.7667997 -0.9465043 0.85682609    q     357
> 6   0.3174604  0.6416281 0.13683264    y     357
> 9  -0.4690151  0.4413502 0.13968804    m     357
> 10  0.4765213  0.6068021 0.97144627    o     357