1

Having some difficulties converting a generated Raster to a PolygonShapeFile. In R, one can use rasterToPolygons to convert a raster set to a set of polygons. Usually, this takes care of everything I need, but I'm converting a very large raster to a set of Polygons. Following the success of this user's function in dealing with large rasters, I tried to use gdal_polygonizeR to create a polygon from my raster. However, it crashes with the following error message when I try to construct the raster:

 Error in readOGR(dirname(outshape), layer = basename(outshape), verbose = !quiet) : 
 no features found In addition: Warning messages:
1: In .local(x, filename, ...) : all cell values are NA
2: In ogrFIDs(dsn = dsn, layer = layer) : no features found

Below I provide two code examples, the first which produces the correct object with R's native rasterToPolygons, and the second with gdal_polygonizeR which fails with the error message above.

library(tigris)
library(rgdal)
library(rgeos)
library(maptools)

states_sp = counties(state = 'RI', cb = TRUE)
states_sp2 = spTransform(states_sp, CRS("+init=epsg:2163 +units=ft"))
grid = raster(extent(states_sp2))
res(grid) = 2640
proj4string(grid) = "+init=epsg:2163 +units=ft"
gridpolygon = rasterToPolygons(grid)
states_sp2.grid = intersect(states_sp2, gridpolygon)
states_sp3.grid = spTransform(states_sp2.grid, 
    CRS("+proj=longlat +no_defs +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0"))
plot(states_sp3.grid)

The second example, below:

library(tigris)
library(rgdal)
library(rgeos)
library(maptools)

## Define the function
gdal_polygonizeR <- function(x, outshape=NULL, gdalformat = 'ESRI Shapefile',
                             pypath=NULL, readpoly=TRUE, quiet=TRUE) {
  if (isTRUE(readpoly)) require(rgdal)
  if (is.null(pypath)) {
    pypath <- Sys.which('gdal_polygonize.py')
  }
  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='.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.')
  print(outshape)
  system2('python', 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=!quiet)
    return(shp)
  }
  return(NULL)
}

states_sp = counties(state = 'RI', cb = TRUE)
states_sp2 = spTransform(states_sp, CRS("+init=epsg:2163 +units=ft"))
grid = raster(extent(states_sp2))
res(grid) = 2640
proj4string(grid) = "+init=epsg:2163 +units=ft"
gridpolygon = gdal_polygonizeR(grid)
states_sp2.grid = intersect(states_sp2, gridpolygon)
states_sp3.grid = spTransform(states_sp2.grid, 
                              CRS("+proj=longlat +no_defs +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0"))
plot(states_sp3.grid)

With the help of the generous answerer, some workarounds were created. In particular, I used

states_sp = counties(state = 'RI', cb = TRUE)
states_sp2 = spTransform(states_sp, CRS("+init=epsg:2163 +units=ft"))
grid = raster(extent(states_sp2))
res(grid) = 2640
proj4string(grid) = "+init=epsg:2163 +units=ft"
# Fills the raster with values to create the grids
vals = 1:ncell(grid)
grid = setValues(grid, vals)
gridpolygon = gdal_polygonizeR(grid)
states_sp2.grid = intersect(states_sp2, gridpolygon)
states_sp3.grid = spTransform(states_sp2.grid, 
                              CRS("+proj=longlat +no_defs +datum=NAD83 +ellps=GRS80 +towgs84=0,0,0"))
plot(states_sp3.grid)

but the other solutions were equally valid.

2

You give as example a dataset with no data. I've used a lot of times gdal_polygonizeR with excellent results (it's really fast), even converting cloud masks of entire Landsat scenes to polygons.

Checking your example, I got this:

!is.na(grid)

class       : RasterLayer 
dimensions  : 123, 88, 10824  (nrow, ncol, ncell)
resolution  : 2640, 2640  (x, y)
extent      : 7459371, 7691691, -66019.53, 258700.5  (xmin, xmax, ymin, ymax)
coord. ref. : +init=epsg:2163 +units=ft +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +no_defs 
data source : in memory
names       : layer 
values      : 0, 0  (min, max)

You can't convert a data without values to polygon. Try this example and check the results:

r <- raster()
r <- setValues(r, values = rep(x=1:2, each=length(r)/2, length.out=length(r)))
plot(r)

raster

poly <- gdal_polygonizeR(r)
plot(poly)

plot

  • "You can't convert a data without values to polygon" - Can you explain why this statement doesn't hold for the rasterToPolygons function? Why is gdal_polygonizeR unable to accomplish the same? – genap Nov 13 '16 at 3:20
  • gdal_polygonizeR converts raster values to polygons. If two cells with the same values are neighbor, you'll get one polygon. Is the best way to polygonize categorical rasters instead two. rasterToPolygons converts each cell to polygon. Check help, if you have some cells with NAvalues, those cells aren't converted... But in this case your entire raster has no values, so, all cells are converted to polygons without data.frame. Check: gridpolygon = rasterToPolygons(grid) and gridpolygon@data – aldo_tapia Nov 13 '16 at 3:39
  • Okay. So, to approach my problem then, I should set each of my cells a unique color, which will be picked up by gdal_polygonizeR as a unique polygon. Is there any way to iterate through the cells of a raster and set the colors? – genap Nov 13 '16 at 3:43
  • The simplest syntax is poly <- as(extent(grid),"SpatialPolygons"), you'll get the extent as a polygon. Before this, you can add information: poly <- SpatialPolygonsDataFrame(poly,data=data.frame(ID=1,Obs="From raster")). Also, this is the faster way. – aldo_tapia Nov 13 '16 at 3:51
  • 1
    Alright. I used the following vals = 1:ncell(grid) and grid = setValues(grid, vals) to fill the raster, and gdal_polygonizeR was able to accept it correctly. Thank you your help! – genap Nov 13 '16 at 3:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.