When faced with changing climate conditions, an organism is expected to track it’s climatic niche through space. Climate velocities are useful metrics that provide ecologists with a way to know how climatic conditions are shifting across a landscape. In this tutorial, we will learn how to calculate different climate velocity metrics using the VoCC package.

Types of climate velocity metrics

There are two types of climate velocities, each having slightly different meanings. Gradient-based velocities are calculated based on local climatic temporal and spatial gradients, while distance-based velocities are calculated based on distance to the nearest analogue (similar) climate.

Below we will discuss each in detail, calculating them as we go and discussing their ecological meanings and exploring their limitations.

1. Gradient-based climate velocities

Gradient-based climate velocities (gVoCC) are a way to describe the propensity of a species to move in response to climate change. They represent the ratio between the long-term temporal trend in climate conditions by the spatial gradient in climate conditions:

gVoCC = long-term trend / spatial gradient

A high gVoCC indicates that a species in that location would have a high propensity to move in response to climate change. When the gVoCC is high, it indicates that the climate is changing quickly over time relative to the amount of spatial variation in climate. High gVoCCs tend to occur in areas with heterogeneous climates that are experiencing high rates of temporal change, such as deserts. In areas with high gVoCCs, an organism would have to shift quickly and move far in order to reach an analogous climate.

On the other hand, a low gVoCC indicates that a species in that location would have a high propensity to move in response to climate change. In these areas, the climate is changing slowly relative to the amount of spatial variation in climate conditions. Areas like mountainous regions, where spatial variation in climates is high, tend to have low gVoCCs. In these areas, on organism would be less inclined to move quickly because a similar climate is closer by.

Let’s try calculating the gVoCC of temperature since 1880. We will use a raster stack of annual average air temperatures from Berkeley Earth:

## read in the raster stack:
temps = readRDS(paste(path, "/BerkEarth_rstack_mean.rds", sep = ""))

## plot the first layer (average air temperatures in 1880):
plot(temps[[1]], main = names(temps)[1])

## plot the last layer (average air temperatures in 2021):
plot(temps[[nlayers(temps)]], main = last(names(temps)))

First, we need to calculate the temporal trend in temperature across the period 1880-2021 for each raster cell. This is done by fitting a linear regression to the temperature time series from each grid cell. The function tempTrend allows us to do this easily:

## use function to calculate temporal trend in annual air temperature for each location:
## th = minimum number of observations in the series needed to calculate the trend at each cell
ttrend = tempTrend(r = temps,
                   th = 0.25*nlayers(temps) ## set minimum # obs. to 1/4 time series length
                   )

## plot it:
plot(ttrend)

You can see here that the temporal trends (slope of the regressions, slpTrends) show that the climate has warmed faster towards the poles.

Next, we need to calculate the spatial temperature gradient. This is done using the spatialGrad function. To use this function, we need to collapse our raster stack into a stack of annual temperature values for each year. I did this for you because it takes a while, so all we have to do is read in the raster:

## read in raster of mean annual temperatures:
mean_temps = readRDS(paste(path, "/BerkEarth_rstack_mean.rds", sep = ""))

## use function to calculate spatial gradient in mean daily air temperature for each location:
spgrad = spatGrad(r = mean_temps, 
                  projected = TRUE) ## our raster is projected to a coordinate system

## plot it: 
plot(spgrad)

Spatial resolution really matters for this step. Let’s see what happens when we resample our raster to a coarser resolution before calculating the spatial gradient:

## resample the raster stack to a coarser grain size:
rsamp_temps = aggregate(mean_temps, fact = 6, fun = mean)
plot(rsamp_temps)

## recalculate the spatial gradient
spgrad_coarse = spatGrad(r = rsamp_temps,
                   projected = TRUE) ## our raster is projected to a coordinate system


## plot it:
plot(spgrad_coarse)

max(spgrad$Grad)
## Warning in max(): Nothing to summarize if you provide a single RasterLayer; see
## cellStats
## class      : RasterLayer 
## dimensions : 180, 360, 64800  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : Grad 
## values     : 0, 0.1417279  (min, max)
max(spgrad_coarse$Grad)
## Warning in max(): Nothing to summarize if you provide a single RasterLayer; see
## cellStats
## class      : RasterLayer 
## dimensions : 30, 60, 1800  (nrow, ncol, ncell)
## resolution : 6, 6  (x, y)
## extent     : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : Grad 
## values     : 0, 0.04589104  (min, max)
## spatial averaging makes gradients far less step

Now that we have both the temporal trends and spatial gradients, we can compute their ratio to find the gradient-based climate velocity. The function gVoCC does this for us:

## calculate gradient based climate velocity:
gvocc = gVoCC(tempTrend = ttrend, spatGrad = spgrad)

## plot it:
plot(gvocc)

hist(gvocc$voccMag)

hist(log(gvocc$voccMag))

## plot on a log scale
plot(log(gvocc[[1]]))

Let’s make a prettier and more informative plot now.

## make a dataframe:
gvocc_df = data.frame(rasterToPoints(gvocc[[1]]))

gvocc_df %>%
  filter(!is.infinite(voccMag)) %>%
  filter(voccMag < 150) %>% ## remove outliers for better visualization
  ggplot(aes(x = x, y = y, fill = voccMag)) + 
  geom_raster() + 
  coord_fixed() +
  labs(fill = "Velocity of temperature",
       x = "Latitude", 
       y = "Longitude") +
  scale_fill_gradient(low = "yellow", high = "red", 
                      trans = "log", breaks = c(0, 1, 10, 100),
                      labels = c("0 km/year", "1 km/year", 
                                 "10 km/year", "100 km/year")) +
  theme_light()
## Warning: Raster pixels are placed at uneven horizontal intervals and will be shifted
## ℹ Consider using `geom_tile()` instead.


Gradient-based climate velocities can provide insights into climate connectivity, as well. gVoCC can be used to calculate climate trajectories, which show paths connecting present local climates with their future locations through the most direct route. Climate trajectories can be a useful tool to analyze how an organism would need to move across the landscape to track its climatic niche across space.

The function VoCCTraj allows us to easily calculate them:

## get velocity and angle separately
vel = gvocc[[1]]
ang = gvocc[[2]]

## calculate the mean temperatureover the period
mn <- mean(temps, na.rm = T)

## get the set of starting cells for the trajectories
lonlat <- na.omit(data.frame(xyFromCell(vel, 1:ncell(vel)), vel[], ang[], mn[]))[,1:2]

## calculate trajectory:
traj = voccTraj(lonlat, # starting lat lon coords
                vel, # velocity 
                ang, # angle
                mn, # mean temp 
                tyr = 71 # length of period of interest, years
                )
## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |===================================                                   |  50%
  |                                                                            
  |=================================================                     |  70%
  |                                                                            
  |========================================================              |  80%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |=============================================================         |  87%
  |                                                                            
  |==============================================================        |  89%
  |                                                                            
  |===============================================================       |  91%
  |                                                                            
  |================================================================      |  92%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |=================================================================     |  93%
  |                                                                            
  |==================================================================    |  94%
  |                                                                            
  |==================================================================    |  95%
  |                                                                            
  |===================================================================   |  95%
  |                                                                            
  |===================================================================   |  96%
  |                                                                            
  |====================================================================  |  97%
  |                                                                            
  |====================================================================  |  98%
  |                                                                            
  |===================================================================== |  98%
  |                                                                            
  |===================================================================== |  99%
  |                                                                            
  |======================================================================|  99%
  |                                                                            
  |======================================================================| 100%
## plot it:

## create a spatial line data frame from traj
lns <- trajLine(x = traj)
plot(lns)

## zoom in:
sp = st_as_sf(lns)
sp_cropped = st_crop(sp, xmin = -20, xmax = 45,
                                    ymin = 30, ymax = 73)
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
sp_cropped %>%
  ggplot(aes()) + geom_sf() + theme_light()

2. Distanced-based climate velocities

Distanced-based climate velocities (dVoCC) identify local environments that have future climatic conditions analogous to the baseline conditions at a location of interest, and then estimate the distance and direction between the location of interest and its future analogue. They are calculated as the geographical distance to the closest climate analogue for divided by the time elapsed between baseline and future periods:

dVoCC = distance to climate analogue / time between analogous climates

dVoCC measurements depend on both the algorithm used to identify climate analogues, and the threshold used to define an “analogous climate”. The algorithm decides a climate is analogous if deviations from the baseline climate remain within the threshold (defined by you). Defining these thresholds can be quite subjective, and can be especially difficult when multiple climate variables (e.g., temperature and rainfall) are used together. In these cases, algorithms might found that no climate analogues actually exist. A logical way to define the thresholds is based on natural variation in baseline conditions. For example, it would make sense to say a climate is not an analogue to another if its conditions have deviated more than one standard deviation from those in a baseline period of the climate of interest.

Let’s calculate dVoCC. To find analogous climates, we need data with high spatial resolution. For this, we will use a data set that is built in to the VoCC package.

## learn about data
?JapTC

Now, we need to make a very specifically structured data frame as input. The first column must be a climate variable at the present time, the second column must be the same variable at a future time. For now, let’s only use mean annual precipitation as our variable of interest (AnMn):

## extract values from cells and format
clim_df <- na.omit(data.frame(AnMn_hist = getValues(JapTC$AnMnPpr60_69), 
                           AnMn_fut = getValues(JapTC$AnMnPpr08_17), 
                           cid = 1:ncell(JapTC)))
clim_df[,c("x","y")] <- xyFromCell(JapTC, clim_df$cid)

Now we can run the dVoCC function:

## now we can run the dVoCC function 
dvocc  = dVoCC(clim = clim_df,
               n = 1, # number of climate variables (for now we are only using precipitation)
               tdiff = 40, # number of years between time periods 
               method = "Single", # specificy that all cells will have the same threshold to be analgous
               distfun = "GreatCircle", # specify function used to calculate distances to climate analogues
               lonlat = TRUE, # not projected coordinates
               geoTol = 160, # set max distance to analogue
               climTol = 10 # set threshold to 10 mm deviation in precipitation  
)

## plot 
r <- raster(JapTC)
r[dvocc$focal] <- dvocc$vel
plot(r)

Let’s try making the analogue threshold location-specific. To do this, we need to add a third column describing the threshold for each cell and change the method argument to "Variable". In this case, we can use the standard deviation of precipitation:

## extract values from cells and format
clim_df <- na.omit(data.frame(AnMn_hist = getValues(JapTC$AnMnPpr60_69), 
                           AnMn_fut = getValues(JapTC$AnMnPpr08_17), 
                           AnMn_sd = getValues(JapTC$AnMnSDPpr60_69),
                           cid = 1:ncell(JapTC)))
clim_df[,c("x","y")] <- xyFromCell(JapTC, clim_df$cid)

## now we can run the dVoCC function 
dvocc_thresh  = dVoCC(clim = clim_df,
               n = 1, # number of climate variables (for now we are only using precipitation)
               tdiff = 40, # number of years between time periods 
               method = "Variable", # specificy that all cells will have different thresholds 
               distfun = "GreatCircle", # specify function used to calculate distances to climate analogues
               lonlat = TRUE, # not projected coordinates
               geoTol = 160, # set max distance to analogue
               climTol = NA 
)

## plot 
r_thresh <- raster(JapTC)
r_thresh[dvocc_thresh$focal] <- dvocc_thresh$vel
plot(r_thresh)

Now, let’s add more variables. We also have maximum and minimum temperature (Tmax and Tmin).

## extract values from cells and format
clim_df <- na.omit(data.frame(AnMn_hist = getValues(JapTC$AnMnPpr60_69), 
                           AnMn_fut = getValues(JapTC$AnMnPpr08_17), 
                           Tmax_hist = getValues(JapTC$AnMnTmax60_69),
                           Tmax_fut = getValues(JapTC$AnMnTmax08_17),
                           Tmin_hist = getValues(JapTC$AnMnTmin60_69),
                           Tmin_fut = getValues(JapTC$AnMnTmin08_17),
                           cid = 1:ncell(JapTC)))
clim_df[, c("x","y")] <- xyFromCell(JapTC, clim_df$cid)

## now we can run the dVoCC function 
dvocc_multi  = dVoCC(clim = clim_df,
               n = 3, # number of climate variables: now 3
               tdiff = 40, # number of years between time periods 
               method = "Single", # specificy that all cells will have same threshold 
               distfun = "GreatCircle", # specify function used to calculate distances to climate analogues
               lonlat = TRUE, # not projected coordinates
               geoTol = 160, # set max distance to analogue
               climTol = c(10, 0.1, 0.1) # set 10mm threshold for precip, 0.1C threshold for temp
)

## plot 
r_multi = raster(JapTC)
r_multi[dvocc_multi$focal] <- dvocc_multi$vel
plot(r_multi)

The function’s argument geoTol can be used to specific a certain distance that a climate analogue cell must be within from the focal cell. This argument can be useful in understanding how organisms with different dispersal distances might be able to find. An organism with a small dispersal distance (low geoTol) might not be able to find a climate analog, whereas an organism with a large dispersal distance (high geoTol) in the same location might.

Let’s compare those situations:

## low dispersal distance:
## use multi variable climate data
dvocc_lowdisp  = dVoCC(clim = clim_df,
               n = 3, # number of climate variables: now 3
               tdiff = 40, # number of years between time periods 
               method = "Single", # specificy that all cells will have same threshold 
               distfun = "GreatCircle", # specify function used to calculate distances to climate analogues
               lonlat = TRUE, # not projected coordinates
               geoTol = 30, # set max distance to analogue
               climTol = c(10, 0.1, 0.1) # set 10mm threshold for precip, 0.1C threshold for temp
)
## high dispersal distance:
dvocc_highdisp  = dVoCC(clim = clim_df,
               n = 3, # number of climate variables: now 3
               tdiff = 40, # number of years between time periods 
               method = "Single", # specificy that all cells will have same threshold 
               distfun = "GreatCircle", # specify function used to calculate distances to climate analogues
               lonlat = TRUE, # not projected coordinates
               geoTol = 200, # set max distance to analogue
               climTol = c(10, 0.1, 0.1) # set 10mm threshold for precip, 0.1C threshold for temp
)

## plot 
r_lowdisp = raster(JapTC)
r_lowdisp[dvocc_lowdisp$focal] <- dvocc_lowdisp$vel

r_highdisp = raster(JapTC)
r_highdisp[dvocc_highdisp$focal] <- dvocc_highdisp$vel

disp = stack(r_lowdisp, r_highdisp)
names(disp) = c("low_dispersal", "high_dispersal")

plot(disp)

Now let’s plot all the places where a species with poor dispersal ability would not be able to reach a climate analogue, but a species with good dispersal ability would:

## plot places where r_lowdisp is NA but r_highdisp is not 
plot(is.na(r_lowdisp) & !is.na(r_highdisp))

This tutorial was developed by Nikki Moore for the BIOS2 working group “Assessing the potential for climate-driven range shifts through multiple landscapes features across the Canada-US border” based on information found in Molinos et al. 2019 (DOI: 10.1111/2041-210X.13295).