Simple tools to use gdal to polygonize and rasterize

rasterTool.R

library(cdlTools)
library(rgdal)
library(raster)
library(rgeos)


# rasterize in the R 'raster' package is really slow, depending on your needs this may work better
# Author: Jonathan Lisic
# License BSD
polyExtract <- function(x) {

  tempTif <- tempfile(fileext=".tif")
  tempShp <- tempfile(fileext=".shp")

  writeRaster(x,tempTif,overwrite=TRUE)

  # use gdalPolygonize
  system( sprintf( "gdal_polygonize.py %s -b 1 -f 'ESRI Shapefile' %s", tempTif, tempShp))

  r.out <- readOGR(dirname(tempShp), gsub("[.]shp","",basename(tempShp)))

  # free up the data
  unlink(tempTif)
  unlink(tempShp)

  return(r.out)
}


# note that since we delete the temp file, r.out must be in memory only
rasterExtract <- function(x,r,variable) {

  tempTif <- tempfile(fileext=".tif")
  tempShp <- tempfile(fileext=".shp")

  # write polygon
  writeOGR(x,dirname(tempShp),gsub("[.]shp","", basename(tempShp)),driver="ESRI Shapefile",overwrite=TRUE)

  # use gdalPolygonize
  if( missing(variable) ) {
    system( 
        sprintf( "gdal_rasterize  -burn 1 -ot UInt32 -a_nodata 0 -a_srs '%s' -tr %f %f -te %f %f %f %f %s %s", 
                  projection(r), xres(r), yres(r), xmin(r), ymin(r), xmax(r), ymax(r), 
                  tempShp, tempTif) 
         )
  } else {
    system( 
        sprintf( "gdal_rasterize  -a '%s' -ot UInt32 -a_nodata 0 -a_srs '%s' -tr %f %f -te %f %f %f %f %s %s", 
                  variable, projection(r), xres(r), yres(r), xmin(r), ymin(r), xmax(r), ymax(r), 
                  tempShp, tempTif) 
         )
  }

  # create a new raster object
  r.out <- raster(extent(r), ncols=ncol(r), nrows=nrow(r), crs=projection(r))
  values(r.out) <- values(raster(tempTif))

  # free up the data
  unlink(tempTif)
  unlink(tempShp)

  return(r.out)
}

Thanks for this, it helped me a lot.