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

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

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.

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)

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

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)

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

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.