polygonizer.R

polygonizer <- 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
  # 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')
  }
  ## The line below has been commented:
  # 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='.asc')})
    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
  #system2('python',
  system2('C:\\OSGeo4W64\\OSGeo4W.bat',
    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=!quietish)
    return(shp) 
  }
  return(NULL)
}

Hi Luis,

I haven't used this function for a few years so it may be outdated. Note it was integrated in the APfun package, which probably still works: https://rdrr.io/cran/APfun/man/APpolygonize.html. Check out also further changes in the original function @johnbaums: https://gist.github.com/johnbaums/26e8091f082f2b3dd279

I'd probably try using terra::as.polygons, as shown here: https://stackoverflow.com/a/65335878. It can work with very large rasters and should be fast.

Other alternatives:

• Using sf and stars: https://gis.stackexchange.com/a/357792 & https://gis.stackexchange.com/a/313550
• spex package: https://mdsumner.github.io/spex/reference/polygonize.html
• fasterVectorize: https://rdrr.io/github/adamlilith/fasterRaster/man/fasterVectorize.html

Hope this helps

@luisrua I posted an alternative approach, using sf: https://gist.github.com/johnbaums/26e8091f082f2b3dd279#gistcomment-3142562

Dear @johnbaums and @Pakillo.
Thanks so much for your help. I will try the different options proposed.

Regards

Luis