Convert raster to polygon using GDAL

polygonizer.R

# Convert raster to polygon using GDAL
## edited 2019-03-27 to use sf::st_read
polygonise <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile', 
                        pypath=NULL, readpoly=TRUE, quietish=TRUE) {
  ## x: an R Raster layer, or the file path to a raster file recognised by GDAL
  ## outshape: the path to the output shapefile (if NULL, a temporary file will be created)
  ## gdalformat: the desired OGR vector format. Currently only supports ESRI Shapefile!
  ## pypath: the path to gdal_polygonize.py (if NULL, an attempt will be made to determine the location
  ## readpoly: should the polygon shapefile be read back into R, and returned by this function? (logical)
  ## quietish: should (some) messages be suppressed? (logical)
  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.")
  }
  ## Now 'python' has to be substituted by OSGeo4W
  ## hardcoded to default OSGeo4W install location
  system2("C:/OSGeo4W64/OSGeo4W.bat",
    args = (sprintf(
      '"%1$s" "%2$s" -f "%3$s" "%4$s.shp"',
      pypath, rastpath, gdalformat, outshape
    ))
  )
  if (isTRUE(readpoly)) {
    require(sf)
    shp <- st_read(paste0(outshape, ".shp"), quiet = quietish)
    return(shp)
  }
  return(NULL)
  }