arcgis desktop - How to speed up raster to polygon conversion in R?

I've decided to process my Landsat data in R instead or ArcGIS - due to my missing knowledge of python and 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 I 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...

Answer

Gdal Installation

Install Gdal command line tools and check to see if its binaries are added to path environment variable. e.g. in windows: open Run and type:

rundll32.exe sysdm.cpl,EditEnvironmentVariables 

Then follow the screenshot

Download and install gdal python bindings from here according to your python and OS.

install it using:

pip.exe install GDAL-2.0.2-cp27-none-win32.whl

You may encounter issues while installing gdal. Please see Installing GDAL with Python on windows?

In that thread users have suggested that Gdal binaries from gisinternals installs both the command line tools and the python bindings. Try to install it from there. Thus none of the above steps would be relevant.

To Make sure Gdal is installed open a command prompt and type:

ogrinfo

and check it works and you don't get:

'ogrinfo' is not recognized as an internal or external command, operable program or batch file.

To check gdal python bindings is installed open command prompt and type:

python
import gdal

if you get the following error then it is not installed:

Traceback (most recent call last): File "", line 1, in ImportError: No module named gdal

R side

Source the following r function from John Baumgartner blog:

gdal_polygonizeR <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
                             pypath=NULL, readpoly=TRUE, quiet=TRUE) {
  if (isTRUE(readpoly)) require(rgdal)
  if (is.null(pypath)) {
    pypath <- Sys.which('gdal_polygonize.py')

  }
  if (!file.exists(pypath)) stop("Can't find gdal_polygonize.py on your system.")
  owd <- getwd()
  on.exit(setwd(owd))
  setwd(dirname(pypath))
  if (!is.null(outshape)) {
    outshape <- sub('\\.shp$', '', outshape)
    f.exists <- file.exists(paste(outshape, c('shp', 'shx', 'dbf'), sep='.'))
    if (any(f.exists))
      stop(sprintf('File already exists: %s',

                   toString(paste(outshape, c('shp', 'shx', 'dbf'),
                                  sep='.')[f.exists])), call.=FALSE)
  } else outshape <- tempfile()
  if (is(x, 'Raster')) {
    require(raster)
    writeRaster(x, {f <- tempfile(fileext='.tif')})
    rastpath <- normalizePath(f)
  } else if (is.character(x)) {
    rastpath <- normalizePath(x)
  } else stop('x must be a file path (character string), or a Raster object.')

  system2('python', args=(sprintf('"%1$s" "%2$s" -f "%3$s" "%4$s.shp"',
                                  pypath, rastpath, gdalformat, outshape)))
  if (isTRUE(readpoly)) {
    shp <- readOGR(dirname(outshape), layer = basename(outshape), verbose=!quiet)
    return(shp)
  }
  return(NULL)
}

use this command to do the job:

r.to.poly<-gdal_polygonizeR(r1,pypath = "C:\\Program Files\\GDAL\\gdal_polygonize.py")#, dissolve = T)

Change pypath parameter according to your system. Be cautious that gdal_plygonize creates huge shapefiles. My 1 MB tif converted to a 128 MB shapefile. R needs a lot of memory to open this shapefile. Although the conversion was very fast. Thanks to python and gdal!

Another option would be Esri r-bridge to do the computation in Arcgis and return the output to R. However r-bridge doesn't support raster layers yet. (Thanks to @JeffreyEvans)

However the Gdal method is commercial free!