This is an R Markdown document containing code for the workshop “Accessing the USGS National Map and making 3D maps with terrainr”, held virtually on 2021-05-28. If you want to follow along with the workshop as we go, click this link to download the R Markdown notebook.

If you’re not familiar with R Markdown, this document lets us write both plain text and code in a single document. Document sections inside three ` marks are called “code chunks”:

plot(Orange)

A single line in a code chunk can be run by pressing Control+Enter with your cursor on the line. The entire chunk can be run by pressing Control+Shift+Enter with your cursor inside the chunk.

Part 1: Installing Packages

For this workshop, we first need to install the packages that we’ll be using. This chunk will install any missing packages, and then load them all:

if (!require("sf")) install.packages("sf")
## Loading required package: sf
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
if (!require("terrainr")) install.packages("terrainr")
## Loading required package: terrainr
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
if (!require("scico")) install.packages("scico")
## Loading required package: scico
if (!require("raster")) install.packages("raster")
## Loading required package: raster
## Loading required package: sp
if (!require("progressr")) install.packages("progressr")
## Loading required package: progressr
if (!require("progress")) install.packages("progress")
## Loading required package: progress
if (!require("tiff")) install.packages('tiff')
## Loading required package: tiff
library("sf")
library("terrainr")
library("ggplot2")
library("scico")
library("raster")
library("progressr")
handlers("progress")

Part 2: Retrieving Data

Up next, let’s download some data! For the purposes of today’s workshop, we’ll work with campsite locations within Bryce Canyon National Park. We can download this data from the National Park Service using the read_sf function from sf:

campsites <- read_sf("https://opendata.arcgis.com/datasets/06993a2f20bf42d382c8ce6bd54027c4_0.geojson")

This dataset can be found online at https://public-nps.opendata.arcgis.com/datasets/bryce-canyon-national-park-backcountry-campsites-?geometry=-112.494%2C37.471%2C-111.908%2C37.566

We now have an sf object containing a handful of campsite locations! We can plot our data using ggplot2 to get a sense of what the spatial distribution of the campsites looks like:

ggplot(campsites) + 
  geom_sf(shape = 4, color = "red") + 
  theme_void()

This is the area we want to download data for from the National Map. Note that everything we do from here on out works with pretty much any sf object; for instance, the terrainr documentation (at https://docs.ropensci.org/terrainr/) often starts with a single latitiude/longitude point and then works from there.

By default, terrainr will download data for the bounding box of whatever sf object you provide it. That is, data will only be downloaded for the smallest rectangle of area that includes your entire sf object.

In this case, we want a little bit of a buffer around our campsite locations so that our points aren’t at the exact edge of the map. To make that happen, we can use set_bbox_side_length to specifically control how much data to download. In this case, a 16 kilometer side length works pretty well:

campsite_bbox <- set_bbox_side_length(campsites, 16, "km")

We can add this to our plot to make sure we’re giving a good buffer to our furthest points:

ggplot() + 
  geom_sf(data = st_set_crs(campsite_bbox, 4326), color = "black", fill = NA) + 
  geom_sf(data = campsites, shape = 4, color = "red") + 
  theme_void()

Since all of our points fall comfortably within this box, we can use it to download our data. We’ll use the get_tiles function to query the National Map to download elevation data (from the USGS 3D Elevation Program DEM) and orthoimagery (from the USDA National Agricultural Imagery Program). This function will break our query into smaller parts and then save the tiles as separate files; we can control where the files get saved by using the output_prefix argument.

We’ll wrap this query in the with_progress function in order to have a progress bar display while our data downloads. While for the purposes of this workshop we’ll be downloading 30 meter rasters, which download pretty quickly (and only require one tile per type of data), this can be very useful for monitoring downloads when using higher-resolution data.

with_progress(
  output_tiles <- get_tiles(campsite_bbox, 
                            output_prefix = "bryce_canyon_np_30m",
                            services = c("elevation", "ortho"),
                            resolution = 30)
)

By assigning the output from get_tiles to output_tiles, we create a list of all of the map tiles we’ve downloaded:

output_tiles
## $elevation
## [1] "bryce_canyon_np_30m_3DEPElevation_1_1.tif"
## 
## $ortho
## [1] "bryce_canyon_np_30m_USGSNAIPPlus_1_1.tif"

Part 3: Plotting Data

This can be really helpful when trying to work with these files inside R. For instance, to see our elevation raster, we can plot our data easily enough using functions from ggplot2. First, we’ll load the raster, then convert it to a data frame:

elevation <- raster(output_tiles$elevation)
elevation <- as.data.frame(elevation, xy = TRUE)
names(elevation) <- c("x", "y", "z")

And then plot the data using geom_raster:

ggplot(elevation, aes(x, y, fill = z)) +
  geom_raster() + 
  scale_fill_gradientn(colors = terrain.colors(255)) + 
  coord_sf(crs = 4326) + 
  theme_void()