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

17

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
Roman Luštrik
la source
1
Merci d'avoir fourni un exemple reproductible. Aide toujours à essayer de comprendre un problème!
fdetsch

Réponses:

21

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

Spacedman
la source
Génial, cela fonctionne également pour plusieurs polygones (j'ai ajouté un exemple pour jouer avec dans la réponse de Joshua).
Roman Luštrik
Avec de grands ensembles de données de polygones, la coercition dans un objet SpatialPolygons représente beaucoup de surcharge et n'est pas nécessaire. L'application "sur" à un SpatialPolygonsDataFrame renvoie l'index de ligne qui peut être utilisé pour sous-définir les points. Voir mon exemple ci-dessous.
Jeffrey Evans
Un grand nombre des frais généraux? Il s'agit essentiellement de prendre l'emplacement @polygons de l'objet SpatialPolygonsDataFrame. Vous pouvez même le «simuler» en réaffectant la classe d'un SpatialPolygonsDataFrame à «SpatialPolygons» (bien que ce soit hacky et déconseillé). Tout ce qui va utiliser la géométrie devra de toute façon avoir cet emplacement à un moment donné, donc relativement parlant, il n'y a pas de surcharge du tout. C'est insignifiant de toute façon dans n'importe quelle application du monde réel où vous allez ensuite faire une charge de tests de polygones de points.
Spacedman
Il y a plus que des considérations de vitesse dans la prise en compte des frais généraux. En créant un nouvel objet dans l'espace de noms R, vous utilisez la RAM nécessaire. Lorsque ce n'est pas un problème dans un petit ensemble de données, cela affectera les performances avec de grandes données. R présente une disparition linéaire des performances. À mesure que les données augmentent, les performances prennent un coup. Si vous n'avez pas besoin de créer un objet supplémentaire, pourquoi le feriez-vous?
Jeffrey Evans
1
Nous ne le savions pas jusqu'à ce que je le teste tout à l'heure.
Spacedman
13

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.

Edzer Pebesma
la source
Bon à savoir! Je n'étais pas au courant de la méthode de la géométrie.
Jeffrey Evans
2
Bienvenue sur notre site, Edzer - c'est un plaisir de vous voir ici!
whuber
1
Merci Bill - ça devient plus silencieux sur stat.ethz.ch/pipermail/r-sig-geo , ou peut-être devrions-nous développer un logiciel qui cause plus de problèmes! ;-)
Edzer Pebesma
6

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) 
Jeffrey Evans
la source
Faux - les noms de domaine de l'objet retourné correspondent à l'index de ligne in_this_case - en général, les noms de ligne semblent être les noms de ligne des points - qui pourraient même ne pas être numériques. Vous pouvez modifier votre solution pour faire une correspondance de caractères qui pourrait la rendre un peu plus robuste.
Spacedman
@Sapcedman, ne sois pas si dogmatique. La solution n'est pas incorrecte. Si vous souhaitez sous-définir des points sur un ensemble de polygones ou affecter des valeurs de polygone à des points, la fonction Over fonctionne sans contrainte. Il y en a plusieurs était d'effectuer le passage pour piétons une fois que vous avez l'objet résultant résultant. Votre solution de contrainte sur un objet SpatialPolygon crée une surcharge considérable car cette opération peut être effectuée directement sur l'objet SpatialPolygonDataFrame. Au fait, avant de modifier un article, assurez-vous d'avoir raison. Le terme bibliothèque et package sont utilisés de manière interchangeable dans R.
Jeffrey Evans
J'ai ajouté quelques repères à mon message et repéré un autre problème avec votre fonction. Aussi "Les packages sont des collections de fonctions R, de données et de code compilé dans un format bien défini. Le répertoire où les packages sont stockés s'appelle la bibliothèque"
Spacedman
Bien que vous soyez techniquement correct en ce qui concerne "package" et "bibliothèque", vous discutez de sémantique. Je viens de recevoir une demande d'un éditeur de modélisation écologique pour que nous modifions notre utilisation de "package" (qui est en fait ma préférence) en "bibliothèque". Mon point est qu'ils deviennent des termes interchangeables et une question de préférence.
Jeffrey Evans
1
«techniquement correct», comme l'a fait remarquer le Dr Sheldon Cooper, «est la meilleure sorte de correct». Cet éditeur est techniquement mauvais, ce qui est le pire type de mal.
Spacedman
4

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)
Josh O'Brien
la source
Échoue horriblement si plypossède plusieurs fonctionnalités, car gIntersectsrenvoie une matrice avec une ligne pour chaque fonctionnalité. Vous pouvez probablement balayer les lignes pour une valeur VRAIE.
Spacedman
@Spacedman - Bingo. Besoin de faire apply(gIntersects(pts, ply, byid=TRUE), 2, any). En fait, je vais aller de l'avant et changer la réponse à cela, car elle englobe également le cas d'un seul polygone.
Josh O'Brien
Ah any. Cela pourrait être légèrement plus rapide que la version que je viens de comparer.
Spacedman
@Spacedman - D'après mes tests rapides, il ressemble obrienet rowlings2fonctionne au coude à coude, avec obrien peut-être 2% plus rapide.
Josh O'Brien
@ JoshO'Brien comment peut-on utiliser cette réponse sur de nombreux polygones? Cela ppdevrait avoir un IDqui indique dans quel polygone les points sont situés.
code123
4

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
fdetsch
la source
1
Devrait l' tmpêtre pts.intersect? En outre, l'analyse des noms de domaine renvoyés comme cela dépend du comportement non documenté.
Spacedman
@Spacedman, vous avez raison, vous avez tmpoublié de le supprimer lorsque vous avez terminé le code. En outre, vous avez raison d'analyser le dimnames. C'était une sorte de solution rapide pour fournir à l'interrogateur une réponse rapide, et il existe sûrement de meilleures approches (et plus universelles), par exemple la vôtre :-)
fdetsch
1

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
rafa.pereira
la source