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?


  • 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 rasterToContour result

ArcGIS result ArcGIS result


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
  • Did that, see Update 3. Oct 16, 2015 at 13:06
  • Can you provide that test dataset, so we can see if proposed alternatives are faster and/or produce the required results?
    – Kersten
    Oct 17, 2015 at 11:35
  • For such a huge raster, you'd be way better using C/C++ with gdal library.
    – Rodrigo
    Sep 17, 2019 at 20:41

4 Answers 4


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.


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


enter image description here

plot(x, add = TRUE)

enter image description here


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.

  • 1
    I'll give it I try. Last time I've used gdal_polygonize.py ArcGIS was still faster. Oct 16, 2015 at 10:12
  • I didn't expect that ArcGis can be faster that gdal. @Generic Militzer
    – Iris
    Oct 16, 2015 at 10:26
  • Ah wait, this will create polygons but I need polylines. Oct 16, 2015 at 10:28
  • If you put your data into a File Geodatabase it is pretty fast. But still not fast enough. That's why I'm searching for alternatives. Oct 16, 2015 at 10:30
  • 2
    It's not necessarily a problem that you get polygons, you can always convert them to polylines (although, with that many it could of course take a while as well).
    – Martin
    Oct 16, 2015 at 11:40

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)
> [1] "SpatialLinesDataFrame"
> attr(,"package")
> [1] "sp"

plot(x, add=TRUE)

enter image description here

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.

writeOGR(x, dsn = getwd(), layer = "coastlines", driver = "ESRI Shapefile")
  • I'll try and keep fingers crossed it will not kill my RAM. Even though I have 16GB, which hopefully is enough, R sometimes is not so efficient with big raster files. But let's see. Oct 16, 2015 at 9:06
  • Conversion worked somehow, but I wasn't able to check in detail. As I am usually more into raster data processing, can you tell me how I can transfer the SpatialLineDataFrame into a shapefile or something comparable. I've googled and still struggling, as I don't know the layer name (OGRwrite). Oct 16, 2015 at 10:10
  • Haha, I definitely see your point. See above update.
    – fdetsch
    Oct 16, 2015 at 10:25
  • 2
    Another hint: try to set 'maxpixels' in rasterToContour to some higher value, e.g. 1e+9. You'll end up with more details then. The default setting creates quite generalized contour lines.
    – fdetsch
    Oct 16, 2015 at 10:46
  • 1
    If you're not willing to resample your data to a coarser spatial resolution, the only solution I can imagine then would be to split your data into multiple tiles (e.g. 16 sub-rasters), then perform rasterToContour on each tile separately in an iterative manner and, finally, merge the resulting shapefiles into one huge shapefile. In case you are interested, our working group's package Rsenal offers a function called splitRaster to create multiple sub-rasters from one huge raster.
    – fdetsch
    Oct 16, 2015 at 11:42

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.

  • Unfortunately it failed with the error message: "Segmentation fault (core dumped)". I am guessing my file is too big or more precise it will produce too many small polygons. Oct 19, 2015 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.