Multi-layer rasters

ImportantSummary

This tutorial explore how to handle multi-layer rasters in R with terra package:

TipThe ecologist mind

How are otters influenced by climate? It would be ineresting to get the climate at the time of the observations of otter. To answer this question, we will need to discover another type of spatial data: multi-layer rasters.

Setup

Follow the setup instructions if you haven’t followed previous tutorials

If haven’t done it already, please follow the setup instructions.

Let’s start with loading the required packages.

suppressPackageStartupMessages({
  library(mapview)
  library(here)
  library(terra)
  # library(sf)
  # library(exactextractr)
})

Now load the observations of otters recorded in 2021 within a 50km buffer from Montpellier, France.

pt_otter <- vect(here("data", "gbif_otter_2021_mpl50km.gpkg"))
pt_otter <- vect(
  "https://github.com/FRBCesab/spatial-r/raw/main/data/gbif_otter_2021_mpl50km.gpkg"
)

Load monthly climate with CHELSA

We will use the monthly average temperature data from CHELSA which has a spatial resolution of 1km. To see how we created the data for our case study, have a look at this tutorial.

Import

The function to read most kind of raster file is terra::rast()

temperature <- rast(here("data", "CHELSA_monthly_tas_2015_2021.tif"))
url_github <- "https://github.com/FRBCesab/spatial-r/raw/main/data/"
temperature <- rast(paste0(url_github, "CHELSA_monthly_tas_2015_2021.tif"))
NoteYour turn
  1. What are the dimensions of the loaded temperature raster?
  2. What is the resolution of rasters? In which unit? Why?
  3. How can we get the name of the different layers?
Click to see the answer
  1. There are 116 rows, 165 columns and 84 layers in the raster. You can access it with dim(temperature) or just by typing temperature in the console.
  2. The spatial resolution of raster is 0.0083, 0.0083 decimal degrees (equal to 29.9999999 arc seconds). The resolution is expressed in degrees because the projection system is WGS 84 (EPSG 4326). You can access these information with res(temperature) and crs(temperature, describe = TRUE).
  3. To access the name of the layers, type names(temperature), they corresponds to monthly averaged temperature.

Rename

To simplify the next steps, we will rename the layers with only their corresponding month and year.

# shorten names
time_chelsa <- substr(names(temperature), 12, 18)
names(temperature) <- time_chelsa

Visualize

Let’s visualize the monthly temperature of 2021

# select the layers of 2021
in_2021 <- grep("2021", time_chelsa)
# plot the 2021 temperature layers
plot(temperature, in_2021)

NoteYour turn
  • what is the unit of the temperature values?
  • mask the sea from the temperature raster. addref to raster tutorial
Click to see the answer 1 The temperature are expressed in Kelvin, as indicated in the documentation of Chelsa-monthly dataset
Click to see the answer 2 
# get the border of the country (level = 0)
france_border <- geodata::gadm("FRA", level = 0, path = here("data"))

# mask (=set to NA) the pixels that are not in the polygon
temperature <- mask(temperature, france_border)

plot(temperature, in_2021)

Transform in degree Celsius

We can make algebra operations in rasters. For instance, we can transform Kelvin to degree celsius.

temp_C <- temperature - 273.15
# plot the 2021 temperature layers
plot(temp_C, in_2021)

Extract information on points

Same as with simple raster, the function terra::extract() will get the temperature values of all months at the location of the observations.

pt_temp <- extract(temp_C, pt_otter, ID = FALSE)
dim(pt_temp)
[1] 83 84
summary(unlist(pt_temp))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   0.95    9.55   13.75   14.69   20.55   26.95

The extracted temperature (in pt_temp) are stored in a data.frame with rows corresponding to the otter observations, and column corresponding to the month.

We can see the seasonal variation of the temperature with a simple boxplot()

# plot the seasonal variation of the mean temperature
boxplot(pt_temp, las = 2, ylab = "Temperature (*C)")

Average temperature

With all these monthly temperature values, we can calculate the average of monthly temperature for the period 2015-2021.

mean_temp <- apply(pt_temp, 1, mean, na.rm = TRUE)
# distribution of average values
summary(mean_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  10.74   14.04   14.74   14.69   15.76   16.31 
# attach the temperature as an attribute
pt_otter$mean_temp <- mean_temp

# let's map the average temperature
plot(pt_otter, "mean_temp")
plot(france_border, add = TRUE)

Seasonal variations

Another interesting indicators is the seasonality of monthly temperature measured as the coefficient of variation.

cv_temp <- apply(pt_temp, 1, sd, na.rm = TRUE) / mean_temp
summary(cv_temp)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.3334  0.3823  0.4326  0.4185  0.4445  0.5517 
pt_otter$cv_temp <- cv_temp
plot(pt_otter, "cv_temp")
plot(france_border, add = TRUE)

Observations made closer to the sea have warmer climate and lower seasonal variations in temperature.

Temperature at the time of observation

We need to identify which temperature layer correspond to the time of the observations. For this we will use the date formatting in R. Make sure to have a look at the documentations of as.Date() and strptime() if you want more details.

pt_otter$time <- as.Date(pt_otter$eventDate) |> format("%m_%Y")

table(pt_otter$time)

01_2021 02_2021 03_2021 04_2021 05_2021 06_2021 07_2021 08_2021 09_2021 10_2021 
     12      11      11       8       3       6       4       3       6       7 
11_2021 12_2021 
      5       7

Now that we have the same format for the name of the temperature layer and the time of the observations, we can get the temperature value at the time of observation.

# get temperature at the month of observation
t_obs <- match(pt_otter$time, names(pt_temp))
xy <- cbind(1:nrow(pt_otter), t_obs)
pt_otter$temp_obs <- pt_temp[xy]
summary(pt_otter$temp_obs)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   4.55    7.75   10.35   12.44   15.75   23.85 
plot(pt_otter, "temp_obs")
plot(france_border, add = TRUE)

NoteYour turn
  1. Extracting values for 1km buffer around the observation. Repeat what we learnt for a single raster, but with the temperature raster.
  2. (advanced) Get the range of monhly temperature within the 12 months before the observation.
Click to see the answer 1

Temperature are slow changing variables, so no real need for buffer, but this is a good exercice.

# create buffers
poly_otter <- buffer(pt_otter, 1000) # buffer of 1km
# no need to project because everything is in WGS84
# extract values for each layer
mean_buffer_temp <- extract(
  temp_C,
  poly_otter,
  fun = mean,
  exact = TRUE,
  ID = FALSE
)
boxplot(mean_buffer_temp, las = 2, ylab = "Temperature (*C)")

Click to see the answer 2 
# we start again from the matching time
t_obs <- match(pt_otter$time, names(pt_temp))
# then we want to get 11 month before the observation
t_past <- t_obs - 11 # ideally check if not negative

# calculate the difference of min-max value within he time period
temp_past12m <- sapply(1:nrow(pt_temp), function(i) {
  diff(range(pt_temp[i, t_past[i]:t_obs[i]]))
})

summary(temp_past12m)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  16.00   17.35   18.00   17.81   18.60   20.00 
# attach to the attribute table
pt_otter$temp_range_12m <- temp_past12m

# visualize the results
plot(pt_otter, "temp_range_12m")
plot(france_border, add = TRUE)