A gentle introduction to terrainr

The goal of the terrainr package is to make it easier to visualize landscapes, both by providing functions to download elevation data and base maps from the USGS National Map and by adding utilities to manipulate base maps for use in visualizations using ggplot2 or the freely available Unity 3D rendering engine. This vignette will walk through the core functions available in terrainr and how they interact.

Let’s load the terrainr package to get started:

library(terrainr)

We’re going to work with data for Mount Elbert today, the highest point in the Rocky Mountain range. I’m just choosing this location for the dramatic scenery; the National Map can be used to retrieve data for the entire United States and much of Canada. Let’s simulate some data for the area right around Mt. Elbert, such as the point data we might get from some field collection:

mt_elbert_points <- data.frame(
  lat = runif(100, min = 39.11144, max = 39.12416),
  lng = runif(100, min = -106.4534, max = -106.437)
)

terrainr is built to play nicely with functions from the sf and raster packages. In order to get our simulated points into the right format, we need to use the st_as_sf function from the sf package:

mt_elbert_points <- sf::st_as_sf(mt_elbert_points,
                                 coords = c("lng", "lat"))
mt_elbert_points <- sf::st_set_crs(mt_elbert_points, 4326)

Now that we’ve got our data in the right format, it’s time to retrieve our data. terrainr currently supports downloading DEMs from the USGS 3D Elevation Program as well as orthoimages from the National Agricultural Imagery Program, in addition to other base map images from the National Map. These programs each have slightly different APIs and different restrictions on file types and the size of image you can download at once. Rather than make you think about this, terrainr handles all the edges of making API requests for you, including splitting your request into tiles and formatting the query.

For this vignette, we’ll retrieve both elevation and orthoimagery using the get_tiles function. We can either use the generic “elevation” and “ortho” shorthands to get our data, or we can specify “3DEPElevation” and “USGSNAIPPlus” to make sure we’re using the same specific service – the short codes aren’t guaranteed to download data from the same service between releases!

One last note – all the longer-running terrainr functions can print out progress bars, if the user requests them via the progressr package. We’ll demonstrate that syntax here:

library(progressr)
handlers("progress")
with_progress(
  output_files <- get_tiles(mt_elbert_points,
                            output_prefix = tempfile(),
                            services = c("elevation", "ortho"))
  )

And just like that, we have our data tiles! To make multi-step processing easier, terrainr functions which deal with these tiles typically return lists of the file paths they saved your data to.

output_files
#> $elevation
#> [1] "/tmp/RtmphTFQvZ/file65e5d859e628_3DEPElevation_1_1.tif"
#> 
#> $ortho
#> [1] "/tmp/RtmphTFQvZ/file65e5d859e628_USGSNAIPPlus_1_1.tif"

If we were requesting more data than we can download at once, each element of the list would be a character vector containing the file paths for all of our downloaded tiles. Since we’re sticking with a relatively small area for this example, we only have one tile for each service.

As a quick aside, note that you can control where these files save to via the output_prefix argument (which appends the suffix servicename_xindex_yindex.tif to each tile it downloads) – you don’t need to save them to a temporary directory (and redownload every time you launch R) as we’re doing here!

If all you want is to access these endpoints to download data, this is probably the only terrainr function you’ll need – the files produced by this function can be processed just like any other spatial data:

raster::plot(raster::raster(output_files[[1]]))

raster::plotRGB(raster::brick(output_files[[2]]), scale = 1)

In addition to the regular methods for plotting rasters in R, terrainr makes it a bit easier to use ggplot2 for plotting the data returned by get_tiles. Plotting single-band rasters, like our elevation file, is already well-supported in base ggplot2:

library(ggplot2)

elevation_raster <- raster::raster(output_files[[1]])
elevation_df <- as.data.frame(elevation_raster, xy = TRUE)
elevation_df <- setNames(elevation_df, c("x", "y", "elevation"))

ggplot() + 
  geom_raster(data = elevation_df, aes(x = x, y = y, fill = elevation)) + 
  scale_fill_distiller(palette = "BrBG") + 
  coord_sf(crs = 4326)

terrainr adds the ability to plot using multi-band RGB rasters, like the tiles downloaded for non-elevation endpoints, using the new geom_spatial_rgb function (or its partner, stat_spatial_rgb):

ortho_raster <- raster::stack(output_files[[2]])
ortho_df <- as.data.frame(ortho_raster, xy = TRUE)
ortho_df <- setNames(ortho_df, c("x", "y", "red", "green", "blue", "alpha"))

ggplot() + 
  geom_spatial_rgb(data = ortho_df,
                   # Required aesthetics r/g/b specify color bands:
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  coord_sf(crs = 4326)
knitr::include_graphics("ortho_ggplot.jpg")

Note that geom_spatial_rgb is a little different from other ggplot2 geoms in that it can also accept RasterStack objects directly:

ggplot() + 
  geom_spatial_rgb(data = ortho_raster,
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  coord_sf(crs = 4326)

Or length 1 character vectors with a path to a file that can be read by raster::stack:

ggplot() + 
  geom_spatial_rgb(data = output_files[[2]],
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  coord_sf(crs = 4326)

You can then use these multi-band rasters as base maps for further plotting as desired.

ggplot() + 
  geom_spatial_rgb(data = output_files[[2]],
                   aes(x = x, y = y, r = red, g = green, b = blue)) + 
  geom_sf(data = mt_elbert_points)

In case you find this visualization falls a little bit flat, terrainr also provides the ability to bring your landscapes into Unity to visualize in 3D. Our first step in this process is going to be replicating our image of base-map-plus-field-sites (above) in a format that we can import into Unity directly. First, we’ll need to use the vector_to_overlay function to create an image overlay from our point data:

mt_elbert_overlay <- vector_to_overlay(mt_elbert_points,
                                       output_files[[2]],
                                       size = 15,
                                       color = "red")
knitr::include_graphics(mt_elbert_overlay)

Note that vector_to_overlay can be used with any sf object, not just point data.

These overlays may be stacked on top of one another or downloaded imagery using the combine_overlays function:

ortho_with_points <- combine_overlays(
  # Overlays are stacked in order, with the first file specified on the bottom
  output_files[[2]],
  mt_elbert_overlay,
  output_file = tempfile(fileext = ".png")
  )
knitr::include_graphics(ortho_with_points)

Unfortunately, this image processing strips the georeferencing on the image. We can restore the original georeferencing via the georeference_overlay function:

georef_overlay <- georeference_overlay(
  ortho_with_points,
  output_files[[2]]
)

We’ve been working so far with a single tile, but Unity is able to handle much, much larger rasters than we would normally work with in R. In order to create overlays for these larger rasters, it’s usually best to create an overlay for smaller image tiles which can then be joined back together with merge_rasters:

tile_overlays <- lapply(output_files[[2]],
                        function(x) vector_to_overlay(mt_elbert_points, 
                                                      x, 
                                                      size = 15, 
                                                      color = "red", 
                                                      na.rm = TRUE))

combined_tiles <- mapply(function(x, y) {
  combine_overlays(x, y, output_file = tempfile(fileext = ".png"))
  },                            
  output_files[[2]],
  tile_overlays)

georef_tiles <- mapply(georeference_overlay, combined_tiles, output_files[[2]])

merged_tiles <- merge_rasters(georef_tiles)

Of course, since we’re only working with a single tile, georef_tiles is identical to merged_tiles. But when working with larger areas, merged_tiles is particularly useful for joining the separate tiles downloaded by get_tiles into a single raster file.

In particular, having a single joined raster is necessary for the function make_manifest, which is designed to turn these larger rasters into tiles in a format that can be imported into the Unity 3D rendering engine. You can find more information about that process in the Unity vignette.

elevation_tile <- output_files[[1]]
make_manifest(elevation_tile, georef_tiles)

After that function runs, it’s a matter of minutes to create beautiful – and fully physically-simulated – landscape visualizations of your area of interest in Unity:

For more instructions on how to create these 3D simulations in Unity, check out the Unity vignette.