Reparieren von verwaisten Löchern in R

Ich versuche, eine Vereinigung auf einem gemeinsamen Feld durchzuführen, nachdem zwei benachbarte Shapefiles zusammengeführt wurden. Zwischen den Shapefiles befindet sich mindestens ein dünner Streifen. Wenn ich versuche, eine Gewerkschaft zu gründen, wird der folgende Fehler angezeigt:

Fehler in createPolygonsComment (p): rgeos_PolyCreateComment: Verwaiste Bohrung, enthält kein Polygon für Bohrung bei Index 17

Ich habe unter diesem Link ein reproduzierbares Beispiel zu Dropbox hochgeladen .

Hier ist der Code, um das Problem neu zu erstellen:

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

Kehrt zurück:

Fehler in createPolygonsComment (p): rgeos_PolyCreateComment: Verwaiste Bohrung, enthält kein Polygon für Bohrung bei Index 17

Probieren Sie die hier und hier vorgeschlagene Lösung aus :

slot(example, "polygons") <- lapply(slot(example, "polygons"), checkPolygonsHoles)

Dies gibt den gleichen Fehler zurück, der vom Vereinigungsversuch stammt, jedoch mit einer anderen Indexnummer:

rgeos_PolyCreateComment: Verwaiste Bohrung, enthält kein Polygon für Bohrung bei Index 30

Versuchen Sie, das in Roger Bivands hilfreichem Tutorial vorgeschlagene Problem zu beheben

fix <- slot(example, "polygons")
fixa <- lapply(fix, checkPolygonsHoles)

Gibt den gleichen Fehler wie oben bei Index 30 zurück.

Andere haben dieses Problem hier und hier angesprochen , und obwohl die oben genannten Lösungen in einigen Fällen zu funktionieren scheinen, sind andere Fälle nicht gelöst. Ein Benutzer verwendete QGIS, um das Problem zu beheben, und der andere hatte 2 von 3 Punkten behoben, aber keine Lösung für den letzten.

Es scheint, dass die Leute weiterhin Probleme haben, obwohl dieser Code von Zeit zu Zeit funktioniert. Hat jemand eine Lösung in R gefunden?

Ich habe das Tool "Geometrie reparieren" in ArcGIS ausgeführt und das Problem behoben. Es scheint jedoch, dass in R eine Korrektur erforderlich ist.

Ich habe die Geometrieprobleme in den angehängten Daten analysiert und es scheint, dass dies nicht NUR der Fall ist, orphaned holessondern auch geometry validity issues. Es ist wahr, dass es sich um ein orphaned holeProblem mit der Gültigkeit der Geometrie handelt, aber RGEOS behandelt es nicht auf die gleiche Weise, wie bei verwaisten Löchern ein Fehler anstelle einer einfachen Warnung ausgegeben wird. Wie Sie angeben, handelt es sich hierbei um Hinweise zum Überprüfen von Polygonlöchern, die jedoch beim Anwenden zum Beheben von verwaisten Löchern nicht immer erfolgreich sind.

Also lasst uns:

1. Bereinigen Sie Ihre Daten (dies ist erforderlich, wenn Sie eine Geoverarbeitung wie eine Vereinigung durchführen möchten).

2. Verwenden Sie die bereinigten Daten für Ihren Gewerkschaftsprozess

1. Geometrie bereinigen Das Bereinigen von Geometrien in R kann manchmal schwierig sein. Daher habe ich versucht, ein experimentelles R-Paket (siehe https://github.com/eblondel/cleangeo ) zu erstellen , das das Bereinigen von spObjekten erleichtern soll ( derzeit nur eingeschränkt verfügbar ) polygonale Formen). Sie können das Paket installieren mit:

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

Zu Beginn ist es gut, dass Sie die Geometrieprobleme mit Ihren Quelldaten erkennen. Dazu können Sie Folgendes ausführen (Ihre Daten sind groß, sodass es einige Zeit dauern kann):

#get a report of geometry validity & issues for a sp spatial object
report <- clgeo_CollectionReport(sp)
summary <- clgeo_SummaryReport(report)
issues <- report[report$valid == FALSE,]

Damit werden Sie feststellen, dass Ihre Daten zwei Arten von Problemen aufweisen: orphaned holesund geometry validity issues. Beide (und nicht nur die verwaisten Löcher) führen wahrscheinlich dazu, dass der unionProzess fehlschlägt. Daher sollten die Daten nach Möglichkeit vorab automatisch bereinigt werden. Für eine schnelle Reproduktion wird im folgenden ersten Beispielcode nur die Teilmenge der als verdächtig gekennzeichneten Funktionen verwendet (mit Ausnahme der neuesten mit dem Index = 9002 in den Originaldaten - siehe meinen Hinweis unten).

#get suspicious features (indexes)
nv <- clgeo_SuspiciousFeatures(report)
mysp <- sp[nv[-14],]

#try to clean data
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

Wenn clgeo_Cleandas gut geht, sollten Sie alle Geometrien jetzt gültig bekommen. Sie können dies auf den gesamten Datensatz anwenden (außer Merkmalsindex = 9002).

#try to clean data
mysp <- sp[-9002,]
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

2. Vereinigungsprozess Nun wollen wir sehen, ob der Vorgangunion mit diesem Datensatz funktioniert:

#Attempting a UnionSpatialPolygons based on the COUNTY field
mysp.df <- as(mysp, "data.frame")
countycol <- mysp.df$COUNTY
mysp.diss <- unionSpatialPolygons(mysp.clean, countycol)

Hinweis: Wie bereits erwähnt, habe ich ein Feature entfernt (Index = 9002). Wenn Sie es zeichnen:, werden Sie feststellen, plot(sp[9002,])dass dieses Feature sehr (sehr) komplex ist. Ich habe es nur aus der Stichprobe ausgeschlossen, weil das Überprüfen der Löcher zu lange gedauert hat. Mal sehen, ob das gleiche Problem beim Lesen der Daten mit readShapePoly(from maptools) auftritt ...

3. Wechseln Sie zu readShapePoly vs. readOGR, um Daten zu lesen (UPDATE).

readOGRist nicht die einzige Funktion zum Lesen von Shapefiles. Sie können auch verwenden , readShapePolyvon maptoolsPaket, in der Regel performanter als die erste:

require(maptools)
mysp <- readShapePoly("ReproducibleExample.shp")

Abgesehen von schneller laufen:

• Wenn Sie den obigen Code verwenden, der auf basiert clgeo_CollectionReport, gibt es kein Problem mit verwaisten Löchern, aber immer noch Probleme mit der Geometrie.

• Das Bereinigen der Geometrie mit clgeo_Cleanläuft ebenfalls gut und bleibt nun nicht mehr beim Feature-Index 9002 hängen

• Und ... der Gewerkschaftsprozess funktioniert.

Siehe unter dem Diagrammergebnis:

#plot the result
plot(mysp, border= "lightgray")
plot(mysp.diss, border="red", add = TRUE)

Fazit : lieber MapTools Ihre Shape - Datei Daten lesen und prüfen , mit cleangeo Ihre Daten vor jedem Geoprocessing zu reinigen.

Eine praktische Lösung, die für mich in R weiterhin funktioniert, ist das Anwenden eines Puffers mit der Breite Null :

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#project your data (I'm using California Albers here) and apply a zero-width buffer
example <- spTransform(example, CRS("+init=epsg:3310"))
example <- gBuffer(example, byid = T, width = 0)

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

unionSpatialPolygons dauert eine Weile mit diesem Datensatz, scheint aber gut zu funktionieren.