Introduction to elevatr

Elevation data is used for a wide array of applications, including, for example, visualization, hydrology, and ecological modelling. Gaining access to these data in R has not had a single interface, is made available through functions across many packages, or requires local access to the data. This is no longer required as a variety of APIs now exist that provide programmatic access to elevation data. The elevatr package was written to standarize access to elevation data from web APIs. This introductory vignette provides details on how to use elevatr to access elevation data and provides a bit of detail on the source data it accesses.

As of version 0.4.2, there are several endpoints that elevatr accesses. For point elevation data it uses USGS Elevation Point Query Service (United States only) as well as the Amazon Web Services (AWS) Terrain Tiles from which point elevations are extracted. Raster elevation data (i.e., Digital Elevation Models or DEMs) are available from the AWS Terrain Tiles or from the OpenTopography Global DEM API (https://portal.opentopography.org/apidocs/#/Public/getGlobalDem) . Currently, elevatr supports the SRTMGL3, SRTMGL1, AW3D30, and SRTM15Plus datasets.

Get Point Elevation Data

Point elevation is accessed from get_elev_point(). This function takes either a data.frame with x (longitude) and y (latitude) locations as the first two columns or a SpatialPoints/SpatialPointsDataFrame as input and then fetches the reported elevation for that location. As mentioned there is one service that provides this information. Details are provided below.

USGS Elevation Point Query Service

The USGS Elevation Point Query Service is accessible from elevatr. It is only available for the United States (including Alaska and Hawaii). Points that fall within the United States but are not on land return a value of zero. Points outside the United States boundaries return a value of -1000000.

Using get_elev_point() to Access The USGS Elevation Point Query Service

Usage of get_elev_point() requires an input SpatialPoints, SpatialPointsDataFrame, or a two-column data frame with column one containing the x (e.g. longitude) coordinates and the second column containing the y coordinates (e.g. latitude). The source data are global and also include estimates of depth for oceans.

Example usage of each is included below. For these examples, we can create a dataset to use.

# Create an example data.frame
set.seed(65.7)
examp_df <- data.frame(x = runif(3, min = -73, max = -72.5), y = runif(3, min = 42,
    max = 43))
prj_dd <- "EPSG:4326"

# Create and example data.frame with additional columns
cats <- data.frame(category = c("H", "M", "L"))

examp_df2 <- data.frame(examp_df, cats)

# Create an example SpatialPoints
examp_sp <- SpatialPoints(examp_df, proj4string = CRS(prj_dd))

# Create an example SpatialPointsDataFrame
examp_spdf <- SpatialPointsDataFrame(examp_sp, data = cats)

If a data frame is used it may have additional columns beyond the first two, which must contain the coordinates. The additional columns, along with the returned elevation, will be part of the output SpatialPointsDataFrame. Similarly, an elevation column is added to the data slot of a SpatialPointsDataFrame.

The USGS Elevation Point Query Service returns a single point at a time. The implemntation in get_elev_point() will loop through each point, thus can be slow for large number of requests.

Accessing data from this service is done by setting the src to "epqs". No API key is required and there are no rate limits.

df_elev_epqs <- get_elev_point(examp_df, prj = prj_dd, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
data.frame(df_elev_epqs)
##   elevation elev_units coords.x1 coords.x2 optional
## 1    207.81     meters -72.94216  42.04268     TRUE
## 2    198.95     meters -72.65599  42.60974     TRUE
## 3    306.15     meters -72.91468  42.26707     TRUE
df2_elev_epqs <- get_elev_point(examp_df2, prj = prj_dd, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
data.frame(df2_elev_epqs)
##   category elevation elev_units coords.x1 coords.x2 optional
## 1        H    207.81     meters -72.94216  42.04268     TRUE
## 2        M    198.95     meters -72.65599  42.60974     TRUE
## 3        L    306.15     meters -72.91468  42.26707     TRUE
sp_elev_epqs <- get_elev_point(examp_sp, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
sp_elev_epqs
## class       : SpatialPointsDataFrame 
## features    : 3 
## extent      : -72.94216, -72.65599, 42.04268, 42.60974  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 2
## names       : elevation, elev_units 
## min values  :    198.95,     meters 
## max values  :    306.15,     meters
spdf_elev_epqs <- get_elev_point(examp_spdf, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
spdf_elev_epqs
## class       : SpatialPointsDataFrame 
## features    : 3 
## extent      : -72.94216, -72.65599, 42.04268, 42.60974  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 3
## names       : category, elevation, elev_units 
## min values  :        H,    198.95,     meters 
## max values  :        M,    306.15,     meters

Point elevation from Amazon Web Service Terrain Tiles

Since version 0.2.0, elevatr has also provided access to point elevations from the AWS Terrain Tiles. This is not a separate service from the raster DEM service (see below). It is provided as a convenience so users don’t need to download the raster DEMs from AWS and perform their own point data extraction. The added benefit is that points outside of the United States may be used with the AWS source. To access the the point elevations using “aws”:

df_elev_aws <- get_elev_point(examp_df, prj = prj_dd, src = "aws")
## Mosaicing & Projecting
## Note: Elevation units are in meters

An important thing to note, that the elevations will differ, and the prime reason is the resolution of the AWS tiles at the specified zoom. The default zoom of 5 (i.e., z=5) is rather coarse and that is reflected in the elevations.

df_elev_aws$elevation
## [1] 313 221 348
df_elev_epqs$elevation
## [1] 207.81 198.95 306.15

A larger zoom results in a smaller pixel size and the two sources converge.

df_elev_aws_z12 <- get_elev_point(examp_df, prj = prj_dd, src = "aws", z = 12)
## Mosaicing & Projecting
## Note: Elevation units are in meters
df_elev_aws_z12$elevation
## [1] 209 198 307
df_elev_epqs$elevation
## [1] 207.81 198.95 306.15

Determining the correct zoom is a function of the needs of the user and represents a trade off between higher accuracy/longer download.

Lastly, an example using locations outside of the United States.

mt_everest <- data.frame(x = 86.925, y = 27.9881)
everest_aws_elev <- get_elev_point(mt_everest, prj = prj_dd, z = 14, src = "aws")
## Mosaicing & Projecting
## Note: Elevation units are in meters
everest_aws_elev
## class       : SpatialPointsDataFrame 
## features    : 1 
## extent      : 86.925, 86.925, 27.9881, 27.9881  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 2
## names       : elevation, elev_units 
## value       :   8730.61,     meters

Get Raster Elevation Data

While point elevations are useful, they will not provide the information required for most elevation based analysis such as hydrologic modeling, viewsheds, etc. To do that requires a raster digital elevation model (DEM). There are several sources for digital elevation models such as the Shuttle Radar Topography Mission (SRTM), the USGS National Elevation Dataset (NED), Global DEM (GDEM), and others. Each of these DEMs has pros and cons for their use. Prior to its closure in January of 2018, Mapzen combined several of these sources to create a synthesis elevation product that utilizes the best available elevation data for a given region at given zoom level. Additionally, the elevation data are enhanced with the inclusion of bathymetry in oceans from ETOPO1. Although closed, these data compiled by Mapzen are made available through two separate APIs: the Nextzen Terrain Tile Service and the Terrain Tiles on Amazon Web Services. Only the Amazon tiles are currently accessible via elevatr.

The input for get_elev_raster() is a data.frame with x (longitude) and y (latitude) locations as the first two columns, any sp object, or any raster object and it returns a RasterLayer of the tiles that overlap the bounding box of the input. If multiple tiles are retrieved, the resultant output is a merged Raster Layer. Details for each service and their usage via get_elev_raster() are provided below.

Using get_elev_raster() to access the Terrain Tiles on AWS.

As mentioned a data frame with x and y columns, a sp object, or a raster object needs be the input and the src needs to be set to “mapzen” (this is the default).

There is no difference in using the sp and raster input data types. The data frame requires a prj. We show examples using a SpatialPolygonsDataFrame and a data frame. The zoom level (z) defaults to 9 (a trade off between resolution and time for download), but different zoom levels are often desired. For example:

# SpatialPolygonsDataFrame example
data(lake)
elevation <- get_elev_raster(lake, z = 9)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation)
plot(lake, add = TRUE)

# data.frame example
elevation_df <- get_elev_raster(examp_df, prj = prj_dd, z = 5)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation_df)
plot(examp_sp, add = TRUE)

The zoom level determines the resolution of the output raster. More details on resolution and zoom level is still available in the Mapzen Documentation on ground resolution.

In addition the the required arguments (locations, z, and prj for data frames), several additional arguments may be passsed to get_elev_raster(). First, the expand argument is provided to expand the size of the bounding box by a given value in map units. This is useful when bounding box coordinates are near the edge of an xyz tile. For example:

# Bounding box on edge
elev_edge <- get_elev_raster(lake, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elev_edge)
plot(lake, add = TRUE)

# Use expand to grab additional tiles
elev_expand <- get_elev_raster(lake, z = 10, expand = 15000)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elev_expand)
plot(lake, add = TRUE)

Lastly, ... provides the ability to pass additional arguments to httr::GET which is used to access the API endpoints. While any httr::GET arguments may be used, this will most likely be used to pass on configuration arguments such as httr::timeout() or httr::verbose() via a named argument, config to httr::GET. The httr::timeout() can be used to increase the timeout if downloads are timing out. For instance:

library(httr)
# Increase timeout:
get_elev_raster(lake, z = 5, config = timeout(5))
## Mosaicing & Projecting
## Note: Elevation units are in meters.
## class      : RasterLayer 
## dimensions : 640, 597, 382080  (nrow, ncol, ncell)
## resolution : 1731.916, 1731.916  (x, y)
## extent     : 1188272, 2222226, 239538.6, 1347965  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
## source     : file38704629611e.tif 
## names      : file38704629611e 
## values     : -32768, 32767  (min, max)

Lastly, multiple configurations may be passed. Below is an example combining httr::timeout() with httr::verbose().

library(httr)
# Increase timeout:
get_elev_raster(lake, z = 5, config = c(verbose(), timeout(5)))
## Mosaicing & Projecting
## Note: Elevation units are in meters.
## class      : RasterLayer 
## dimensions : 640, 597, 382080  (nrow, ncol, ncell)
## resolution : 1731.916, 1731.916  (x, y)
## extent     : 1188272, 2222226, 239538.6, 1347965  (xmin, xmax, ymin, ymax)
## crs        : +proj=aea +lat_0=40 +lon_0=-96 +lat_1=20 +lat_2=60 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs 
## source     : file3870687f4c17.tif 
## names      : file3870687f4c17 
## values     : -32768, 32767  (min, max)

Access OpenTopography API with get_elev_raster()

As of version 0.3.1, the OpenTopography API (https://portal.opentopography.org/apidocs/#/Public/getGlobalDem) has been available as another source of data from get_elev_raster(). To access this data you need to specify which of the available global DEMs you would like to access. The currently available options are: gl3“,”gl1“,”alos“,”srtm15plus".

Starting in January of 2022, all OpenTopography API requests will require an API key. Version 0.4.3 and greater of elevatr supports this.

To create a key, visit https://portal.opentopography.org/myopentopo. If you do not already have an account, you can create one here: https://portal.opentopography.org/newUser. Once you have your account and are logged into your myOpenTopo Workbench, you can create a new key by scrolling to My Account and selecting “myOpenTopoAuthorization and API Key”. Once there scroll down to the API Key section to create your key.

Within elevatr the API key is expected to be stored as an R environment variable. The easiest way to set this is to use the set_opentopo_key() function. It has a single argument for your API key. Once this is set, restart R and elevatr will use this key for subsequent OpenTopography API requests.

Below is an example for grabbing the OpenTopography SRTM data.

lake_srtmgl1 <- get_elev_raster(lake, src = "gl1", clip = "bbox", expand = 1000)
## Downloading OpenTopography DEMs
## 
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 8.2 kB     
Downloading: 16 kB     
Downloading: 16 kB     
Downloading: 16 kB     
Downloading: 16 kB     
Downloading: 25 kB     
Downloading: 25 kB     
Downloading: 25 kB     
Downloading: 25 kB     
Downloading: 25 kB     
Downloading: 25 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 33 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 41 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 49 kB     
Downloading: 57 kB     
Downloading: 57 kB     
Downloading: 57 kB     
Downloading: 57 kB     
Downloading: 57 kB     
Downloading: 57 kB     
Downloading: 65 kB     
Downloading: 65 kB     
Downloading: 65 kB     
Downloading: 65 kB     
Downloading: 74 kB     
Downloading: 74 kB     
Downloading: 74 kB     
Downloading: 74 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 82 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 90 kB     
Downloading: 98 kB     
Downloading: 98 kB     
Downloading: 98 kB     
Downloading: 98 kB     
Downloading: 98 kB     
Downloading: 98 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 110 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 120 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 130 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 140 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 150 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 160 kB     
Downloading: 170 kB     
Downloading: 170 kB     
Downloading: 170 kB     
Downloading: 170 kB     
Downloading: 180 kB     
Downloading: 180 kB     
Downloading: 180 kB     
Downloading: 180 kB     
Downloading: 180 kB     
Downloading: 180 kB     
Downloading: 190 kB     
Downloading: 190 kB     
Downloading: 190 kB     
Downloading: 190 kB     
Downloading: 190 kB     
Downloading: 190 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 200 kB     
Downloading: 210 kB     
Downloading: 210 kB     
Downloading: 210 kB     
Downloading: 210 kB
## 
## Mosaicing & Projecting
## Clipping DEM to bbox
## Note: Elevation units are in meters.
plot(lake_srtmgl1)
plot(lake, add = TRUE)