13

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

ArcGIS 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
  • Did that, see Update 3. – Generic Wevers Oct 16 '15 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 '15 at 11:35
  • For such a huge raster, you'd be way better using C/C++ with gdal library. – Rodrigo Sep 17 at 20:41
7

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. – Generic Wevers Oct 16 '15 at 10:12
  • I didn't expect that ArcGis can be faster that gdal. @Generic Militzer – Iris Oct 16 '15 at 10:26
  • Ah wait, this will create polygons but I need polylines. – Generic Wevers Oct 16 '15 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. – Generic Wevers Oct 16 '15 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 '15 at 11:40
5

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)

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.

library(rgdal)
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. – Generic Wevers Oct 16 '15 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). – Generic Wevers Oct 16 '15 at 10:10
  • Haha, I definitely see your point. See above update. – fdetsch Oct 16 '15 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 '15 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 '15 at 11:42
5

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

plot(x)

enter image description here

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

enter image description here

2

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. – Generic Wevers Oct 19 '15 at 11:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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