# How to speed up raster to polygon conversion in R?

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:

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 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... • 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. May 8, 2016 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 ! ;) May 9, 2016 at 18:06 • Which step is taking the most time? Converting to polygons, or the intersection? May 23, 2016 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. May 23, 2016 at 4:55 ## 5 Answers 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). ### 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))
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. May 23, 2016 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... May 23, 2016 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 ! May 23, 2016 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. May 23, 2016 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 !! ;) May 27, 2016 at 18:19 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',
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)

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

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.