Adding column with data extracted from raster stack with each row from different band depending on data frame characteristic using R

Question

I have a data frame with observations, each of them has, among others, coordinates where the observation was taken, and the year where it was taken.

I want to add a column with the temperature at these coordinates and year, and I have a raster file with the mean temperature of the region at each year, therefore 40 rasters, per 40 years.

I am quite noob at working with R and GIS. For now, I would group the observations per year, and then extract the temperature per each coordinates of each year. But since I have over 40 years, I was wandering if there is a faster way to do this. Here a bit of code to have a better idea of the problem, to create this data frame:

            Species Year Longitud Latitude
1      S.pilchardus 2014        8       52
2     T. silvestris 2010        8       54
3       C. monogyna 2014        9       58
4       C. monogyna 2010        8       55
5     T. silvestris 2015        4       55
6  M.novaezelandiae 2014        5       51
7     T. silvestris 2011        5       52
8      S.pilchardus 2015        3       50
9       C. monogyna 2014        1       60
10 M.novaezelandiae 2013        4       58    

species = c("C. monogyna", "T. silvestris", 
           "M.novaezelandiae", "S.pilchardus")

database = data.frame(Species = sample(species, 10, replace = T),
  Year = sample(2010:2015,10, replace = T), 
  Longitud = sample(1:10, 10, replace = T), 
  Latitude = sample(50:60, 10, replace = T))

raster = stack(temperature_2010_2015) 

Assuming that raster is a stack group of 6 rasters, one per each year from 2010 to 2015, how could I add a column to the "database" object, that contains the value of temperature for the "Longitude, Latitude" coordinates in the year of the row?

Answer 1

This is a fairly common question on this forum and, in the future, I would encourage you to spend some time with search terms exploring previous posts. I am addressing this, rather than flagging it as a duplicate and providing a link to a previous answer, because the terra package is functionally starting to replace the raster package and this is a good opportunity to start demonstrating terra workflows.

First, you need to add packages that support spatial classes in R.

library(terra)
library(sf)

We can now create your data.frame and coerce it into a sf POINT object.

database = data.frame(Species = sample(c("C. monogyna", "T. silvestris", 
  "M.novaezelandiae", "S.pilchardus"), 10, replace = T),
  Year = sample(2010:2015,10, replace = T), 
  Longitude = sample(1:10, 10, replace = T), 
  Latitude = sample(50:60, 10, replace = T))

( spp <- sf::st_as_sf(database, coords = c("Longitude", "Latitude"), 
              crs = "EPSG:4326", agr = "constant") )

Since I obviously do not have access to your climate data we can create some random rasters that will act as climate data (big thank you to Robert Hijmans for adding the nlyrs argument to rast) and plot the points on one of the rasters. Please note that I am adding names that reflect a set of years (dates) so that the data can be queried later.

clim <- terra::rast(terra::ext(spp), nrow=500, ncol=1000, nlyrs=20)
  clim[] <- runif(250000, 0, 32)

plot(clim[[1]])
  plot(sf::st_geometry(spp), pch=20, add=TRUE) 

Now, for extracting all of the climate data, it is as simple as using the terra::extract function. One big difference between raster and terra is that sp and sf vector class objects need to be coerced into a terra vector object using terra::vect.

( spp <- cbind(spp, terra::extract(clim, terra::vect(spp))) )

For conditional extraction say, by year, we need to iterate. For simplicity sake we will use a for loop but, one could write a function and pass it to lapply to act as the iterator.

This may look a bit nasty but, we are creating a vector that is the same length as our point data, iterating through the Years column and subsetting the points and climate raster by the i year. Both which and grep return an index (position in the data.frame or vector) based on a query. So, grep("2011", names(clim)) would return us the index of the 2011 raster which, we then index using a double bracket.

    climate.year <- rep(NA, nrow(spp))
      for(i in unique(spp$Year)) {
        y <- which(spp$Year == i)  
        climate.year[y] <- as.numeric(terra::extract(clim[[grep(i, 
          names(clim))]], vect(spp[y,]))[,2])
      }

This vector can now be added to the spp data spp$clim <- climate.year

Answer 2

Building on @Jeffery Evans' example data:

library(terra)
#terra version 1.3.2
set.seed(20210615)

database = data.frame(Species = sample(c("C. monogyna", "T. silvestris", 
  "M.novaezelandiae", "S.pilchardus"), 10, replace = T),
  Year = sample(2010:2015,10, replace = T), 
  Longitude = sample(1:10, 10, replace = T), 
  Latitude = sample(50:60, 10, replace = T))

# create a SpatVector
spp <- vect(database, c("Longitude", "Latitude"), crs = "EPSG:4326")

# example raster data
clim <- terra::rast(ext(spp)+2, nrows=50, ncols=100, nlyrs=15)
values(clim) <- runif(prod(dim(clim)), 0, 32)
names(clim) <- 2005:2019

I agree with Jeffrey that this is a fairly common problem, and this has now been addressed in extract, with argument layer. This argument that can either be a name in the SpatVector, or a vector of layer-names or layer-numbers. With that, you can do this:

extract(clim, spp, layer="Year")
#   ID  Year     value
#1   1  2013 11.741345
#2   2  2015  8.828645
#3   3  2014 31.127033
#4   4  2013 12.325310
#5   5  2015 24.297580
#6   6  2014  7.626536
#7   7  2012 11.511305
#8   8  2012 26.142886
#9   9  2013 24.423164
#10 10  2010 14.945456

or like this

extract(clim, spp, layer=as.character(spp$Year))

or like this

i <- match(spp$Year, names(clim))
extract(clim, spp, layer=i)