elevatr
elevatr
get_elev_point()
to Access The USGS Elevation
Point Query Serviceget_elev_raster()
to access the Terrain Tiles on
AWS.get_elev_raster()
elevatr
Several major changes have been made to elevatr
in
response to the retirement of legacy spatial packages (see https://r-spatial.org/r/2023/05/15/evolution4.html for
details). Version 0.99.0 has switched to using sf
and
terra
for all data handling; however, in this version a
raster RasterLayer
is still returned from
get_elev_raster()
. Additional changes are planned for
version 1+, most notably the return for get_elev_raster()
will be a terra SpatRaster
. Please plan accordingly for
your analyses and/or packages account for this change.
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.
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 sf
, 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.
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.
get_elev_point()
to Access The USGS Elevation
Point Query ServiceUsage 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))
crs_dd <- 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
examp_sf <- sf::st_as_sf(examp_df2, coords = c("x", "y"), crs = crs_dd)
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
POINT
or MULTIPOINT
sf
object.
Similarly, an elevation and units column is added to the data frame.
The USGS Elevation Point Query Service returns a single point at a
time. The implementation 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 = crs_dd, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
df_elev_epqs
## Simple feature collection with 3 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -72.94216 ymin: 42.04268 xmax: -72.65599 ymax: 42.60974
## Geodetic CRS: WGS 84
## geometry elevation elev_units
## 1 POINT (-72.94216 42.04268) 208.44 meters
## 2 POINT (-72.65599 42.60974) 199.09 meters
## 3 POINT (-72.91468 42.26707) 306.40 meters
df2_elev_epqs <- get_elev_point(examp_df2, prj = crs_dd, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
df2_elev_epqs
## Simple feature collection with 3 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -72.94216 ymin: 42.04268 xmax: -72.65599 ymax: 42.60974
## Geodetic CRS: WGS 84
## category geometry elevation elev_units
## 1 H POINT (-72.94216 42.04268) 208.44 meters
## 2 M POINT (-72.65599 42.60974) 199.09 meters
## 3 L POINT (-72.91468 42.26707) 306.40 meters
sf_elev_epqs <- get_elev_point(examp_sf, src = "epqs")
## Downloading point elevations:
## Note: Elevation units are in meters
sf_elev_epqs
## Simple feature collection with 3 features and 3 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -72.94216 ymin: 42.04268 xmax: -72.65599 ymax: 42.60974
## Geodetic CRS: WGS 84
## category geometry elevation elev_units
## 1 H POINT (-72.94216 42.04268) 208.44 meters
## 2 M POINT (-72.65599 42.60974) 199.09 meters
## 3 L POINT (-72.91468 42.26707) 306.40 meters
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 = crs_dd, src = "aws")
## Mosaicing & Projecting
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
## (status 2 uses the sf package in place of rgdal)
## 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] 208.44 199.09 306.40
A larger zoom results in a smaller pixel size and the two sources converge.
df_elev_aws_z12 <- get_elev_point(examp_df, prj = crs_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] 208.44 199.09 306.40
Determining the correct zoom is a function of the needs of the user and represents a trade off between higher accuracy/longer downloads.
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 = crs_dd, z = 14, src = "aws")
## Mosaicing & Projecting
## Note: Elevation units are in meters
everest_aws_elev
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 86.925 ymin: 27.9881 xmax: 86.925 ymax: 27.9881
## Geodetic CRS: WGS 84
## geometry elevation elev_units
## 1 POINT (86.925 27.9881) 8730.61 meters
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, any sf
object, any
terra
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 RasterLayer
. Details for each
service and their usage via get_elev_raster()
are provided
below.
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:
# sf POLYGON example
data(lake)
elevation <- get_elev_raster(lake, z = 9)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation)
plot(st_geometry(lake), add = TRUE, col = "blue")
# data.frame example
elevation_df <- get_elev_raster(examp_df, prj = crs_dd, z = 5)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
plot(elevation_df)
plot(examp_sf, add = TRUE, col = "black", pch = 19, max.plot = 1)
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 passed 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(st_geometry(lake), add = TRUE, col = "blue")
# 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(st_geometry(lake), add = TRUE, col = "blue")
Second, the clip
argument provides some control over the
extent and shape of the elevation raster that is returned. The default
value returns the entire tile for the “aws” src
. The
default value for the OpenTopography is also tile, but as these datasets
are not served by tile, the “tile” and “bbox” clip return the same
elevation raster. The “locations” clip option will clip the elevation
raster to the locations themselves. If the input locations are points,
this is no different that “bbox”, however if the input locations are a
polygon, the elevation raster will be clipped to the boudary of that
polygon. For example:
lake_buffer <- st_buffer(lake, 1000)
lake_buffer_elev <- get_elev_raster(lake_buffer, z = 9, clip = "locations")
## Mosaicing & Projecting
## Clipping DEM to locations
## Note: Elevation units are in meters.
plot(lake_buffer_elev)
plot(st_geometry(lake), add = TRUE, col = "blue")
plot(st_geometry(lake_buffer), 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(100))
## 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 : file60187ccf4d3a.tif
## names : file60187ccf4d3a
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 : file60182aac318b.tif
## names : file60182aac318b
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)
plot(lake_srtmgl1)
plot(st_geometry(lake), add = TRUE, col = "blue")