How do you extract elevation values from WorldClim 2.1?

52 views Asked by At

I have been using this helpful link to learn to extract bioclimatic variables from WorldClim 2.1 in R: https://www.benjaminbell.co.uk/2018/01/extracting-data-and-making-climate-maps.html?lr=1

When extracting values for elevation from 14 GPS points, I have found that the elevations WorldClim 2.1 finds do not match what I have obtained on GPSs or on Google Earth. They are sometimes off by several hundred meters. Below is my code.

# Set working directory to look at elevation (this contains "wc2.1_30s_elev.tif" from WorldClim 2.1):
setwd("/Users/XXXX/Desktop/WorldClim 2.1")

# Create raster layers for elevation
elev <- raster("wc2.1_30s_elev.tif")

# Create a data frame with sample site coordinates for Costa Rica and Panama
site <- c("PN Volcan Tenorio", "Rara Avis", "RB El Copal", "PN La Cangreja", "La Chaqueta", "San Isidro", "Las Cruces", "Risco Abajo", "La Fortuna", "Santa Fe", "Altos de Maria", "Cerro Jefe", "Chucantí", "Altos de Campana")
lon <- c(-84.99, -84.04, -83.75, -84.36, -83.91, -83.62, -82.97, -82.50, -82.20, -81.14, -80.07, -79.40, -78.45, -79.92)
lat <- c(10.70, 10.28, 9.78, 9.71, 9.50, 9.43, 8.78, 9.19, 8.72, 8.52, 8.64, 9.23, 8.79, 8.69)

elevationsamples <- data.frame(site, lon, lat, row.names="site")

# Here, we created three vector objects containing our data, then we created the data frame object by telling R which vectors to use to construct the data frame. 
# The vectors will become columns in the data frame. 
# Specifying row.names="site" tells R that the "site" vector should be used for the row names in the data frame. 
elevationsamples

# Extract data from WorldClim for your sites
elev.data <- elevationsamples 
elev.data$elevation <- extract(elev, elevationsamples)

# Check to see what it looks like
elev.data

I am wondering if I am doing something wrong while extracting, or if anyone knows of issues with the elevation values obtained by WorldClim 2.1.

1

There are 1 answers

0
Robert Hijmans On

This is not a coding question (and therefore inappropriate for this site) as the below (your code with some example data) works as expected.

library(terra)
x <- rast(names="elevation")
values(x) <- 1:ncell(x)

# Create a data frame with sample site coordinates for Costa Rica and Panama
site <- c("PN Volcan Tenorio", "Rara Avis", "RB El Copal", "PN La Cangreja", "La Chaqueta", "San Isidro", "Las Cruces", "Risco Abajo", "La Fortuna", "Santa Fe", "Altos de Maria", "Cerro Jefe", "Chucantí", "Altos de Campana")
lon <- c(-84.99, -84.04, -83.75, -84.36, -83.91, -83.62, -82.97, -82.50, -82.20, -81.14, -80.07, -79.40, -78.45, -79.92)
lat <- c(10.70, 10.28, 9.78, 9.71, 9.50, 9.43, 8.78, 9.19, 8.72, 8.52, 8.64, 9.23, 8.79, 8.69)

elev <- data.frame(site, lon, lat, row.names="site")
elev$elevation <- extract(x, elev, ID=FALSE)

head(elev)
#                     lon   lat elevation
#PN Volcan Tenorio -84.99 10.70     28536
#Rara Avis         -84.04 10.28     28536
#RB El Copal       -83.75  9.78     28897
#PN La Cangreja    -84.36  9.71     28896
#La Chaqueta       -83.91  9.50     28897
#San Isidro        -83.62  9.43     28897

Elevation data not matching is to be expected as there can be a lot of variation in elevation in grid cells of almost 1 by 1 km. In steep mountain areas a couple of 100m can certainly occur.