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

17

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
Roman Luštrik
quelle
1
Vielen Dank für die Bereitstellung eines reproduzierbaren Beispiels. Hilft immer beim Versuch, ein Problem zu verstehen!
fdetsch

Antworten:

21

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

Raumfahrer
quelle
Großartig, das funktioniert auch für mehrere Polygone (ich habe ein Beispiel hinzugefügt, mit dem ich in Joshuas Antwort spielen kann).
Roman Luštrik
Bei großen Polygon-Datasets ist die Umwandlung in ein SpatialPolygons-Objekt sehr aufwändig und nicht erforderlich. Wenn Sie "over" auf einen SpatialPolygonsDataFrame anwenden, wird der Zeilenindex zurückgegeben, mit dem die Punkte untergeordnet werden können. Siehe mein Beispiel unten.
Jeffrey Evans
Eine Menge von Overhead? Es wird im Wesentlichen nur der Slot @polygons aus dem SpatialPolygonsDataFrame-Objekt verwendet. Sie können es sogar "vortäuschen", indem Sie die Klasse eines SpatialPolygonsDataFrame "SpatialPolygons" zuweisen (obwohl dies hacky ist und nicht empfohlen wird). Alles, was die Geometrie verwenden wird, muss sowieso irgendwann diesen Schlitz bekommen, also ist es relativ gesehen überhaupt kein Overhead. Es ist ohnehin unbedeutend in jeder realen Anwendung, in der Sie dann eine Menge von Punkt-Polygon-Tests durchführen werden.
Spacedman
Bei der Berechnung des Overheads sind mehr als nur Geschwindigkeitsaspekte zu berücksichtigen. Beim Erstellen eines neuen Objekts im R-Namespace verwenden Sie den erforderlichen RAM. Wenn dies in kleinen Datenmengen kein Problem darstellt, wird die Leistung bei großen Datenmengen beeinträchtigt. R zeigt ein lineares Leistungsabsterben. Wenn die Daten größer werden, nimmt die Leistung ab. Wenn Sie kein zusätzliches Objekt erstellen müssen, warum dann?
Jeffrey Evans
1
Das wussten wir erst, als ich es gerade getestet habe.
Spacedman
13

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.

Edzer Pebesma
quelle
Gut zu wissen! Die Geometriemethode war mir nicht bekannt.
Jeffrey Evans
2
Willkommen auf unserer Seite, Edzer - schön, Sie hier zu sehen!
whuber
1
Danke Bill - auf stat.ethz.ch/pipermail/r-sig-geo wird es ruhiger , oder vielleicht sollten wir Software entwickeln, die mehr Ärger macht! ;-)
Edzer Pebesma
6

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) 
Jeffrey Evans
quelle
Falsch - die Rownamen des zurückgegebenen Objekts entsprechen dem Zeilenindex in_diesem_Fall - im Allgemeinen scheinen die Zeilennamen die Zeilennamen der Punkte zu sein - die möglicherweise nicht einmal numerisch sind. Sie können Ihre Lösung dahingehend modifizieren, dass eine Zeichenübereinstimmung durchgeführt wird, die sie möglicherweise etwas robuster macht.
Spacedman
@Sapcedman, sei nicht so dogmatisch. Die Lösung ist nicht falsch. Wenn Sie Punkte auf eine Menge von Polygonen setzen oder Punkte mit Polygonwerten versehen möchten, funktioniert die Over-Funktion ohne Zwang. Es gibt mehrere Möglichkeiten, den Zebrastreifen auszuführen, sobald Sie das resultierende Überobjekt haben. Ihre Lösung für das Erzwingen eines SpatialPolygon-Objekts verursacht einen erheblichen erforderlichen Overhead, da dieser Vorgang direkt für das SpatialPolygonDataFrame-Objekt ausgeführt werden kann. Übrigens, bevor Sie einen Beitrag bearbeiten, stellen Sie sicher, dass Sie Recht haben. Der Begriff Bibliothek und Paket werden in R.
Jeffrey Evans am
Ich habe meinem Beitrag einige Benchmarks hinzugefügt und ein weiteres Problem mit Ihrer Funktion entdeckt. "Pakete sind Sammlungen von R-Funktionen, Daten und kompiliertem Code in einem genau definierten Format. Das Verzeichnis, in dem Pakete gespeichert werden, wird als Bibliothek bezeichnet"
Spacedman,
Während Sie in Bezug auf "Paket" und "Bibliothek" technisch korrekt sind, argumentieren Sie mit der Semantik. Ich hatte gerade die Bitte eines Herausgebers von Ecological Modeling, die Verwendung von "package" (was eigentlich meine Präferenz ist) in "library" zu ändern. Mein Punkt ist, dass sie austauschbare Begriffe und eine Frage der Präferenz werden.
Jeffrey Evans
1
"Technisch korrekt", wie Dr. Sheldon Cooper einmal bemerkte, "ist die beste Art von richtig". Dieser Editor ist technisch falsch, was das Schlimmste ist.
Spacedman
4

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)
Josh O'Brien
quelle
Schlägt schrecklich fehl, wenn plymehr als ein Feature vorhanden ist, da gIntersectsfür jedes Feature eine Matrix mit einer Zeile zurückgegeben wird. Wahrscheinlich können Sie die Zeilen nach einem WAHREN Wert durchsuchen.
Spacedman
@Spacedman - Bingo. Müssen tun apply(gIntersects(pts, ply, byid=TRUE), 2, any). Tatsächlich werde ich die Antwort darauf umstellen, da sie auch den Fall eines einzelnen Polygons umfasst.
Josh O'Brien
Ah, any. Das ist vielleicht etwas schneller als die Version, die ich gerade getestet habe.
Spacedman
@Spacedman - Von meinen schnellen Tests, es sieht aus wie obrienund rowlings2Lauf Hals und Nacken, mit obrien vielleicht 2% schneller.
Josh O'Brien
@ JoshO'Brien wie kann man diese antwort auf viele polygone anwenden? Das heißt, ppsollte ein haben, IDdas angibt, in welchem ​​Polygon die Punkte liegen.
Code123
4

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
fdetsch
quelle
1
Sollte tmpsein pts.intersect? Das Parsen der zurückgegebenen Dimnames hängt auch vom undokumentierten Verhalten ab.
Spacedman
@Spacedman, Sie haben Recht tmp, vergessen, es zu entfernen, wenn der Code fertig ist. Sie haben auch Recht mit dem Parsen von dimnames. Dies war eine schnelle Lösung, um dem Fragesteller eine schnelle Antwort zu geben, und es gibt sicherlich bessere (und universellere) Ansätze, zum Beispiel Ihren :-)
fdetsch
1

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
rafa.pereira
quelle