ArcGIS Desktop – How to Speed Up Raster to Polygon Conversion in R?

arcgis-desktoppolygonizerraster

I've decided to process my Landsat data in R instead of ArcGIS – due to my missing knowledge of python and because of (assumed) high computation capacities of R. I want to :

  1. import r1 raster to R,
  2. import shp1 convert raster r1 to shp r.to.poly (dissolve = TRUE)
  3. intersect converter raster r.to.poly with my polygon shp1
  4. calculate area of every created polygon of intersected shp

Thus:

# read shp
shp <-readOGR(dsn = "C://...",
    layer = "m")

#read raster
r1<-raster("r1.tif")

# convert raster to polygon, dissolved neighboring same values
r.to.poly<-rasterToPolygons(r1, dissolve = T)

# define the same projection 
proj4string(shp) <- proj4string(r.to.poly)  

# use intersection from raster package
int.r <-raster::intersect(r.to.poly,shp)  

# calculate area per polygon
int.r$area <-gArea(int.r, byid = T) 
  
# export shapefile
writeOGR(int.r, dsn = "C:/...",
          layer = "...", driver="ESRI Shapefile", overwrite = TRUE)

So far, so good, but it takes about an hour to run the single conversion! moreover, when I tried FOR loop, my R on Windows crashed twice… It runs on mac, for the moment. Where the problem could be and how can I increase computation speed? Am I running out of R memory? The raster size on my disk is only 779 580 byte, size of shp is 1 729 532 bytes, thus are small. Also, make the same task in ArcGIS takes only couple seconds.

I've found some related discussion here: Increasing speed of crop, mask, & extract raster by many polygons in R? but as I have only about 10 rasters to process I don't want to start with parallel processing…

Best Answer

There is a "new" method from the stars package, which revolutionized the workflow for me (I was using the gdal_polygonizeR function previously). It has been faster than the John Baumgartner solution for all complicated rasters I have tried it on. Further, it does not require the gdal_polygonize.py script, which can be difficult to install on some machines. Try:

r.to.poly <- sf::as_Spatial(sf::st_as_sf(stars::st_as_stars(r1), 
                            as_points = FALSE, merge = TRUE)
                            ) # requires the sf, sp, raster and stars packages

See also this thread.

PS. If your polygons are complex, you probably need to validate them, before you can use them to calculate area:

rgeos::gIsValid(r.to.poly) # FALSE here means that you'll need to run the buffer routine:
r.to.poly <- rgeos::gBuffer(r.to.poly, byid = TRUE, width = 0)

If the polygons are very complex, you may need to add buffer width to get rid of all crossing edges:

r.to.poly <- rgeos::gBuffer(r.to.poly, byid = TRUE, width = 1000)
r.to.poly <- rgeos::gBuffer(r.to.poly, byid = TRUE, width = -1000)

Further, raster::area is quicker also for polygons than rgeos::gArea, I think (have not tested this properly, though).

Related Question