bioticVelocity: Velocity of shifts in densities across a series of rasters

Description

Calculates metrics of "movement" of cell densities across a time series of rasters. Rasters could represent, for example, the probability of presence of a species through time. In this case, velocities would indicate rates and directions of range shift. The simplest metric is movement of the density-weighted centroid (i.e., range "center"), but many more are available to provide a nuanced indicator of velocity. See Details for the types of metrics that can be calculated.

Usage

bioticVelocity(
  x,
  times = NULL,
  atTimes = NULL,
  elevation = NULL,
  metrics = c("centroid", "nsCentroid", "ewCentroid", "nCentroid", "sCentroid",
    "eCentroid", "wCentroid", "nsQuants", "ewQuants", "similarity", "summary"),
  quants = c(0.05, 0.1, 0.5, 0.9, 0.95),
  onlyInSharedCells = FALSE,
  cores = 1,
  warn = TRUE,
  longitude = NULL,
  latitude = NULL,
  paths = NULL,
  ...
)

Value

A data frame with biotic velocities and related values. Fields are as follows:

timeFrom: Start time of interval
timeTo: End time of interval
timeMid: Time point between timeFrom and timeTo
timeSpan: Duration of interval

Depending on metrics that are specified, additional fields are as follows. All measurements of velocity are in distance units (typically meters) per time unit (which is the same as the units used for times and atTimes). For example, if the rasters are in an Albers equal-area projection and times are in years, then the output will be meters per year.

If metrics has 'centroid': Columns named centroidVelocity, centroidLong, centroidLat -- Speed of weighted centroid, plus its longitude and latitude (in the timeTo period of each time step). Values are always >= 0.
If metrics has 'nsCentroid': Columns named nsCentroid and nsCentroidLat -- Velocity of weighted centroid in north-south direction, plus its latitude (in the timeTo period of each time step). Positive values connote movement north, and negative values south.
If metrics has 'ewCentroid': ewCentroid and ewCentroidLong -- Velocity of weighted centroid in east-west direction, plus its longitude (in the timeTo period of each time step). Positive values connote movement east, and negative values west.
If metrics has 'nCentroid', 'sCentroid', 'eCentroid', and/or 'wCentroid': Columns named nCentroidVelocity and nCentroidAbund, sCentroid and sCentroidAbund, eCentroid and eCentroidAbund, and/or wCentroid and wCentroidAbund -- Speed of weighted centroid of all cells that fall north, south, east, or west of the landscape-wide centroid, plus a column indicating the total weight (abundance) of all such populations. Values are always >= 0.
If metrics contains any of nsQuants or ewQuants: Columns named nsQuantVelocity_quantQ and nsQuantLat_quantQ, or ewQuantVelocity_quantQ and ewQuantLat_quantQ: Velocity of the Qth quantile weight in the north-south or east-west directions, plus the latitude or longitude thereof (in the timeTo period of each time step). Quantiles are cumulated starting from the south or the west, so the 0.05th quantile, for example, is in the far south or west of the range and the 0.95th in the far north or east. Positive values connote movement north or east, and negative values movement south or west.
If metrics contains similarity, metrics of similarity are calculated for each pair of successive landscapes, defined below as x1 (raster in timeFrom) and x2 (raster in timeTo), with the number of shared non-NA cells between them being n:
- A column named simpleMeanDiff: sum(x2 - x1, na.rm = TRUE) / n
- A column named meanAbsDiff: sum(abs(x2 - x1), na.rm = TRUE) / n
- A column named rmsd (root-mean square difference): sqrt(sum((x2 - x1)^2, na.rm = TRUE)) / n
- A column named godsoeEsp: 1 - sum(2 * (x1 * x2), na.rm=TRUE) / sum(x1 + x2, na.rm = TRUE), values of 1 ==> maximally similar, 0 ==> maximally dissimilar.
- A column named schoenersD: 1 - (sum(abs(x1 - x2), na.rm = TRUE) / n), values of 1 ==> maximally similar, 0 ==> maximally dissimilar.
- A column named warrensI: 1 - sqrt(sum((sqrt(x1) - sqrt(x2))^2, na.rm = TRUE) / n), values of 1 ==> maximally similar, 0 ==> maximally dissimilar.
- A column named cor: Pearson correlation between x1 and x2.
- A column named rankCor: Spearman rank correlation between x1 and x2.
If metrics contains elevCentroid: Columns named elevCentroidVelocity and elevCentroidElev -- Velocity of the centroid in elevation (up or down) and the elevation in the "to" timestep. Positive values of velocity connote movement upward, and negative values downward.
If metrics contains elevQuants: Columns named elevQuantVelocity_quantQ and elevQuantVelocityElev_quantQ -- Velocity of the Nth quantile of mass in elevation (up or down) and the elevation of this quantile in the "to" timestep. Positive values of velocity connote movement upward, and negative values downward.
If metrics contains summary:
- A column named propSharedCellsNotNA: Proportion of cells that are not NA in both the "from" and "to" time steps. The proportion is calculated using the total number of cells in a raster as the denominator (i.e., not total number of cells across two rasters).
- Columns named timeFromPropNotNA and timeToPropNotNA: Proportion of cells in the "from" time and "to" steps that are not NA.
- A column named mean: Mean weight in timeTo time step. In the same units as the values of the cells.
- Columns named quantile_quantQ: The Qth quantile(s) of weight in the timeTo time step. In the same units as the values of the cells.
- A column named prevalence: Proportion of non-NA cells with weight >0 in the timeTo time step relative to all non-NA cells. Unitless.

Arguments

x

Either a SpatRaster or a 3-dimensional array. Values should really be be either NA or >= 0.

If x is a SpatRaster, then each layer is assumed to represent a time slice. Rasters must be in an equal-area projection. They must also be ordered temporally, with the raster "on top" assumed to represent the starting time.
If x is an array then each "layer" in the third dimension is assumed to represent a map at a particular time slice in an equal-area projection. Note that if this is an array you should probably also specify the arguments longitude and latitude.

times

Numeric vector with the same number of layers in x or NULL (default). This specifies the time represented by each layer in x from beginning of the time series (top layer) to the end (bottom layer). Times must appear in sequential order. For example, if time periods are 24 kybp, 23 kybp, 22 kybp, use c(-24, -23, -22), not c(24, 23, 22). If NULL (default), values are assigned starting at 1 and ending at the total number of layers in x.

atTimes

Numeric, values of times across which to calculate biotic velocity. You can use this to calculate biotic velocities across selected time periods (e.g., just the first and last time periods). Note that atTimes must be the same as or a subset of times. The default is NULL, in which case velocity is calculated across all time slices (i.e., between times 1 and 2, 2 and 3, 3 and 4, etc.).

elevation

Either NULL (default) or a raster or matrix representing elevation. If this is supplied, changes in elevation are incorporated into all velocity and speed metrics. Additionally, you can also calculate the metrics elevCentrioid and elevQuants.

metrics

Biotic velocity metrics to calculate (default is to calculate them all). All metrics ignore NA cells in x. Here, "starting time period" represents one layer in x and "end time period" the next layer.

centroid: Speed of mass-weighted centroid (directionless).
nsCentroid or ewCentroid: Velocity in the north-south or east-west directions of the mass-weighted centroid.
nCentroid, sCentroid, eCentroid, and wCentroid: Speed of mass-weighted centroid of the portion of the raster north/south/east/west of the landscape-wide weighted centroid of the starting time period.
nsQuants or ewQuants: Velocity of the location of the Qth quantile of mass in the north-south or east-west directions. The quantiles can be specified in quants. For example, this could be the movement of the 5th, 50th, and 95th quantiles of population size going from south to north. The 0th quantile would measure the velocity of the southernmost or easternmost cell(s) with values >0, and the 100th quantile the northernmost or westernmost cell(s) with non-zero values.
similarity: Metrics of similarity between each time period. Some of these make sense only for cases where values in x are in the range [0, 1], but not if some values are outside this range. See nicheOverlapMetrics for more details. The metrics are:
- Simple mean difference
- Mean absolute difference
- Root-mean squared difference
- Expected Fraction of Shared Presences or ESP (Godsoe, W. 2014. Ecography 37:130-136 tools:::Rd_expr_doi("10.1111/j.1600-0587.2013.00403.x"))
- D statistic (Schoener, T.W. 1968. Ecology 49:704-726. tools:::Rd_expr_doi("10.2307/1935534"))
- I statistic (Warren, D.L., et al. 2008. Evolution 62:2868-2883 tools:::Rd_expr_doi("10.1111/j.1558-5646.2008.00482.x"))
- Pearson correlation
- Spearman rank correlation
summary: This calculates a series of measures for each "starting time period" raster. None of these are measures of velocity:
- Mean: Mean value across all cells.
- Sum: Total across all cells.
- Quantiles: Qth quantile values across all cells. Quantiles are provided through argument quants.
- Prevalence: Number of cells with values > 0.
elevCentroid: Velocity of the centroid of mass in elevation (up or down). A raster or matrix must be supplied to argument elevation.
elevQuants: Velocity of the Qth quantile of mass in elevation (up or down). The quantiles to be evaluated are given by quants. The lowest elevation with mass >0 is the 0th quantile, and the highest elevation with mass >0 is the 100th. Argument elevation must be supplied.

quants

Numeric vector indicating the quantiles at which biotic velocity is calculated for the "quant" and "Quants" metrics. Default quantiles to calculate are c(0.1, 0.9).

onlyInSharedCells

If TRUE, calculate biotic velocity using only those cells that are not NA in the start and end of each time period. This is useful for controlling for shifting land mass due to sea level rise, for example, when calculating biotic velocity for an ecosystem or a species. The default is FALSE, in which case velocity is calculated using all cells in each time period, regardless of whether some become NA or change from NA to not NA.

cores

Positive integer. Number of processor cores to use. Note that if the number of time steps at which velocity is calculated is small, using more cores may not always be faster. If you have issues when cores > 1, please see the troubleshooting_parallel_operations guide.

warn

Logical, if TRUE (default) then display function-specific warnings.

longitude

Numeric matrix or NULL (default):

If x is a SpatRaster, then this is ignored (longitude is ascertained directly from the rasters, which must be in equal-area projection for velocities to be valid).
If x is an array and longitude is NULL (default), then longitude will be ascertained from column numbers in x and velocities will be in arbitrary spatial units (versus, for example, meters). Alternatively, this can be a two-dimensional matrix whose elements represent the longitude coordinates of the centers of cells of x. The matrix must have the same number of rows and columns as x. Coordinates must be from an equal-area projection for results to be valid.

latitude

Numeric matrix or NULL (default):

If x is a SpatRaster, then this is ignored (latitude is obtained directly from the rasters, which must be in equal-area projection for velocities to be valid).
If x is an array and latitude is NULL (default), then latitude will be obtained from row numbers in x and velocities will be in arbitrary spatial units (versus, for example, meters). Alternatively, this can be a two-dimensional matrix whose elements represent the latitude coordinates of the centers of cells of x. The matrix must have the same number of rows and columns as x. Coordinates must be from an equal-area projection for results to be valid.

paths

This is used internally and rarely (never?) needs to be defined by a user (i.e., leave it as NULL). Valid values are a character vector or NULL (default). If a character vector, it should give the values used by .libPaths.

...

Other arguments (not used).

Details

Attention:

This function may yield erroneous velocities if the region of interest is near or spans a pole or the international date line. Results using the "Quant" and "quant" metrics may be somewhat counterintuitive if just one cell is >0, or one row or column has the same values with all other values equal to 0 or NA because defining quantiles in these situations is not intuitive. Results may also be counterintuitive if some cells have negative values because they can "push" a centroid away from what would seem to be the center of mass as assessed by visual examination of a map.

Note:

For the nsQuants and ewQuants metrics it is assumed that the latitude/longitude assigned to a cell is at its exact center. For calculating the position of a quantile, density is interpolated linearly from one cell center to the center of the adjacent cell. If a desired quantile does not fall exactly on the cell center, it is calculated from the interpolated values. For quantiles that fall south/westward of the first row/column of cells, the cell border is assumed to be at 0.5 * cell length south/west of the cell center.

Examples

Run this code

# \donttest{
# NB These examples can take a few minutes to run.
# To illustrate calculation and interpretation of biotic velocity,
# we will calibrate a SDM for the Red-Bellied Lemur and project
# the model to the present and successive future climates. The time series
# of rasters is then used to calculate biotic velocity.

library(sf)
library(terra)

### process environmental rasters
#################################

# get rasters
rastFile <- system.file('extdata/madClim.tif', package='enmSdmX')
madClim <- rast(rastFile)

rastFile <- system.file('extdata/madClim2030.tif', package='enmSdmX')
madClim2030 <- rast(rastFile)

rastFile <- system.file('extdata/madClim2050.tif', package='enmSdmX')
madClim2050 <- rast(rastFile)

rastFile <- system.file('extdata/madClim2070.tif', package='enmSdmX')
madClim2070 <- rast(rastFile)

rastFile <- system.file('extdata/madClim2090.tif', package='enmSdmX')
madClim2090 <- rast(rastFile)

# The bioticVelocity() function needs rasters to be in equal-area
# projection, so we will project them here.
madAlbers <- getCRS('madAlbers') # Albers projection for Madagascar
madClim <- project(madClim, madAlbers)
madClim2030 <- project(madClim2030, madAlbers)
madClim2050 <- project(madClim2050, madAlbers)
madClim2070 <- project(madClim2070, madAlbers)
madClim2090 <- project(madClim2090, madAlbers)

# Coordinate reference systems:
wgs84 <- getCRS('WGS84') # WGS84
madAlbers <- getCRS('madAlbers') # Madagascar Albers

# lemur occurrence data
data(lemurs)
occs <- lemurs[lemurs$species == 'Eulemur fulvus', ]
occs <- vect(occs, geom=c('longitude', 'latitude'), crs=wgs84)
occs <- project(occs, madAlbers)

# eliminate cell duplicates
occs <- elimCellDuplicates(occs, madClim)

# extract environment at occurrences
occEnv <- extract(madClim, occs, ID = FALSE)
occEnv <- occEnv[complete.cases(occEnv), ]
	
# create background sites (using just 1000 to speed things up!)
bgEnv <- terra::spatSample(madClim, 3000)
bgEnv <- bgEnv[complete.cases(bgEnv), ]
bgEnv <- bgEnv[sample(nrow(bgEnv), 1000), ]

# collate occurrences and background sites
presBg <- data.frame(
   presBg = c(
      rep(1, nrow(occEnv)),
      rep(0, nrow(bgEnv))
   )
)

env <- rbind(occEnv, bgEnv)
env <- cbind(presBg, env)

### calibrate model
###################

predictors <- c('bio1', 'bio12')

# MaxEnt
mx <- trainMaxEnt(
	data = env,
	resp = 'presBg',
	preds = predictors,
	regMult = 1, # too few values for reliable model, but fast
	cores = 2
)

### project to present and future climate
#########################################

predPresent <- predictEnmSdm(mx, madClim)
pred2030 <- predictEnmSdm(mx, madClim2030)
pred2050 <- predictEnmSdm(mx, madClim2050)
pred2070 <- predictEnmSdm(mx, madClim2070)
pred2090 <- predictEnmSdm(mx, madClim2090)

plot(predPresent, main = 'Present Suitability')

# plot change in suitability between present and 2090s
delta <- pred2090 - predPresent
plot(delta, main = 'Change in Suitability')

### calculate biotic velocity
#############################

series <- c(
	predPresent,
	pred2030,
	pred2050,
	pred2070,
	pred2090
)

names(series) <- c('present', 't2030', 't2050', 't2070', 't2090')
plot(series)

times <- c(1985, 2030, 2050, 2070, 2090)
quants <- c(0.10, 0.90)

bv <- bioticVelocity(
	x = series,
	times = times,
	quants = quants,
	cores = 2
)
 
bv

### centroid velocities

# centroid (will always be >= 0)
# fastest centroid movement around 2060
plot(bv$timeMid, bv$centroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Speed (m / y)', main = 'Centroid Speed')
  
# velocity northward/southward through time
# shows northward shift because always positive, fastest around 2060
plot(bv$timeMid, bv$nsCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)', main = 'Centroid N/S Velocity')
  
# velocity eastward (positive)/westward (negative) through time
# movement eastward (positive) first, then westward (negative)
plot(bv$timeMid, bv$ewCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)', main = 'Centroid E/W Velocity')

### map of centroid location through time
# shows centroid moves slightly northward through time
plot(delta, main = 'Centroid Location &\nChange in Suitability')
points(bv$centroidLong[1], bv$centroidLat[1], pch = 1)
points(bv$centroidLong[4], bv$centroidLat[4], pch = 16)
lines(bv$centroidLong, bv$centroidLat)
legend(
  'bottomright',
  legend = c(
    'start (~1985)',
	'stop (~2090)',
	'trajectory'
  ),
  pch = c(1, 16, NA),
  lwd = c(NA, NA, 1)
)

### velocities of portions of range north/south/east/west of centroid
# positive ==> northward shift
# negative ==> southward shift
  
# portion of range north of centroid
# shows northward expansion because always positive
plot(bv$timeMid, bv$nCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Northern Part of Range')
  
# portion of range south of centroid
# shows northward contraction because always positive
plot(bv$timeMid, bv$sCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Southern Part of Range')
  
# portion of range east of centroid
# shows eastern portion moves farther east
plot(bv$timeMid, bv$eCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Eastern Part of Range')
  
# portion of range west of centroid
# shows western portion moves east
plot(bv$timeMid, bv$wCentroidVelocity, type = 'l',
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Western Part of Range')

### velocities of range margins

# from south to north, 10th and 90th quantiles of density
# positive ==> northward shift
# negative ==> southward shift
# shows both northern and southern range margins shift northward
# because always positive... northern margin shift usually slower
ylim <- range(bv$nsQuantVelocity_quant0p1, bv$nsQuantVelocity_quant0p9)

plot(bv$timeMid, bv$nsQuantVelocity_quant0p1, type = 'l', ylim = ylim,
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Northern/Southern Range Margins')
lines(bv$timeMid, bv$nsQuantVelocity_quant0p9, lty = 'dashed')
legend(
  'bottomright',
  legend = c('Southern Margin', 'Northern Margin'),
  lty = c('solid', 'dashed')
)
  
# from east to west, 10th and 90th quantiles of density
# positive ==> eastward shift
# negative ==> westward shift
ylim <- range(bv$ewQuantVelocity_quant0p1, bv$ewQuantVelocity_quant0p9)

plot(bv$timeMid, bv$ewQuantVelocity_quant0p1, type = 'l', ylim = ylim,
  xlab = 'Year', ylab = 'Velocity (m / y)',
  main = 'Eastern/Western Range Margins')
lines(bv$timeMid, bv$ewQuantVelocity_quant0p9, lty = 'dashed')
legend(
  'bottomright',
  legend = c('Eastern Margin', 'Western Margin'),
  lty = c('solid', 'dashed')
)
  
  
### summary statistics

# mean density across cells through time
plot(bv$timeMid, bv$mean, type = 'l',
  xlab = 'Year', ylab = 'Mean Density',
  main = 'Mean Density')

# sum of density across cells through time
plot(bv$timeMid, bv$sum, type = 'l',
  xlab = 'Year', ylab = 'Sum of Density',
  main = 'Sum of Density')

### change metrics

# average change in suitability from one time period to next
# shows average conditions getting worse
plot(bv$timeMid, bv$simpleMeanDiff, type = 'l',
  xlab = 'Year', ylab = 'Mean Change in Suitability')
  
# average absolute change in suitability from one time period to next
# shows average absolute change declining
plot(bv$timeMid, bv$meanAbsDiff, type = 'l',
  xlab = 'Year', ylab = 'Mean Absolute Change in Suitability')
  
# root-mean square difference from one time period to the next
# shows difference between successive rasters declines through time
plot(bv$timeMid, bv$rmsd, type = 'l',
  xlab = 'Year', ylab = 'RMSD')
  
### raster similarity
# most indicate that successive rasters are similar through time
ylim <- range(bv$godsoeEsp, bv$schoenerD, bv$warrenI, bv$cor, bv$warrenI)
plot(bv$timeMid, bv$godsoeEsp, type = 'l', lty = 1, col = 1,
  xlab = 'Year', ylab = 'Raster similarity', ylim = ylim)
lines(bv$timeMid, bv$schoenerD, lty = 2, col = 2)
lines(bv$timeMid, bv$warrenI, lty = 3, col = 3)
lines(bv$timeMid, bv$cor, lty = 4, col = 4)
lines(bv$timeMid, bv$rankCor, lty = 5, col = 5)

legend(
  'right',
  legend = c(
    'Godsoe\'s ESP',
    'Schoener\'s D',
    'Warren\'s I',
    'Correlation',
    'Rank Correlation'
  ),
  col = 1:5,
  lty = 1:5
)
  

# values of 10th and 90th quantiles across cells through time
# shows most favorable cells becoming less favorable
# least favorable cells remain mainly unchanged
ylim <- range(bv$quantile_quant0p1, bv$quantile_quant0p9)

plot(bv$timeMid, bv$quantile_quant0p1, type = 'l', ylim = ylim,
  xlab = 'Year', ylab = 'Quantile Value',
  main = 'Quantiles across Cells')
lines(bv$timeMid, bv$quantile_quant0p9, lty = 'dashed')

legend(
  'topright',
  legend = c('10th quantile', '90th quantile'),
  lty = c('solid', 'dashed')
)

### map of northern/southern range margins through time

# range of longitude shown in plot
madExtent <- ext(madClim)
xExtent <- as.vector(madExtent)[1:2]

plot(predPresent, main = 'North/South Range Margin Location')
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p9[1], bv$nsQuantLat_quant0p9[1]))
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p9[2], bv$nsQuantLat_quant0p9[2]), lty = 'dashed')
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p9[3], bv$nsQuantLat_quant0p9[3]), lty = 'dotdash')
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p9[4], bv$nsQuantLat_quant0p9[4]), lty = 'dotted')

lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p1[1], bv$nsQuantLat_quant0p1[1]))
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p1[2], bv$nsQuantLat_quant0p1[2]), lty = 'dashed')
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p1[3], bv$nsQuantLat_quant0p1[3]), lty = 'dotdash')
lines(c(xExtent[1], xExtent[2]),
  c(bv$nsQuantLat_quant0p1[4], bv$nsQuantLat_quant0p1[4]), lty = 'dotted')

legend(
  'bottomright',
  legend = c(
    '1980s',
	'2030s',
	'2050s',
	'2070s',
	'2090s'
  ),
  lty = c('solid', 'dashed', 'dotdash', 'dotted')
)

### map of eastern/western range margins through time

# range of longitude shown in plot
madExtent <- ext(madClim)
yExtent <- as.vector(madExtent)[3:4]

plot(predPresent, main = 'North/South Range Margin Location')
lines(c(bv$ewQuantLong_quant0p9[1], bv$ewQuantLong_quant0p9[1]),
  c(yExtent[1], yExtent[2]))
lines(c(bv$ewQuantLong_quant0p9[2], bv$ewQuantLong_quant0p9[2]),
  c(yExtent[1], yExtent[2]), lty = 'dashed')
lines(c(bv$ewQuantLong_quant0p9[3], bv$ewQuantLong_quant0p9[3]),
  c(yExtent[1], yExtent[2]), lty = 'dotdash')
lines(c(bv$ewQuantLong_quant0p9[4], bv$ewQuantLong_quant0p9[4]),
  c(yExtent[1], yExtent[2]), lty = 'dotted')

lines(c(bv$ewQuantLong_quant0p1[1], bv$ewQuantLong_quant0p1[1]),
  c(yExtent[1], yExtent[2]))
lines(c(bv$ewQuantLong_quant0p1[2], bv$ewQuantLong_quant0p1[2]),
  c(yExtent[1], yExtent[2]), lty = 'dashed')
lines(c(bv$ewQuantLong_quant0p1[3], bv$ewQuantLong_quant0p1[3]),
  c(yExtent[1], yExtent[2]), lty = 'dotdash')
lines(c(bv$ewQuantLong_quant0p1[4], bv$ewQuantLong_quant0p1[4]),
  c(yExtent[1], yExtent[2]), lty = 'dotted')

legend(
  'bottomright',
  legend = c(
    '1980s',
	'2030s',
	'2050s',
	'2070s',
	'2090s'
  ),
  lty = c('solid', 'dashed', 'dotdash', 'dotted')
)

# }

Run the code above in your browser using DataLab