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

17

Tengo un SpatialPointsDataFrame con algunos datos adicionales. Me gustaría extraer esos puntos dentro de un polígono y, al mismo tiempo, conservar el objeto SPDF y sus datos correspondientes.

Hasta ahora he tenido poca suerte y he recurrido a la combinación y fusión a través de una ID común, pero esto solo funciona porque tengo datos de cuadrícula 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 utilizar 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
    
pregunta Roman Luštrik 18.06.2013 - 16:51

6 respuestas

20

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

    
respondido por el Spacedman 19.06.2013 - 11:51
11

sp proporciona una forma más corta para seleccionar entidades basadas en intersecciones espaciales, siguiendo el ejemplo de OP:

pts[ply,]

a partir de:

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

Detrás de las escenas esto es corto para

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

Lo que hay que tener en cuenta es que existe un método geometry que elimina los atributos: over cambia el comportamiento si su segundo argumento tiene atributos o no (esta fue la confusión de OP). Esto funciona en todas las clases de Spatial * en sp , aunque algunos métodos de over requieren rgeos , consulte this vignette para más detalles, por ejemplo, el caso de coincidencias múltiples para polígonos superpuestos.

    
respondido por el Edzer Pebesma 21.02.2015 - 10:25
6

Estabas en el camino correcto con más. Los rownames 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) 
    
respondido por el Jeffrey Evans 18.06.2013 - 19:33
4

¿Esto es lo que buscas?

Una nota, en la edición: la llamada envolvente a apply() es necesaria para que esto funcione con objetos arbitrarios SpatialPolygons , que posiblemente contengan más de una entidad poligonal. Gracias a @Spacedman por ayudarme 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)
    
respondido por el Josh O'Brien 19.06.2013 - 01:15
4

Aquí hay un posible enfoque utilizando el paquete rgeos . Básicamente, utiliza la función gIntersection que le permite intersectar dos objetos sp . Al extraer las ID de esos puntos que se encuentran dentro del polígono, luego podrá subordenar su SpatialPointsDataFrame original, manteniendo todos los datos correspondientes. El código es casi autoexplicativo, 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([email protected])[[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
    
respondido por el fdetsch 18.06.2013 - 17:53
1

Hay una solución extremadamente simple que utiliza la biblioteca spatialEco .

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
    
respondido por el rafa.pereira 27.10.2015 - 22:33

Lea otras preguntas en las etiquetas