How to extract the mean raster value for a polygon based on matching a layer from a raster stack in R

Question

Apologies for my naivety (and general incompetence), first time asking a question here and I'm not particularly good at R. I have tried to provide a reproducible example below.

I have a stack of 365 rasters with each raster containing temperature values for a day of the year. For ease I have simplified everything in this example:

library(raster)
library(sf)

# Create raster brick
Jan <- raster(nrows = 50, ncols = 50, res = 0.5, xmn = -10, xmx = 10, ymn = -10, 
             ymx = 10, vals = 0.5)
             
Feb <- raster(nrows = 50, ncols = 50, res = 0.5, xmn = -10, xmx = 10, ymn = -10, 
              ymx = 10, vals = 1)
r_brick <- brick(Jan,Feb)
names(r_brick) <- c("Jan","Feb")

I also have a set of polygons which are of the class SF. Each polygon has a date:

# Create polygons
x = rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))
m = matrix(c(3, 0), 5, 2, byrow = TRUE)
z = x + m

p1 = Polygons(list( Polygon(x[5:1,])), "Jan")
p2 = Polygons(list( Polygon(z[1:5,])), "Feb")
r = SpatialPolygons(list(p1,p2))
a = (st_as_sf(r))

I want to get the mean raster values for each polygon from the rasterlayer within the brick that has the corresponding date (or month in this case) to that polygon.

I have successfully managed to do this by extracting raster data using the centroid of the polygon:

# Add a layer index by matching the month to the month of the raster
a$layerindex = match(a$month, names(r_brick))

# Add x and y informaiton of the centroid 
centroid <- st_centroid(a)
plot(centroid, add=T)
centroid <- st_coordinates(centroid)
a <- cbind(a,centroid)

# Take a sample using the central coordinate
a$sample <- sapply(1:nrow(a), function(i){
  raster::extract(r_brick, cbind(a$X[i],a$Y[i]), 
                  layer=a$layerindex[i], nl=1, na.rm=T)})

but I can't seem to write it correctly to take the mean for the entire layer. I tried:

# Add a layer index by matching the month to the month of the raster
a$layerindex = match(a$month, names(r_brick))

# Add x and y informaiton of the centroid 
centroid <- st_centroid(a)
plot(centroid, add=T)
centroid <- st_coordinates(centroid)
a <- cbind(a,centroid)

# Take a sample using the central coordinate
a$sample <- sapply(1:nrow(a), function(i){
  raster::extract(r_brick, a$geometry[i], 
                  layer=a$layerindex[i], nl=1, fun = mean, na.rm=T)})

but just get the following error.

Error: unexpected '}' in:
"  raster::extract(r_brick, a$geometry[i], 
                  layer=a$layerindex[i], nl=1, na.rm=T, fun=mean)}"

If anyone has any suggestions then I'd be very grateful. In an ideal world I would use exactextractr because I've found it to be much faster than using raster::extract so if anyone has any suggestions with that it would be even better but just working it out using raster::extract would be great for now.

Thanks in advance!

Robert Hijmans · Accepted Answer · 2023-02-03T03:12:28.717

"terra" has replaced "raster" so let's use that.

Example data

library(terra)
r <- rast(nrows=50, ncols=50, ext=c(-10, 10, -10, 10), nlyr=2,
          vals=0.5, names=c("Jan","Feb")) * 1:2

# creating a SpatVector, but it also works with sf polygons
x <- rbind(c(-1,-1), c(1,-1), c(1,1), c(-1,1), c(-1,-1))
z <- x + matrix(c(3, 0), 5, 2, byrow = TRUE)
xz <- rbind(cbind(1, 1, x), cbind(2, 1, z))
p <- vect(xz, "polygons", data.frame(month=c("Feb", "Jan")))

Match each polygon to a layer and use that in extract

i <- match(p$month, names(r))

e <- extract(r, p, layer=i, fun=mean)
e
#  ID layer value
#1  1   Feb   1.0
#2  2   Jan   0.5

Also see these posts: #1, #2

How to extract the mean raster value for a polygon based on matching a layer from a raster stack in R

1 Answers1