Shapefile to Raster Conversion in R

Converting shapefile to raster

You assumed that the polygons were in lon/lat coordinates, but they are not:

library(raster)
library(rgdal)
p <- shapefile('Global_Threats/Integrated_Future/rf_int_2030_poly.shp')
p

#class       : SpatialPolygonsDataFrame 
#features    : 63628 
#extent      : -18663508, 14601492, -3365385, 3410115  (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=cea +lon_0=-160 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#variables   : 3
#names       :    ID, THREAT, THREAT_TXT 
#min values  :     1,      0,   Critical 
#max values  : 63628,   2000,  Very High

You can either change the projection

pgeo <- spTransform(p, CRS('+proj=longlat +datum=WGS84'))

and then do something like:

ext <- floor(extent(pgeo))
rr <- raster(ext, res=0.5)
rr <- rasterize(pgeo, rr, field=1)

Or keep the orginal CRS and do something like:

ext <- extent(p)
r <- raster(ext, res=50000)  
r <- rasterize(p, r, field=1)
plot(r)

Note that you are rasterizing very small polygons to large raster cells. A polygon is considered 'inside' if it covers the center of a cell (i.e. assuming a case where polygons cover multiple cells). So for these data you would need to use a much higher resolution (and then perhaps aggregate the results). Alternatively you could rasterize polygon centroids.

But none of the above is relevant really, as you are doing this all backwards. The polygons are clearly derived from a raster (look how blocky they are) and the raster is available in the dataset you point to!

So instead of rasterizing, do:

x <- raster('Global_Threats/Integrated_Future/rf_int_2030')
x
#class       : RasterLayer 
#dimensions  : 25456, 80150, 2040298400  (nrow, ncol, ncell)
#resolution  : 500, 500  (x, y)
#extent      : -20037508, 20037492, -6363885, 6364115  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 
#data source : C:\temp\Global_Threats\Integrated_Future\rf_int_2030 
#names       : rf_int_2030 
#values      : 0, 2000  (min, max)
#attributes  :
#   ID  COUNT THREAT_TXT
#    0  80971        Low
#  100 343535     Medium
# 1000 322231       High
# 1500 168518  Very High
# 2000  83598   Critical

Here plotting a part of Palawan:

e <- extent(c(-8990636, -8929268, 1182946, 1256938))
plot(x, ext=e)
plot(p, add=TRUE)

If you need a lower resolution see raster::aggregate. For a different coordinate reference system, see raster::projectRaster.

Shapefile to raster conversion in R?

You are right to think that you should be using raster (rather than the sp raster Spatial classes) for spatial raster data. You should also use rgdal (rather than maptools) for reading, writing, and otherwise manipulating spatial vector data.

This should get you started:

library(rgdal)
library(raster)

## Read in the ecoregion shapefile (located in R's current working directory)
teow <- readOGR(dsn = "official_teow/official", layer = "wwf_terr_ecos")

## Set up a raster "template" to use in rasterize()
ext <-  extent (-95, -50, 24, 63)
xy <- abs(apply(as.matrix(bbox(ext)), 1, diff))
n <- 5
r <- raster(ext, ncol=xy[1]*n, nrow=xy[2]*n)

## Rasterize the shapefile
rr <-rasterize(teow, r)

## A couple of outputs
writeRaster(rr, "teow.asc")
plot(rr)

Sample Image

convert shapefile of multiple polygons to raster in R

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.

shapefile to a raster .tif file R

Here is a minimal, reproducible, self-contained example

library(raster)
p <- shapefile(system.file("external/lux.shp", package="raster"))
r <- raster(p, res=0.01)
p$land <- ifelse(p$ID_2 > 6, "Cultivated Land", "Unused Land")

The error occurs because you are trying to rasterize a character variable. That is not supported. You can use a factor, though

p$land <- as.factor(p$land)
r <- rasterize(p, r, "land")
r
#class      : RasterLayer 
#dimensions : 73, 78, 5694  (nrow, ncol, ncell)
#resolution : 0.01, 0.01  (x, y)
#extent     : 5.74414, 6.52414, 49.45162, 50.18162  (xmin, xmax, ymin, ymax)
#crs        : +proj=longlat +datum=WGS84 +no_defs 
#source     : memory
#names      : layer 
#values     : 1, 2  (min, max)

How to convert shapefile into vector geometry in R so it can be plotted?

It used to be necessary to convert shapefiles to polygons to make them usable for ggplot; the method was ggplot2::fortify(). This method is still valid (in meaning it was not descoped from current ggplot version).

But while this was necessary (note the past tense) and while you can still find a lot of oldish materials describing the process it is no longer required. Nor is it current best practice.

Since ggplot 3.0.0 (released in July 2018) ggplot has support for sf, and the best practice for visualizing spatial data ("shapefiles") in ggplot is ggplot2::geom_sf() which requires no prior conversion.

The best part is that shapefiles in context of {sf} package are also data frames - meaning that you can use standard data frame manipulating tools (such as dplyr::left_join()) to join a dataframe with data information (e.g. data per some administrative units) to a dataframe with spatial information (e.g. polygons of the administrative units)

How to rasterize a vector for area variable?

In cases like that you should rasterize density instead of area

Example data:

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
# all areas have 100 ha of the crop
v$crop_area <- 100
r <- rast(v, res=.01)

Compute density

e <- expanse(v, unit="ha") # or "km" to avoid small numbers
v$density <- v$crop_area / e

Rasterize

x <- rasterize(v, r, "density")

Back to area

ra <- cellSize(r, unit="ha")         
area <- x * ra

Check that the numbers are reasonable (error should be smallest for large areas / small cells). The expected value for each polygon is 100.

extract(area, v, sum, na.rm=TRUE, exact=TRUE) |> round()
#      ID density
# [1,]  1      99
# [2,]  2     101
# [3,]  3      99
# [4,]  4      95
# [5,]  5      99
# [6,]  6      98
# [7,]  7      96
# [8,]  8      99
# [9,]  9      98
#[10,] 10      98
#[11,] 11     101
#[12,] 12     100

Below I show how to do this in a loop over multiple years