[R-sig-Geo] R version of gdal_polygonize.py

[Lyndon Estes]

Hi Jonathan,

With rasterToPolygons, the dissolve = T argument should do what you
want in terms of merging common value polygons, provided you have
rgeos installed.

I also use the a system call to gdal_polygonize.py, which seems to
help for larger polygons.

runGdalPolygonize <- function(inraster, outshape, attname, readin =
TRUE, gdalformat = "ESRI Shapefile") {
# Runs gdal_polygonize.py to create polygons from raster (with limited
options at present)
# Args:
#   inraster: The name of a raster object referring to a permanent,
GDAL-readable raster
#   outshape: The name of the desired output polygon file
#   attname: The name for the column in the output attribute table
#   gdalformat: Defaults to ESRI shapefile at present
#   readin: Option specifies whether result is read into an R object
or not, defaults to TRUE. Set to FALSE
#     (and don't assign to an object) and it will just write to disk
(useful for when polygons are very large)
# Returns:
#   Polygon file, read back in to an R object
# Notes:
#   This needs to be run within the directory containing the raster file

  py.c <- "python2.7
/Library/Frameworks/GDAL.framework/Versions/1.8/Programs/gdal_polygonize.py"
  rast.nm <- unlist(strsplit(inraster at file@name,
"/"))[length(unlist(strsplit(inraster at file@name, "/")))]
  full.c <- paste(py.c, rast.nm, "-f", paste("'", gdalformat, "'", sep = ""),
                  paste(outshape, ".shp", sep = ""), outshape, attname)
  system(full.c, intern = TRUE)
  if(readin == TRUE) {
    shp <- readOGR(dsn = paste(outshape, ".shp", sep = ""), layer = outshape)
    return(shp)
  }
}

Cheers, Lyndon

On Tue, Feb 12, 2013 at 9:40 AM, Jonathan Greenberg <jgrn at illinois.edu> wrote:
> R-sig-geo'ers:
>
> Continuing on my quest to *precisely* reproject/resize a categorical raster
> where the output cell is the category which has the maximum areal cover of
> that cell (this is not the same as nearest neighbor), it is looking more
> and more like the solution is to first convert a class raster to a polygon
> surface, reproject/resize, then re-rasterize the output.
>
> As such, is there a function within R that does this?  I see raster has the
> "rasterToPolygons" but as far as I can tell, this produces one polygon per
> cell (e.g. it doesn't merge the cells if they contain the same values --
> Robert, correct me if I'm wrong here) -- this could cause some memory
> issues, I imagine, on some big files.  I'm looking for something like
> what gdal_polygonize.py
> does.
>
> Any thoughts?
>
> --j
>
> --
> Jonathan A. Greenberg, PhD
> Assistant Professor
> Global Environmental Analysis and Remote Sensing (GEARS) Laboratory
> Department of Geography and Geographic Information Science
> University of Illinois at Urbana-Champaign
> 607 South Mathews Avenue, MC 150
> Urbana, IL 61801
> Phone: 217-300-1924
> http://www.geog.illinois.edu/~jgrn/
> AIM: jgrn307, MSN: jgrn307 at hotmail.com, Gchat: jgrn307, Skype: jgrn3007
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo