-
-
Save cromanpa94/d287fc0fbc8d2982565271b6a985fa73 to your computer and use it in GitHub Desktop.
Convert raster data to a ESRI polygon shapefile and (optionally) a SpatialPolygonsDataFrame
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| 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 <- system("mdfind gdal_polygonize.py", intern=T)[1] ##A bit slower but it ends up finding it in mac | |
| 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)) | |
| } | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Sys.whichwas unable to locategdal_polygonize.pyunder mac. I'm now using a more explicit command:system("mdfind gdal_polygonize.py", intern=T). Note that this code will use the first path among the ones that are found.