Polygon abschneiden und Daten behalten?

13

Ich habe diese zwei Polygone:

library(sp); library(rgeos); library(maptools)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

Bildbeschreibung hier eingeben

Ich möchte Teile davon ausschneiden, die von spdf1geschnitten werden spdf2. Ich möchte spdf1jedoch als SpatialPolygonsDataFrame bleiben und alle darin enthaltenen Informationen behalten spdf1@data.

Ich habe gDifference wie folgt ausprobiert, wobei Teile spdf1, die von geschnitten werden spdf2, ausgeschnitten werden , aber dann spdf1in SpatialPolygons konvertiert werden (dh die darin enthaltenen Informationen verworfen werden spdf1@data).

gDifference(spdf1, spdf2, byid=T)

Wie kann ich schneiden in spdf1mitspdf2 den darin enthaltenen Daten arbeiten, diese aber behalten spdf1@data? Ich habe diese ähnliche Frage geprüft und ausprobiert, ohne ein Polygon über SpatialPointsDataFrame zu legen und die SPDF-Daten beizubehalten.

Luciano
quelle

Antworten:

9

Der einfachste Ansatz wäre

  library(raster)
  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

Weitere Funktionen für die Überlagerung von Polygonen finden Sie unter? 'Raster-package' (Abschnitt XIV). Diese Funktionen verwenden die Basisfunktionen von RGEOS unter der Haube in Funktionen auf Benutzerebene (im Gegensatz zu Funktionen auf Entwicklerebene).

Robert Hijmans
quelle
Was verstehen Sie unter "Funktionen auf Benutzerebene (im Gegensatz zu Funktionen auf Entwicklerebene)"?
Luciano
rgeosBietet geometrische Operationen, befasst sich jedoch nicht mit den Attributen der Daten. Daher erfordert die Verwendung dieser Funktionen viel Arbeit, um alles zusammenzuhalten. Die Rasterfunktionen vereinfachen dies und verhalten sich wie ähnliche Funktionen in der GIS-Software,
Robert Hijmans
Ja, aber dies kann zu einem Fehler in SpatialPolygonsDataFrame (Teil2, x @ data [match (Zeilennamen (Teil2),: Zeilennamen von Daten und Polygon-IDs stimmen nicht überein
jebyrnes
Das wäre ein Fehler. Können Sie ein Beispiel dafür geben?
Robert Hijmans
4

Eine Problemumgehung wäre, die Attribute nach dem Ausführen des Clips erneut hinzuzufügen, während Sie von SpatialPolygonsnach konvertieren SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)

spdf3@data

   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

Bildbeschreibung hier eingeben

Andre Silva
quelle
Das sieht nach einem Problem aus, das ich habe. Nur der Ausgabeclip in meiner speziellen Instanz erzeugt Rownamen aus spdf1, die nicht existieren (als einfaches Gsub, um die 2. Ziffer in der Reihe loszuwerden. Names sollte sich um die Übereinstimmung der Rownamen kümmern, nicht wahr? )
jebyrnes
Fehler in SpatialPolygonsDataFrame (clip, data = as.data.frame (all_spdfs_together @ data)): Objektlängeninkongruenz: clip hat 1718 Polygonobjekte, aber as.data.frame (all_spdfs_together @ data) hat 86 Zeilen
jebyrnes
Klar - ich habe ein paar Polygone von Dingen im Ozean. Einige sind falsch auf dem Land platziert oder überschneiden sich mit dem Land. Ich möchte das ausschneiden, damit ich nur Gebiete im Ozean habe. Ich habe ein Shapefile an der Küste, mit dem ich vergleichen kann. Hier ist der Code - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
jebyrnes
1
nevermind - die raster :: erase lösung funktioniert (es hatte aus irgendeinem
grund