comment superposer un polygone sur SpatialPointsDataFrame et préserver les données SPDF?

J'ai un SpatialPointsDataFrameavec quelques données supplémentaires. Je voudrais extraire ces points à l'intérieur d'un polygone et en même temps, conserver l' SPDFobjet et ses données correspondantes.

Jusqu'à présent, j'ai eu peu de chance et j'ai eu recours à la correspondance et à la fusion via un ID commun, mais cela ne fonctionne que parce que j'ai quadrillé les données avec des IDS individuels.

Voici un exemple rapide, je recherche des points à l'intérieur du carré rouge.

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")

L'approche la plus évidente serait d'utiliser over, mais cela renvoie les données du polygone.

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

De l' sp::overaide:

 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.

Donc, si vous convertissez votre SpatialPolygonsDataFrameen SpatialPolygonsvous récupérez un vecteur d'index et vous pouvez sous-définir vos points sur 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
> 

Pour les sceptiques, voici la preuve que la surcharge de conversion n'est pas un problème:

Deux fonctions - d'abord la méthode de Jeffrey Evans, puis ma version originale, puis ma conversion piratée, puis une version basée sur gIntersectsla réponse 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,]
}

Maintenant, pour un exemple réel, j'ai dispersé quelques points aléatoires sur l' columbusensemble de données:

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

Cela semble bon

plot(columbus)
points(pts)

Vérifiez que les fonctions font la même chose:

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

Et exécutez 500 fois pour l'analyse comparative:

> 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 

Selon mon intuition, ce n'est pas une surcharge importante, en fait, cela pourrait être moins d'une surcharge que la conversion de tous les index de ligne en caractère et retour, ou l'exécution de na.omit pour obtenir les valeurs manquantes. Ce qui conduit d'ailleurs à un autre mode de défaillance de la evansfonction ...

Si une ligne du bloc de données de polygones est tout NA(ce qui est parfaitement valide), la superposition avec SpatialPolygonsDataFramepour les points dans ce polygone produira un bloc de données de sortie avec tous les NAs, qui evans()tombera alors:

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

MAIS gIntersectsest plus rapide, même avec le balayage de la matrice pour vérifier les intersections dans R plutôt que dans le code C. Je soupçonne que ce sont les prepared geometrycompétences de GEOS, créant des index spatiaux - oui, avec prepared=FALSEcela prend un peu plus de temps, environ 5,5 secondes.

Je suis surpris qu'il n'y ait pas de fonction pour renvoyer directement les indices ou les points. Quand j'ai écrit splancsil y a 20 ans, les fonctions de point dans le polygone avaient les deux ...

sp fournit une forme plus courte pour sélectionner les entités en fonction de l'intersection spatiale, en suivant l'exemple OP:

pts[ply,]

au:

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

Dans les coulisses, c'est court pour

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

La chose à noter est qu'il existe une geometryméthode qui supprime les attributs: overchange le comportement si son deuxième argument a des attributs ou non (c'était la confusion de OP). Cela fonctionne dans toutes les classes Spatial * dans sp, bien que certaines overméthodes nécessitent rgeos, voir cette vignette pour plus de détails, par exemple le cas de plusieurs correspondances pour des polygones qui se chevauchent.

Vous étiez sur la bonne voie avec fini. Les noms renommés de l'objet retourné correspondent à l'index de ligne des points. Vous pouvez implémenter votre approche exacte avec seulement quelques lignes de code supplémentaires.

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) 

C'est ce que tu cherches?

Une remarque, lors de la modification: l'appel d'habillage à apply()est nécessaire pour que cela fonctionne avec des SpatialPolygonsobjets arbitraires , contenant éventuellement plusieurs entités surfaciques. Merci à @Spacedman de m'avoir incité à montrer comment l'appliquer au cas plus général.

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)

Voici une approche possible en utilisant le rgeospackage. Fondamentalement, il utilise la gIntersectionfonction qui vous permet d'intersecter deux spobjets. En extrayant les ID de ces points qui se trouvent dans le polygone, vous pouvez ensuite sous-définir votre original SpatialPointsDataFrame, en conservant toutes les données correspondantes. Le code est presque explicite, mais s'il y a des questions, n'hésitez pas à demander!

# 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

Il existe une solution extrêmement simple utilisant la spatialEcobibliothèque.

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)

Vérifiez le résultat:

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