convert shapefile of multiple polygons to raster in R

Question

I am in a big trouble in converting polygons to raster in R. What I wanted to do is: I have shapefile (i.e. polygons) of 574 species. That is in attribute table it has 574 rows (i.e., FID is between 0 to 573). A subset of data can be found here: https://drive.google.com/file/d/1AdTChjerCXopE1-PZIAPp5ZXISdq45i8/view?usp=sharing

I wanted convert it to raster. In the output raster I see that the minimum and maximum values are 1 and 574. What I suspect is: it is getting field ID in the cell as pixel value that should not. Cell values should come from covering polygons. Any help will be highly appreciated. Below is the sample code:

library(raster)
library(rgdal)
library(maptools)


# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

# load species polygons for all 574 species
# already in the home directory

sp <- readShapePoly("AllP574sp.shp", proj4string = projection1)

# load a raster file to use as a mirror for rasterize
raster1 <- raster("/data/projects/MeanTemp2050rcp45_BCC_CSM1_1.tif")


r.sp <- rasterize(sp, raster1) # rasterize our species polygon to the same resoluton of loaded raster

writeRaster(r.sp, "/data/projects/all574spRaster", format = "GTiff", overwrite = TRUE)

The properties of output raster is given below:

class       : RasterLayer 
dimensions  : 18000, 43200, 777600000  (nrow, ncol, ncell)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs 
data source : E:\all574spRaster.tif 
names       : all574spRaster 
values      : 1, 574  (min, max)

Answer 1

When you used the rasterize function, it is important to specify the field argument, otherwise by default it will try to create one for you; it looks like in your case it created one from the FID column.

I made some guesses to regenerate a working set of polygons that might be similar to yours.

library(maptools)
library(rgdal)
library(sp)
library(geosphere)

# set seed for duplicatable results
set.seed(1)

# some data that looks a little like yours
BINOMIAL <- c("Controversial chimneyswift", "Dull dungbeetle", 
"Easternmost eel", "Jumping jaeger", "Qualified queenconch")
FID <- 0:(length(BINOMIAL) - 1)
RANGE <-  runif(length(BINOMIAL), min = 118, max = 3875370)
MyData <- cbind.data.frame(FID, BINOMIAL, RANGE)
row.names(MyData) <- FID

# some semi-random polygons in your extent box
ext <- extent(c(-180, 180, -60, 90))
create_polygon <- function(n = 4, lat, lon, r) {
  lengths <- rnorm(n, r, r/3)
  smoother_lengths <- c(sort(lengths), rev(sort(lengths)))
  lengths <- smoother_lengths[sort(sample(n * 2, n))]
  lengths <- rep(lengths[1], length(lengths))
  directions <- sort(runif(n, 0, 360))
  p <- cbind(lon, lat) 
  vertices <- t(mapply(destPoint, b = directions, 
                       d = lengths, MoreArgs = list(p = p)))
  vertices <- rbind(vertices, vertices[1, ])
  sapply(vertices[,1], min, ext@xmax)
  sapply(vertices[,1], max, ext@xmin)
  sapply(vertices[,2], min, ext@ymax)
  sapply(vertices[,2], max, ext@ymin)
  Polygon(vertices) 
}
rand_lats <- runif(nrow(MyData), min = -50, max = 60)
rand_lons <- runif(nrow(MyData), min = -100, max = 100)
rand_sides <- sample(4:20, nrow(MyData), replace = TRUE)
rand_sizes <- rnorm(nrow(MyData), mean = 5e+06, sd = 1e+06)
make_species_polygon <- function(i) {
  p.i <- list(create_polygon(rand_sides[i], rand_lats[i], 
                             rand_lons[i], rand_sizes[i]))
  P.i <- Polygons(p.i, FID[i])
}

polys <- SpatialPolygons(lapply(1:nrow(MyData), make_species_polygon))
spdf <- SpatialPolygonsDataFrame(Sr = polys, data = MyData)

t.shp <- tempfile(pattern = "MyShapefile", fileext = ".shp")
raster::shapefile(spdf, t.shp)

At this point there is a shapefile written in your temp directory, whose name is stored in the variable t.shp. I intend that shapefile to be a workable duplicate of whatever big shapefile is your true one. So now we can look at your code, what it was doing, and what you would like it to do:

## now we get into your code
library(raster)
library(rgdal)
library(maptools)

# define porjection
projection1 <- CRS ("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

## 
## I don't know what your shapefile looks like exactly, 
## but substituting `t.shp` the tempfile that I created above
## also since the function readShapePoly is deprecated
## instead I use the recommended new function rgdal::readOGR() 
## 
sp <- rgdal::readOGR(t.shp)

##
## I don't know what your tiff file looks like exactly,
## but I can duplicate its characteristics
## for speed I have decreased resolution by a factor of 10
##
raster1 <- raster(nrow = 1800, ncol = 4320, ext)

# rasterize our species polygon to the same resoluton of loaded raster
r.sp <- rasterize(x = sp, y = raster1, field = MyData$RANGE) 

t.tif <- tempfile(pattern = "MyRastfile", fileext = ".tif")
writeRaster(r.sp, t.tif, format = "GTiff", overwrite = TRUE)

With the result as follows:

raster(t.tif)
class       : RasterLayer 
dimensions  : 1800, 4320, 7776000  (nrow, ncol, ncell)
resolution  : 0.08333333, 0.08333333  (x, y)
extent      : -180, 180, -60, 90  (xmin, xmax, ymin, ymax)
coord. ref. : NA 
data source : [["a file name in your temp directory"]] 
names       : MyRastfile1034368f3cec 
values      : 781686.3, 3519652  (min, max)

The result is now showing values that are taken from your RANGE column instead of your FID column.