Fastest way to convert big raster to polyline using R or Python?

Question

I have a big raster file (129600 by 64800 pixel) with global water bodies (1 bit values 0 and 1) and try to extract ocean and inland water shorelines.

I've tried with ArcGIS and QGIS to convert from raster to polyline, but it takes ages.

Does anybody know a better/faster way (Python or R) or a better tool for this task?

Update

• R: rasterToContour might be fast and precise but if you have a very large dataset like mine (8,398,080,000 pixels) you need either a very big amount of RAM (more than 16GB) or you force R to do more processing on the hard drive and it will also take ages.
• Python/GDAL: gdal_poligonize creates polygons instead of polylines

Update 2

• R rasterToContour: rasterToContour does not deliver the wanted results. Compared to ArcGIS (raster to polygon followed by feature to line) it does not extract the exact pixel outline, as shown in the examples below.

rasterToContour result

ArcGIS result

UPDATE 3

Python/GDAL: I've run gdal_polygonize from command line against ArcGIS on a test dataset and the results were extremely clear:

• gdal: 49 seconds
• ArcGIS: 1.84 seconds

Answer 1

Edit: In the latest version of sf (v0.9.1) and stars (v0.4.1), you need to specify the merge = TRUE argument in the call to sf::st_as_stars() function.

For posterity, I've been having success with the stars:: package in R for doing this type of operation quickly.

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(merge = TRUE) %>% # this is the raster to polygons part
  st_cast("MULTILINESTRING") # cast the polygons to polylines

plot(x)

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

Answer 2

I'm working with R and used rasterToPolygons from the raster package in the past, but now I prefer gdal_polygonizeR by John Baumgartner. It bases on gdal_polygonize.py and is much faster. John Baumgartner published the code and gave an example for usage in his blog.

If you are familiar with python you could use gdal_polygonize.py directly of course.

Answer 3

Try rasterToContour from the raster package.

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)

You may then easily write the files to a local folder, e.g. as 'ESRI Shapefile' (.shp), using the below code. Have a look at ogrDrivers from rgdal to find out which drivers your system is compatible with.

library(rgdal)
writeOGR(x, dsn = getwd(), layer = "coastlines", driver = "ESRI Shapefile")

Answer 4

While I'm a big fan of GDAL, the polygonize tool was way too slow for my applications as well.

A fast alternative is gdal_trace_outline from Dans GDAL scripts which also has more options regarding tolerance, donuts, etc.

Like gdal_polygonize this also produces polygons which you'd need to convert afterwards with ogr2ogr -nlt MULTILINESTRING.

Downside to that is you need to compile it yourself, unless you are on a Linux or Mac OsX System.