Skip to content

Instantly share code, notes, and snippets.

@JoshOBrien
Forked from Pakillo/polygonizer.R
Last active February 24, 2016 17:29
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save JoshOBrien/2bf028356f38538b5ddf to your computer and use it in GitHub Desktop.
Save JoshOBrien/2bf028356f38538b5ddf to your computer and use it in GitHub Desktop.
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')
pypath <- tools::file_path_sans_ext(normalizePath(pypath))
}
## 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)
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment