Naprawianie osieroconych otworów w R.

18

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.

Luke Macaulay
źródło
Bez twoich danych trudno jest powiedzieć, gdzie jest problem.
@Pascal, właśnie przesłałem link do skrzynki odbiorczej ze zmniejszonym plikiem shapefile 10mb spakowanym i rozpakowanym 16mb, który odtworzy problem. Nie byłem pewien, jak dostarczyć dane, ponieważ oryginał miał 1,5 GB, ale udało mi się użyć ArcGIS, aby zawęzić problem do mniejszego pliku. Czy istnieje dobry protokół do tworzenia i udostępniania powtarzalnych przykładów w zarządzalnej wielkości?
Luke Macaulay
Próbowanie różnych podejść z R nie działało. Qgis marznie podczas sprawdzania geometrii.

Odpowiedzi:

25

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)

Wynik unijny

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

eblondel
źródło
Dzięki eblondel! Spróbuję tego. Dziękujemy za opracowanie pakietu!
Luke Macaulay
Cześć eblondel, To działało dobrze, ale chciałem poinformować cię, że przy poprawianiu geometrii często tworzyłby bardzo duży wielokąt w przypadku skomplikowanych i złożonych operacji. Na przykład sieć drogowa została skorygowana do dużego wielokąta, który był w zasadzie zasięgiem sieci. Nie jestem pewien, jak łatwo to naprawić, ale chciałem cię poinformować.
Luke Macaulay
Łał. Bardzo imponujący pakiet. To musiało być dużo pracy.
nograpes
3
Dzięki @nograpes za twoją opinię. Zbudowałem ten pakiet od zera, gdy ten problem został opublikowany, również dlatego, że czyszczenie geometrii nie zawsze jest łatwym zadaniem. Jeśli korzystasz z Github, witam twoją „gwiazdkę” :-), chciałbym w przyszłości ulepszyć pakiet i ewentualnie wydać go w CRAN.
eblondel,
7
Aby poinformować, że cleangeo zostało opublikowane w CRAN ( cran.r-project.org/package=cleangeo ), wszystkim osobom, które go używają, możesz zgłaszać prośby o ulepszenia lub błędy w Github.
eblondel
1

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.

FGG
źródło