# 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")

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...

• If the speed isn't acceptable perhaps try GDAL_Polygonize gdal.org/gdal_polygonize.html as a process.. or as you have access to ArcGis but lack experience in arcpy perhaps I could help with that instead. – Michael Stimson May 8 '16 at 23:31
• thank you @MichaelMiles-Stimson, I've solved in with ModelBuilder for a while, it was just 11 rasters. But I will really appreciate your help in the future ! ;) – maycca May 9 '16 at 18:06
• Which step is taking the most time? Converting to polygons, or the intersection? – amball May 23 '16 at 4:47
• Also, are you only interested in the part of your raster that intersects with your shapefile? If so, you might try using raster::intersect to crop your raster BEFORE converting it to polygons. – amball May 23 '16 at 4:55

### 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 (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)))
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!

• Have you used the ESRI R-bridge? In previous versions it was a one-way street from ArcGIS and not designed to call ArcGIS processes from R. – Jeffrey Evans May 23 '16 at 18:51
• no, I didn't use R-bridge. canI use it just like an another library? and, moreover, I can avoid all of the GDAL installation? moreover, I am sure that i have python on my PC, as I hve GRASS and QGIS... – maycca May 23 '16 at 20:20
• I'm sorry @FaridCher.. I can't follow it! where can I put 'path = %path%;C:\Program Files\GDAL'? in R? I/ve downloaded the GDAL-2.0.2-cp35-none-win32.whl, and also I've installed OSGeo4W, but I am sure that 've already had this because I am using QGIS... also, when I tried to install pip.exe, I couldn't find the file ! please, is your code chronological, or there is different order of steps? thank you ! – maycca May 23 '16 at 20:30
• @JeffreyEvans I haven't tried it myself just a suggestion on top of my head. I checked the latest version yet native arcgis raster datasets are not supported. However one could use `Arcpy` and create a python script like `gdal_polygonizeR.py` and then call it from r and use the supported raster formats by R as the output. – Farid Cheraghi May 23 '16 at 22:58
• Thank you a lot @FaridCher, it works well now!! My problem was, that I was ustaling it while RStudio was running, and I was testing it on R. After Closing all programs and assigning path to MAC bash: `echo 'export PATH=/Library/Frameworks/GDAL.framework/Programs:\$PATH' >> ~/.bash_profile` and defining my `pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py"` after sourcing the gdal_polygonizeR it worked well. Thank you again for all help ! My 30 minutes needed to convert raster to polygon takes now only 3 seconds !! ;) – maycca May 27 '16 at 18:19

I used rasterToPolygons from the raster package too 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.

I'm copying these things in case link be broken:

## Installing GDAL on Mac OSX and Polygonize Rasters

GDAL Complete: you must install both gdal.pkg and numpy.pkg

In all cases:

1.1. Double click and you will get a *.dmg file in Downloads (or Descargas). Wait until the download completes (GDAL Complete will take longer).

1.2. Once the dmg is downloaded, it normally gets uncompressed by itself and the folder is automatically opened. Otherwise, go to Downloads and double-click on the dmg file.

1.3. Look for the *.pkg inside the folder that opens, double click, accept the license and proceed. Sometimes you have to click on the path for the installation in order to get the Accept button active.

Once you have all 4 packages installed,

1. Open a Terminal (Applications/Utilities/Terminal) and write export PATH=/Library/Frameworks/GDAL.framework/Programs:\$PATH and hit Return

READY TO GO TO R! :)

``````library(raster)
library(rgeos)
``````

Define the function aavailable from https://johnbaumgartner.wordpress.com/2012/07/26/getting-rasters-into-shape-from-r/

``````gdal_polygonizeR <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
pypath=NULL, readpoly=TRUE, quiet=TRUE) {
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)))
shp <- readOGR(dirname(outshape), layer = basename(outshape), verbose=!quiet)
return(shp)
}
return(NULL)
}
``````

make raster polygonization on dummy raster

## - the way you define the path may be different on R, RStudio and R from Terminal !! (at least on OSX 10.9.5)

``````# read raster
r5<-raster("r5.img")
plot(r5)

# in R + RStudio: define whole path where to find gdal_polygonize.py !!
r.to.poly<-gdal_polygonizeR(r5,pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py")#

#in R run from Terminal:  no need to define the whole path !!
r.to.poly<-gdal_polygonizeR(r5)
``````

to check, if you need to define the path to gdal_polygonize.py or not:

``````> Sys.which("gdal_polygonize.py")
gdal_polygonize.py
"/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py" # means path is defined, you should be able just to run

r.to.poly<-gdal_polygonizeR(r5)
``````

If it looks like this:

``````> Sys.which("gdal_polygonize.py")
gdal_polygonize.py
""
``````

you would need to define the whole path to access gdal_polygonize.py

``````r.to.poly<-gdal_polygonizeR(r5,pypath = "/Library/Frameworks/GDAL.framework/Programs/gdal_polygonize.py")
``````

You can do the analysis described in your post without converting the raster to a polygon. Use the `raster::extract` function to extract the raster values to each polygon. You can then use `lapply` on the resulting list object with `table` to return cell counts of each class. For area of each raster class, you just use a standard conversion of cell area and counts. Here is a quick example.

First lets create a discrete raster with values [5:8] and two polygons.

``````library(raster)

r <- raster(ncol=36, nrow=18)
r[] <- round(runif(ncell(r),5,8),0)

polys <- spPolygons(rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20)),
rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)))

plot(r)
``````

Now we can extract the raster values for each polygon. The `lapply` function is used with `table` to get the cell counts of each raster class.

``````( v <- extract(r, polys) )
( v.counts <- lapply(v, FUN = table) )
``````

We can now merge the raster cell counts with the polygon data. The `do.call` function is used to create a data.frame from the list object. The list is ordered the same as the polygons so it can be directly merged. The columns are cell counts for each raster class value intersecting a given polygon.

``````polys <- SpatialPolygonsDataFrame(polys, data.frame(IDS=1:length(polys),
do.call("rbind", v.counts)))
names(polys@data)[2:ncol(polys)] <- paste("class", 5:8, sep=".")
polys@data
``````