Wie überlagere ich SpatialPointsDataFrame mit einem Polygon und behalte die SPDF-Daten bei?

Ich habe eine SpatialPointsDataFramemit einigen zusätzlichen Daten. Ich möchte diese Punkte in einem Polygon extrahieren und gleichzeitig das SPDFObjekt und die entsprechenden Daten beibehalten .

Bisher hatte ich wenig Glück und habe versucht, eine gemeinsame ID abzugleichen und zusammenzuführen. Dies funktioniert jedoch nur, weil ich Daten mit individuellen IDs gerastert habe.

Hier ist ein kurzes Beispiel, ich suche nach Punkten innerhalb des roten Quadrats.

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

Der naheliegendste Ansatz wäre, overdie Daten aus dem Polygon zurückzugeben.

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

Aus der sp::overHilfe:

 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.

Wenn Sie also Ihre SpatialPolygonsDataFramein konvertieren , erhalten SpatialPolygonsSie einen Vektor von Indizes zurück, und Sie können Ihre Punkte auf NAfolgende Werte setzen :

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

Für die Zweifler ist hier der Beweis, dass der Konvertierungsaufwand kein Problem ist:

Zwei Funktionen - zuerst Jeffrey Evans 'Methode, dann mein Original, dann meine gehackte Konvertierung, dann eine Version, gIntersectsdie auf Josh O'Briens Antwort basiert:

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,]
}

Für ein reales Beispiel habe ich nun einige zufällige Punkte über den columbusDatensatz verteilt:

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

Sieht gut aus

plot(columbus)
points(pts)

Überprüfen Sie, ob die Funktionen dasselbe tun:

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

Und 500 Mal für das Benchmarking ausführen:

> 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 

Meiner Intuition nach ist dies kein großer Aufwand, sondern möglicherweise ein geringerer Aufwand, als alle Zeilenindizes in Zeichen und zurück zu konvertieren oder na.omit auszuführen, um fehlende Werte abzurufen. Was übrigens zu einem anderen Fehlermodus der evansFunktion führt ...

Wenn eine Zeile des Polygondatenrahmens alle ist NA(was vollkommen gültig ist), dann erzeugt die Überlagerung mit SpatialPolygonsDataFramefür Punkte in diesem Polygon einen Ausgabedatenrahmen mit allen NAs, der evans()dann abfällt:

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

ABER gIntersectsist schneller, selbst wenn die Matrix gewobbelt werden muss, um Schnittpunkte in R statt in C-Code zu prüfen. Ich vermute, es liegt an den prepared geometryFähigkeiten von GEOS, räumliche Indizes zu erstellen - ja, prepared=FALSEes dauert etwas länger, ungefähr 5,5 Sekunden.

Ich bin überrascht, dass es keine Funktion gibt, um die Indizes oder die Punkte direkt zurückzugeben. Als ich vor splancs20 Jahren schrieb, hatten die Point-in-Polygon-Funktionen beide ...

sp Bietet eine kürzere Form zum Auswählen von Features basierend auf der räumlichen Überschneidung gemäß dem OP-Beispiel:

pts[ply,]

ab:

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

Hinter den Kulissen ist dies eine Abkürzung für

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

Zu beachten ist, dass es eine geometryMethode gibt, die Attribute löscht: overÄndert das Verhalten, wenn das zweite Argument Attribute enthält oder nicht (dies war die Verwirrung von OP). Dies funktioniert accross alle Spatial * Klassen sp, obwohl einige overMethoden erfordern rgeos, sehen diese Vignette für Details, zB bei mehrfacher Übereinstimmungen für überlappende Polygone.

Sie waren auf dem richtigen Weg mit über. Die Rownamen des zurückgegebenen Objekts entsprechen dem Zeilenindex der Punkte. Sie können Ihren genauen Ansatz mit nur wenigen zusätzlichen Codezeilen implementieren.

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) 

Ist es das was du suchst?

Ein Hinweis zum Bearbeiten: Der Wrapping-Aufruf für apply()wird benötigt, damit dies mit beliebigen SpatialPolygonsObjekten funktioniert , die möglicherweise mehr als ein Polygon-Feature enthalten. Vielen Dank an @Spacedman, der mich dazu gebracht hat, zu zeigen, wie dies auf den allgemeineren Fall angewendet werden kann.

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)

Hier ist ein möglicher Ansatz mit dem rgeosPaket. Grundsätzlich wird die gIntersectionFunktion verwendet, mit der Sie zwei spObjekte schneiden können. Durch Extrahieren der IDs der Punkte, die innerhalb des Polygons liegen, können Sie anschließend Ihr Original SpatialPointsDataFrameunter Beibehaltung aller entsprechenden Daten unterteilen. Der Code ist fast selbsterklärend, aber wenn Sie Fragen haben, können Sie diese gerne stellen!

# 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

Es gibt eine äußerst einfache Lösung für die Verwendung der spatialEcoBibliothek.

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)

Überprüfen Sie das Ergebnis:

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