Najszybszy sposób na konwersję dużego rastra na polilinię za pomocą R lub Pythona?

14

Mam duży plik rastrowy (129600 na 64800 pikseli) z globalnymi częściami wód (1 bit wartości 0 i 1) i próbuję wydobyć linie brzegowe oceanów i wód śródlądowych.

Próbowałem z ArcGIS i QGIS przekonwertować z rastra na polilinię, ale to trwa wieki.

Czy ktoś zna lepszy / szybszy sposób (Python lub R) lub lepsze narzędzie do tego zadania?

Aktualizacja

  • R: rasterToContour może być szybki i precyzyjny, ale jeśli masz bardzo duży zestaw danych, taki jak mój (8 398,080 000 pikseli), potrzebujesz bardzo dużej ilości pamięci RAM (ponad 16 GB) lub zmuszasz R do większego przetwarzania na dysku twardym i to zajmie także wieki.
  • Python / GDAL: gdal_poligonize tworzy wielokąty zamiast polilinii

Aktualizacja 2

  • R rasterToContour: rasterToContour nie zapewnia oczekiwanych rezultatów. W porównaniu do ArcGIS (raster do wielokąta, a następnie funkcja do linii) nie wyodrębnia dokładnego obrysu pikseli, jak pokazano w poniższych przykładach.

wynik rasterToContour wynik rasterToContour

Wynik ArcGIS Wynik ArcGIS

AKTUALIZACJA 3

Python / GDAL: Uruchomiłem gdal_polygonize z wiersza poleceń przeciwko ArcGIS na testowym zbiorze danych, a wyniki były bardzo jasne:

  • gdal: 49 sekund
  • ArcGIS: 1,84 sekundy
python r performance Ogólne Wevers
źródło
Zrobił to, patrz aktualizacja 3.
Generic Wevers
Czy możesz podać ten zestaw danych testowych, abyśmy mogli sprawdzić, czy proponowane alternatywy są szybsze i / lub wygenerować wymagane wyniki?
Kersten
Dla tak ogromnego rastra lepiej byłoby używać C / C ++ z biblioteką gdal.
Rodrigo,

Odpowiedzi:

8

Pracuję z R i korzystałem rasterToPolygonsz rasterpakietu w przeszłości, ale teraz wolę gdal_polygonizeRJohna Baumgartnera. Opiera się na gdal_polygonize.pyi jest znacznie szybszy. John Baumgartner opublikował kod i podał przykład zastosowania na swoim blogu .

Jeśli znasz język Python, możesz gdal_polygonize.pyoczywiście używać go bezpośrednio.

Irys
źródło
1
Dam to, próbuję. Ostatnim razem, gdy użyłem gdal_polygonize.py, ArcGIS był jeszcze szybszy.
Generic Wevers
Nie spodziewałem się, że ArcGis może być szybszy niż ten. @Generic Militzer
Iris
Ach, czekaj, to stworzy wielokąty, ale potrzebuję polilinii.
Generic Wevers
Jeśli umieścisz swoje dane w Geobazie Plików, jest to dość szybkie. Ale wciąż za mało. Dlatego szukam alternatyw.
Generic Wevers
2
Otrzymanie wielokątów niekoniecznie jest problemem, zawsze możesz je przekonwertować na polilinie (chociaż przy tak wielu może to również trochę potrwać).
Martin
7

Jeśli chodzi o potomstwo, odniosłem sukces dzięki stars::pakietowi Rumożliwiającemu szybkie wykonanie tego typu operacji.

library(raster)
library(stars)
library(sf)
library(magrittr)

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r[r[] < 750] <- 0
r[r[] >= 750] <- 1

x <- st_as_stars(r) %>% 
  st_as_sf() %>% # this is the raster to polygons part
  st_cast("MULTILINESTRING") # cast the polygons to polylines

plot(x)

wprowadź opis zdjęcia tutaj

plot(r)
plot(x, add = TRUE)

wprowadź opis zdjęcia tutaj

mikoontz
źródło
5

Spróbuj rasterToContourz pakietu rastrowego .

f <- system.file("external/test.grd", package="raster")
r <- raster(f)
r[r[] < 750] <- 0
r[r[] >= 750] <- 1

x <- rasterToContour(r)
class(x)
> [1] "SpatialLinesDataFrame"
> attr(,"package")
> [1] "sp"

plot(r)
plot(x, add=TRUE)

wprowadź opis zdjęcia tutaj

Następnie możesz łatwo zapisać pliki w folderze lokalnym, np. Jako „Plik kształtu ESRI” (.shp), używając poniższego kodu. Przyjrzeć się ogrDriversz rgdal , aby dowiedzieć się, jakie sterowniki system jest kompatybilny z.

library(rgdal)
writeOGR(x, dsn = getwd(), layer = "coastlines", driver = "ESRI Shapefile")
fdetsch
źródło
Spróbuję trzymać kciuki, aby nie zabić mojej pamięci RAM. Mimo że mam 16 GB, co, mam nadzieję, wystarczy, R czasami nie jest tak wydajny w przypadku dużych plików rastrowych. Ale zobaczmy.
Generic Wevers
Konwersja jakoś działała, ale nie byłem w stanie szczegółowo sprawdzić. Ponieważ zwykle bardziej interesuję się przetwarzaniem danych rastrowych, czy możesz mi powiedzieć, jak mogę przenieść SpatialLineDataFrame do pliku kształtu lub czegoś porównywalnego. Mam google i wciąż walczę, ponieważ nie znam nazwy warstwy (OGRwrite).
Generic Wevers
Haha, zdecydowanie rozumiem twój punkt widzenia. Zobacz powyższą aktualizację.
fdetsch
2
Kolejna wskazówka: spróbuj ustawić „maxpixels” rasterToContourna jakąś wyższą wartość, np. 1e + 9. Otrzymasz wtedy więcej szczegółów. Ustawienie domyślne tworzy dość ogólne linie konturu.
fdetsch
1
Jeśli nie chcesz, aby resampleTwoje dane były z grubszą rozdzielczością przestrzenną, jedynym rozwiązaniem, jakie mogę sobie wyobrazić, byłoby podzielenie danych na wiele kafelków (np. 16 sub-rastrów), a następnie wykonanie rasterToContourkażdej płytki osobno w sposób iteracyjny i , w końcu mergepowstałe pliki kształtów w jeden ogromny plik kształtów. Jeśli jesteś zainteresowany, pakiet naszej grupy roboczej Rsenal oferuje funkcję zwaną splitRastertworzeniem wielu podrzędnych rastrów z jednego ogromnego rastra.
fdetsch
2

Chociaż jestem wielkim fanem GDAL, narzędzie do poligonizacji było zbyt wolne dla moich aplikacji.

Szybkim rozwiązaniem jest gdal_trace_outlineze skryptami Dans gdal który ma również więcej opcji w zakresie tolerancji, pączki itp

W gdal_polygonizeten sposób powstają również wielokąty, które należy później przekonwertować ogr2ogr -nlt MULTILINESTRING.

Minusem tego jest to, że musisz go skompilować samodzielnie, chyba że korzystasz z systemu Linux lub Mac OSX.

Kersten
źródło
Niestety nie powiodło się z komunikatem o błędzie: „Błąd segmentacji (zrzut rdzenia)”. Zgaduję, że mój plik jest zbyt duży lub bardziej precyzyjny, że wytworzy zbyt wiele małych wielokątów.
Generic Wevers,