Naprawianie osieroconych otworów w R.

Próbuję wykonać połączenie na wspólnym polu po scaleniu dwóch sąsiednich plików kształtów. Pliki kształtów kończą się co najmniej jednym cienkim kawałkiem przestrzeni między nimi. Podczas próby połączenia pojawia się następujący błąd osieroconej dziury:

Błąd w createPolygonsComment (p): rgeos_PolyCreateComment: osierocona dziura, nie można znaleźć zawierającego wielokąt dla dziury o indeksie 17

Pod tym linkiem przesłałem do Dropbox odtwarzalny przykład .

Oto kod, aby odtworzyć problem:

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

Zwroty:

Błąd w createPolygonsComment (p): rgeos_PolyCreateComment: osierocona dziura, nie można znaleźć zawierającego wielokąt dla dziury o indeksie 17

Wypróbowanie poprawki zaproponowanej tutaj i tutaj :

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

Zwraca ten sam błąd, który pochodzi z próby połączenia, ale z innym numerem indeksu:

rgeos_PolyCreateComment: osierocona dziura, nie można znaleźć zawierającego wielokąta dla dziury o indeksie 30

Wypróbowanie poprawki zaproponowanej w pomocnym samouczku Rogera Bivanda

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

Zwraca ten sam błąd przy indeksie 30 jak powyżej.

Inni podnieśli ten problem tu i tutaj , i chociaż powyższe rozwiązania wydają się działać w niektórych przypadkach, inne nie zostały rozwiązane. Jeden użytkownik użył QGIS, aby rozwiązać problem, a drugi miał naprawione 2 z 3 elementów, ale brak ostatecznego rozwiązania.

Wygląda na to, że ludzie nadal mają problemy pomimo tego, że kod działa od czasu do czasu. Czy ktoś znalazł rozwiązanie w R?

W ArcGIS wykonałem narzędzie „napraw geometrię” i to rozwiązało problem, ale wygląda na to, że powinna być poprawka w R.

Przeanalizowałem problemy z geometrią w dołączonych danych i wydaje się, że nie tylko, orphaned holesale także geometry validity issues. To prawda, że ​​an orphaned holejest w pewnym sensie kwestią poprawności geometrii, ale rgeos nie traktuje tego w taki sam sposób, jak w przypadku osieroconych otworów, zamiast zwykłego ostrzeżenia pojawia się błąd. Jak wskazujesz, są one wskazówkami, aby sprawdzić otwory wielokąta, ale nie zawsze udaje się, gdy są stosowane w celu naprawy osieroconych otworów.

Więc:

1. wyczyść swoje dane (co jest wymagane, jeśli chcesz wykonywać geoprzetwarzanie jak związek)

2. użyj oczyszczonych danych w procesie związkowym

1. Czyszczenie geometrii Naprawianie geometrii w R może być czasem trudne, dlatego próbowałem zbudować eksperymentalny pakiet R (patrz https://github.com/eblondel/cleangeo ), który ma na celu ułatwienie czyszczenia spobiektów (obecnie ograniczone do kształty wielokątne). Możesz zainstalować pakiet za pomocą:

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

Na początek dobrze jest zobaczyć, jakie są problemy geometrii z danymi źródłowymi. W tym celu możesz wykonać następujące czynności (twoje dane są duże, więc może to zająć trochę czasu):

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

Dzięki temu zobaczysz, że twoje dane mają 2 rodzaje problemów: orphaned holesi geometry validity issues. Obie (i nie tylko osierocone dziury) prawdopodobnie spowodują, że unionproces się nie powiedzie, więc dane powinny być wcześniej wyczyszczone, w miarę możliwości w sposób zautomatyzowany. W celu szybkiego odtworzenia pierwszy przykładowy kod poniżej bierze tylko podzbiór funkcji, które są oznaczone jako podejrzane (z wyjątkiem najnowszej, z indeksem = 9002 w oryginalnych danych - zobacz moją notatkę poniżej na ten temat)

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

Jeśli clgeo_Cleandobrze spełnia swoje zadanie, wszystkie geometrie powinny być teraz prawidłowe. Możesz zastosować to do pełnego zestawu danych (z wyjątkiem indeksu funkcji = 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. Proces Unii Teraz zobaczmy, czy uniondziała na tym zestawie danych:

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

Uwaga: jak powiedziano wcześniej, usunąłem jedną funkcję (indeks = 9002) .Planując ją: plot(sp[9002,])zobaczysz, że ta funkcja jest bardzo (bardzo) złożona. Wyłączyłem go z próbki tylko dlatego, że sprawdzanie otworów zajmowało zbyt dużo czasu. Zobaczmy teraz, czy ten sam problem występuje przy użyciu readShapePoly(z maptools) do odczytu danych ...

3. Przejdź do readShapePoly vs. readOGR do odczytu danych (UPDATE)

readOGRnie jest jedyną funkcją dostępną do odczytu plików kształtów. Możesz także użyć readShapePolyz maptoolspakietu, ogólnie bardziej wydajnego niż pierwszy:

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

Oprócz szybszego działania:

• jeśli użyjesz powyższego kodu opartego na clgeo_CollectionReport, nie ma problemu z osieroconymi otworami, ale nadal problemy z geometrią.

• Czyszczenie geometrii za pomocą clgeo_Cleanrównież działa dobrze, a teraz nie zacina się z indeksem funkcji 9002

• I ... proces związkowy działa.

Zobacz poniżej wynik wykresu:

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

Wniosek : wolą maptools, aby odczytać dane pliku kształtu , i rozważ użycie cleangeo do czyszczenia danych przed jakimkolwiek przetwarzaniem geoprzestrzennym .

Wygodnym rozwiązaniem, które nadal działa dla mnie w R, jest zastosowanie bufora o zerowej szerokości :

#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 zajmuje trochę czasu z tym zestawem danych, ale wydaje się, że działa dobrze.