1

I have some problems to convert a shapefile into a raster. I used the function rasterize (package raster) but all raster values are equal to NA.

Following the answers to this question rasterize produces a raster full of NA values. The raster aligns with the SpatialPolygons:

r.ext <- raster(shp.prj, res=30)
rfn <- rasterize(shp.prj, r.ext, field="land_cover_type")

class(shp.prj$land_cover_type)
# [1] "character"

But I have always NAs in the output raster (rfn).

> rfn
class       : RasterLayer 
dimensions  : 1, 1, 1  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : -74, -44, 15.875, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
data source : in memory
names       : layer 
values      : NA, NA  (min, max)

Here are some information on the input shapefile (shp.prj) and the input raster (r.ext):

> shp.prj
class       : SpatialPolygonsDataFrame 
features    : 78152 
extent      : -74, -71.5, 45, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0 
variables   : 3

> r.ext
class       : RasterLayer 
dimensions  : 1, 1, 1  (nrow, ncol, ncell)
resolution  : 30, 30  (x, y)
extent      : -74, -44, 15.875, 45.875  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0

I also tried to reduce the raster resolution (res(r.ext) <- 0.05) but I obtained this error message:

> rfn <- rasterize(shp.prj, r.ext, field="land_cover_type")
Error in cbind(poly, rbind(poly[-1, ], poly[1, ])) : 
  trying to get slot "coords" from an object of a basic class ("NULL") with no slots

Following the answers to this question r raster package rasterize failing:no slots, I also tested:

rfn <- rasterize(shp.prj, r.ext, shp.prj$land_cover_type)

but I have always the error message:

Error in cbind(poly, rbind(poly[-1, ], poly[1, ])) : 
  trying to get slot "coords" from an object of a basic class ("NULL") with no slots
  • r.ext does not have the the same extent as r.ext <- raster(extent(shp.prj)) - so that presumably is not the actual code involved here. (You could just do raster(shp.prj) to avoid the need for projection()<-). Otherwise I would try rasterize(shp.prj, r.ext, field=shp.prj$land_cover_type) – mdsumner Feb 4 '16 at 5:11
  • Thank you very much midsummer. I have NAs in both the input raster r.ext and the output raster rfn (I added the rfn output in the post). I also tried to do r.ext <- raster(shp.prj) to avoid projection() and rasterize(shp.prj, r.ext, field=shp.prj$land_cover_type) but I have always NAs in the output raster rfn with the raster resolution = 30. With the raster resolution = 0.05, I have always the error message with "no slots". – Nell Feb 4 '16 at 16:10
3

Since the data is in a geographic projection, you cannot assign 30 to the resolution. If you look at the result of doing so, the r.ext raster has 1 row, 1 column and 1 value. So, that problem solved.

Try defining the resolution in the call to raster. Since sometimes nodata begets nodata, I am assigning a value of 1 to all cells in the reference raster "r.ext". You will see that the resulting raster now has 18 rows, 50 columns and 900 values.

library(raster)
r.ext <- raster(extent(-74, -71.5, 45, 45.875), resolution = 0.05)
  r.ext[] <- rep(1,ncell(r.ext))
  projection(r.ext) <- CRS("+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0")

r.ext  
plot(r.ext)

Now (this will make @mdsumner's head explode) try a direct call to the @data slot to access the "land_cover_type" column. I have found this necessary in the past. These packages are ever evolving and what was previously necessary may not be in the next release. However, I have never had success simply passing the column name to the argument (ie., field="land_cover_type"), it always has to be a call the the object with the name (eg., x$land_cover_type or x@data$land_cover_type). I have found that sometimes coercing to a SpatialPolygons object sidesteps issues.

rfn <- rasterize(as(shp.prj, "SpatialPolygons"), r.ext, shp.prj@data[,"land_cover_type"], fun="first")  

Following the example from the functions help, here is a worked example.

library(raster)
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
pols <- spPolygons(p1, p2, p3)
pols <- SpatialPolygonsDataFrame(pols, data.frame(ID=1:3))
ref <- raster(ncol=90, nrow=45)

r <- rasterize(pols, ref, field = pols@data[,"ID"], fun="first")
plot(r)  

I would highly recommend selecting a distance based projection, reprojecting your polygon data and then running this analysis. Then you can define a cell size of 30m.

This is not a failure of the rasterize function but rather a failing of the analysis/code calling the function. Developers work very hard, and for free, to provide these packages. Please do not intimate a bug or failure until your are sure that the fault is not yours.

  • Thank you very much Jeffrey for your answer. I re-projected my polygon data by using a projected coordinate system as follows: shp.prj <- spTransform(shp,CRS("+proj=....")). Then, I re-built my input raster as follows: r.ext <- raster(extent(sief.prj), resolution = 30). Indeed, my raster has now several rows, columns and values. However, I still have a problem with rasterize even if I do shp.prj@data[,"land_cover_type"]. I obtain Error in cbind(poly, rbind(poly[-1, ], poly[1, ])): trying to get slot "coords" from an object of a basic class ("NULL") with no slots . – Nell Feb 4 '16 at 18:19
  • Sorry, I forgot that the expected input to the rasterize function is a SpatialPolygons object. Try my modification to the rasterize syntax and let me know if it works. – Jeffrey Evans Feb 4 '16 at 20:22
  • Thank you very much Jeffrey for your help. It works when the field "Land_cover_type" is a factor. When the field is a character, I still have the error message: Error in cbind(poly, rbind(poly[-1, ], poly[1, ])): trying to get slot "coords" from an object of a basic class ("NULL") with no slots. – Nell Feb 5 '16 at 16:04
  • Since a raster has to be numeric I imagine that this is because it has to coerce the column to numeric before conversion. The underlying code is probably assuming a factor when the column is non-numeric. This is good to know. – Jeffrey Evans Feb 5 '16 at 17:02
2

Please show your sessionInfo() and traceback() after the error occurs. This error message: Error in cbind(poly, rbind(poly[-1, ], poly[1, ])): trying to get slot "coords" from an object of a basic class ("NULL") with no slots is very odd as it would normally only occur when expecting SpatialPoints.

As it seems you have very many, presumably very small, polygons, Consider using rasterize(..., small=TRUE).

Here to illustrate how it should work:

#example data
p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
hole <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))
p1 <- list(p1, hole)
p2 <- rbind(c(-10,0), c(140,60), c(160,0), c(140,-55), c(-10,0))
p3 <- rbind(c(-125,0), c(0,60), c(40,5), c(15,-45), c(-125,0))
p <- spPolygons(p1, p2, p3)
p <- SpatialPolygonsDataFrame(p, data.frame(id = 5:7, 
                            landcover=c('forest','savanna','forest')))

# now rasterize

r <- raster(p, res=2.5)
r <- rasterize(p, r, field='landcover')

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.