Issues converting Raster to Polygons in R

Question

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.

Answer 1

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)

poly <- gdal_polygonizeR(r)
plot(poly)