How to speed up raster to polygon conversion in R?

Question

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

Answer 1

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

Answer 2

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!

Answer 3

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)
  plot(polys, add=TRUE)

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

Answer 4

from : https://sites.google.com/site/eospansite/installinggdalonmacosx

---

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

Installing GDAL on Mac OSX and Polygonize Rasters

• Go to http://www.kyngchaos.com/software/frameworks and install the Unix Compatibility Frameworks in this order:
  • GDAL Complete
  • GSL framework
  • FreeType
  • cairo

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,

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

---

---

READY TO GO TO R! :)

Read libraries:

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 (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)
}

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

Answer 5

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.