elevatrSeveral 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.
elevatrElevation 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 metersdf_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 metersdf2_elev_epqs <- get_elev_point(examp_df2, prj = crs_dd, src = "epqs")## Downloading point elevations:
## Note: Elevation units are in metersdf2_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 meterssf_elev_epqs <- get_elev_point(examp_sf, src = "epqs")## Downloading point elevations:
## Note: Elevation units are in meterssf_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 metersSince 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 metersAn 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 348df_elev_epqs$elevation## [1] 208.44 199.09 306.40A 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 metersdf_elev_aws_z12$elevation## [1] 209 198 307df_elev_epqs$elevation## [1] 208.44 199.09 306.40Determining 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 meterseverest_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 metersWhile 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")