Skip to content

Instantly share code, notes, and snippets.

@johnbaums
Last active February 15, 2021 02:18
Show Gist options
  • Save johnbaums/26e8091f082f2b3dd279 to your computer and use it in GitHub Desktop.
Save johnbaums/26e8091f082f2b3dd279 to your computer and use it in GitHub Desktop.
Convert raster data to a ESRI polygon shapefile and (optionally) a SpatialPolygonsDataFrame
polygonizer <- function(x, outshape=NULL, pypath=NULL, readpoly=TRUE,
fillholes=FALSE, aggregate=FALSE,
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)
# pypath: the path to gdal_polygonize.py or OSGeo4W.bat (if NULL, the function
# will attempt to determine the location)
# readpoly: should the polygon shapefile be read back into R, and returned by
# this function? (logical)
# fillholes: should holes be deleted (i.e., their area added to the containing
# polygon)
# aggregate: should polygons be aggregated by their associated raster value?
# quietish: should (some) messages be suppressed? (logical)
if (isTRUE(readpoly) || isTRUE(fillholes)) require(rgdal)
if (isTRUE(aggregate)) require(rgeos)
if (is.null(pypath)) {
cmd <- Sys.which('OSGeo4W.bat')
pypath <- 'gdal_polygonize'
if(cmd=='') {
cmd <- 'python'
pypath <- Sys.which('gdal_polygonize.py')
if (!file.exists(pypath))
stop("Could not find gdal_polygonize.py or OSGeo4W on your system.")
}
}
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(cmd, args=(
sprintf('"%s" "%s" %s -f "ESRI Shapefile" "%s.shp"',
pypath, rastpath, ifelse(quietish, '-q ', ''), outshape)))
if(isTRUE(aggregate)||isTRUE(readpoly)||isTRUE(fillholes)) {
shp <- readOGR(dirname(outshape), layer=basename(outshape),
verbose=!quietish)
} else return(NULL)
if (isTRUE(fillholes)) {
poly_noholes <- lapply(shp@polygons, function(x) {
Filter(function(p) p@ringDir==1, x@Polygons)[[1]]
})
pp <- SpatialPolygons(mapply(function(x, id) {
list(Polygons(list(x), ID=id))
}, poly_noholes, row.names(shp)), proj4string=CRS(proj4string(shp)))
shp <- SpatialPolygonsDataFrame(pp, shp@data)
if(isTRUE(aggregate)) shp <- aggregate(shp, names(shp))
writeOGR(shp, dirname(outshape), basename(outshape),
'ESRI Shapefile', overwrite=TRUE)
}
if(isTRUE(aggregate) & !isTRUE(fillholes)) {
shp <- aggregate(shp, names(shp))
writeOGR(shp, dirname(outshape), basename(outshape),
'ESRI Shapefile', overwrite=TRUE)
}
ifelse(isTRUE(readpoly), return(shp), return(NULL))
}
@johnbaums
Copy link
Author

@rafaleo you can use sf now:

polygonize <- function(srcfile, dstfile, mask, ogr_format, fieldname, band, connect8=FALSE) {
  options <- if(isTRUE(connect8)) 'CONNECT8=8' else character(0)
  contour_options <- character(0)
  if(missing(mask)) mask <- character(0)
  .Call("_sf_CPL_polygonize", PACKAGE = "sf", srcfile, mask, 
        'GTiff', ogr_format, layer=dstfile, options, 0, #iPixValField, 
        contour_options, use_contours=FALSE, use_integer=TRUE)
}
  
library(sf)
library(raster)
library(gdalUtilities) # for example raster
r <- raster::raster(system.file(package='gdalUtilities', 'extdata', 'tahoe.tif')) > 250
raster::writeRaster(r, f <- tempfile(fileext='.tif'))
polygonize(f, f2 <- tempfile(fileext='.shp'), ogr_format='ESRI Shapefile', connect8=TRUE)
plot(sf::read_sf(f2)) 

@rafaleo
Copy link

rafaleo commented Jan 18, 2020

@johnbaums Thank you for those alternatives. Thankfully I've managed to repair my python and can use that code as well.

@lovalery
Copy link

@johnbaums Dear John, I would like to polygonize a list of rasters using the lapply function but the polygonizer function seems to accept only one raster object (if I'm not wrong).
So, I was wondering if it would be possible for you to slightly modify the script of the function so that a list of rasters can be processed.
Thank you in advance for your answer.

@johnbaums
Copy link
Author

johnbaums commented Jun 30, 2020

@lovalery I don't think you should have any issues using lapply on a list of rasters:

library(raster)
L <- replicate(3, raster(matrix(runif(100), nrow=10)))
polys <- lapply(L, polygonizer)

See also this comment for an approach that uses sf (with a self-contained version of GDAL/OGR).

@lovalery
Copy link

lovalery commented Jul 1, 2020

@johnbaums Thank you very much for your reply.
lapply doesn't work in my case (I'm not sure why because, as you say, it should work).
I've tried the approach using 'sf' object and it works. That's the main thing :-)
Thanks again for your advice.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment