---
title: "Accessing elevation data in R with the elevatr package"
author: "Jeffrey W. Hollister"
date: '`r Sys.Date()`'
output:
  html_document:
    theme: readable
    toc: yes
    toc_float: yes
vignette: >
  %\VignetteIndexEntry{Introduction to elevatr} 
  %\VignetteEngine{knitr::rmarkdown} 
  %\VignetteEncoding{UTF-8}
---


```{r setup, include=FALSE, echo=FALSE}
################################################################################
#Load packages
################################################################################
library("sp")
library("raster")
library("knitr")
library("elevatr")
library("httr")
NOT_CRAN <- identical(tolower(Sys.getenv("NOT_CRAN")), "true")
knitr::opts_chunk$set(purl = NOT_CRAN, 
                      eval = NOT_CRAN,
                      fig.width = 4, 
                      fig.height = 4, 
                      tidy = TRUE,
                      dpi = 100)
```

```{r environ, echo=FALSE}
#key <- readRDS("../tests/testthat/key_file.rds")
#Sys.setenv(mapzen_key=key)
```

# 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.

There are currently two endpoints that `elevatr` accesses.  For point elevation data it uses USGS Elevation Point Query Service and to access raster elevation data (e.g. a DEM) it uses the Amazon Web Services Terrain Tiles.

# 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](http://ned.usgs.gov/epqs/) 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.

```{r example_dataframe}
# Create an example data.frame
set.seed(65.7)
examp_df <- data.frame(x = runif(10, min = -73, max = -71), 
                       y = runif(10, min = 41 , max = 45))
prj_dd <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"

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

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, proj4string = CRS(prj_dd), 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.


```{r}
df_elev_epqs <- get_elev_point(examp_df, prj = prj_dd, src = "epqs")
data.frame(df_elev_epqs)
df2_elev_epqs <- get_elev_point(examp_df2, prj = prj_dd, src = "epqs")
data.frame(df2_elev_epqs)
sp_elev_epqs <- get_elev_point(examp_sp, src = "epqs")
sp_elev_epqs
spdf_elev_epqs <- get_elev_point(examp_spdf, src = "epqs")
spdf_elev_epqs
```

# 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](https://mapzen.com/documentation/terrain-tiles/data-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](https://www.nextzen.org#terrain-tiles) and the [Terrain Tiles on Amazon Web Services](https://aws.amazon.com/public-datasets/terrain/).  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:  

```{r get_raster}
# SpatialPolygonsDataFrame example
data(lake)
elevation <- get_elev_raster(lake,z = 9)
plot(elevation)
plot(lake, add=TRUE)

# data.frame example
elevation_df <- get_elev_raster(examp_df,prj=prj_dd, z = 5)
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](https://mapzen.com/documentation/terrain-tiles/data-sources/#what-is-the-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:

```{r expand}
# Bounding box on edge
elev_edge<-get_elev_raster(lake, z = 10)
plot(elev_edge)
plot(lake, add = TRUE)

# Use expand to grab additional tiles
elev_expand<-get_elev_raster(lake, z = 10, expand = 1500)
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:

```{r timeout}
# Increase timeout:
get_elev_raster(lake, z = 5, config = timeout(5))
```

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

```{r timeout_verbose}
# Increase timeout:
get_elev_raster(lake, z = 5, config = c(verbose(),timeout(5)))
```
