My fifth notebook: Elevation data in R

1. Why this notebook?

This is an R Notebook created using R Studio on an old laptop. It illustrates several functionalities to obtain, process and visualize digital elevation models in R. It aims to help Geomatica Basica students at Universidad Nacional to get started with R geospatial capabilities.

A few tips for writing your own notebook:

• Write and run every chunk step-by-step
• In case of errors, read the message, interpret it and try to fix it
• In case the error persists, write a message on slack (be as precise and descriptive as possible)
• Once your html file is produced, it is time to publish it at RPubs
• After publishing the notebook, review it, fix poor quality wording, and republish it using the update option

In case of trouble, please review the rmarkdown documentation.

Do not forget to look for help at community.rstudio.com/.

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

To access raster elevation data (e.g. a DEM) the elevatr package uses the Amazon Web Services Terrain Tiles. We will explore this functionality in this notebook.

#install the packages from the console not from here
# install.packages("rgdal")
# install.packages("raster")
# install.packages("elevatr")
#install.packages("rasterVis")
#install.packages("rgl")
#install.packages("mapview")

library(rasterVis)
library(raster)
library(rgl)
library(rgdal)
library(elevatr)
library(sf)
library(tidyverse)
library(mapview)

3. Get Raster Elevation Data

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. Although closed, these data compiled by Mapzen are currently accessible through the Terrain Tiles on Amazon Web Services (AWS).

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.

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 CRS (i.e. a projection file). I will show examples using a simple feature. The zoom level (z) defaults to 9 (a trade off between resolution and time for download). I was able to get elevation data at zoom = 10.

First, we need to load a shapefile representing municipalities of our department. In this notebook, we will load a shapefile downloaded from the DANE’s geoportal as illustrated last week in the GB class. This shapefile is part of the marco geostadistico. It is located within a local folder named 15_BOYACA.

Let’s read the shapefile using a function provided by the sf package:

# SpatialPolygonsDataFrame example
(munic <-  st_read("../../datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

## Reading layer `MGN_MPIO_POLITICO' from data source `/Users/ivanlizarazo/Documents/ivan/GB/datos/shapefiles/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 123 features and 2 fields
## geometry type:  POLYGON
## dimension:      XY
## bbox:           xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## geographic CRS: WGS 84

What are the attributes of the munic object:

head(munic)

Let’s download the elevation data for Boyaca.

elevation <- get_elev_raster(munic, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  GEOGCRS["WGS 84 (with axis order normalized for visualization)",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]

The messaged displayed by RStudio should look similar to this one:

The get_elev_raster function downloads elevation data as a raster (i.e. as a raster DEM) from the AWS Open Data Terrain Tiles. This function accepts a data frame of of x (long) and y (lat) or any sp or raster object as input and will return a raster object of the elevation tiles that cover the bounding box of the input spatial data.

Note that elevation is a temporal raster (i.e. it exists only in-memory).However, it is possible to write the DEM to your local directory using the writeRaster function provided by the rgdal library:

# uncomment if needed 
# write to geotiff file (depends on rgdal)
if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("../../datos/boyaca_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}

In case you have already downloaded and saved the elevation data in a previous session, you have to read the TIF file using the raster function provided by the raster package. You only should run the following chunk in case you have previously saved your DEM.

## Run this code in case you need:
#elevation <-  raster("../../datos/boyaca_dem_z10.tif")

Now, let’s check what is inside the elevation object:

elevation

## class      : RasterLayer 
## dimensions : 4602, 5141, 23658882  (nrow, ncol, ncell)
## resolution : 0.000683952, 0.000683952  (x, y)
## extent     : -74.88281, -71.36662, 4.21492, 7.362467  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : /private/var/folders/l_/1nv5l0dx78g_d3k3_t85lrjr0000gn/T/RtmpyuOnGe/file4bb822f30b96.tif 
## names      : file4bb822f30b96 
## values     : -32768, 32767  (min, max)

Notice a few things about this DEM:

• dimensions: the “size” of the file in pixels
• nrow, ncol: the number of rows and columns in the data (imagine a spreadsheet or a matrix).
• ncells: the total number of pixels or cells that make up the raster.
• resolution: the size of each pixel (in decimal degrees in this case). Let’s remind that 1 decimal degree represents about 111.11 km at the Equator.
• extent: the spatial extent of the raster. This value will be in the same coordinate units as the coordinate reference system of the raster.
• crs: the coordinate reference system string for the raster. This raster is in geographic coordinates with a datum of WGS 84.

munic_sp <- as(munic, "Spatial")

plot(elevation, main="This the downloaded DEM [meters]")
plot(munic_sp,col="NA",border="black", add=TRUE)
text(coordinates(munic_sp), labels=as.character(munic$MPIO_CNMBR), 
     col="black", cex=0.20)

4. Crop the elevation data to match the study area extent

Notice that the DEM covers a larger extent than the one we need. This is due to the spatial arrangement of the elevation tiles in AWS.

Anyway, it is a good idea to save the DEM for the future.

# Write the RasterLayer to disk (See datatype documentation for other formats)
# Review the meaning of datatype in the following link
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/dataType
# The filename includes a prefix related to the zoom level (equal to 8 in this case)
# Uncomment and run the following line:
#writeRaster(elevation, filename="./dem/elev_z8.tif", datatype='INT4S', overwrite=TRUE)

Now, let’s crop the elevation data corresponding to the selected municipalities:

# review the raster documentation to understand how the crop function works
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/crop
elev_crop = crop(elevation, munic_sp)
plot(elev_crop, main="Cropped digital elevation model")
plot(munic_sp, add=TRUE)
text(coordinates(munic_sp), labels=as.character(munic_sp$MPIO_CNMBR), cex=0.2)

Let’s check the new object:

elev_crop

## class      : RasterLayer 
## dimensions : 3509, 3971, 13934239  (nrow, ncol, ncell)
## resolution : 0.000683952, 0.000683952  (x, y)
## extent     : -74.66463, -71.94866, 4.655385, 7.055372  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : file4bb822f30b96 
## values     : -270, 5331  (min, max)

5. Reproject the elevation data

When working with DEMs, it is always a good idea to use map coordinates rather than geographic coordinates. This is due to the fact that, in geographic coordinates, horizontal dimensions units are decimal degrees BUT the vertical dimension unit is meters. Let’s reproject the elevation data.

crs(elev_crop)

## CRS arguments: +proj=longlat +datum=WGS84 +no_defs

library(sp)
spatialref <- CRS(SRS_string="EPSG:32618")
##spatialref <-  CRS(SRS_string = "EPSG:9377")
##spatialref <- st_crs(9377)

Now, we can reproject the elevation data from WGS84 geographic coordinates into the UTM 18N CRS.

# Project Raster
# to have control,  we create one raster object
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elev_crop, crs=spatialref)
# Adjust the cell size 
res(pr3) <- 180
# now project
rep_elev <- projectRaster(elev_crop, pr3)

What we got?

rep_elev

## class      : RasterLayer 
## dimensions : 1480, 1675, 2479000  (nrow, ncol, ncell)
## resolution : 180, 180  (x, y)
## extent     : 537037.6, 838537.6, 514573.6, 780973.6  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file4bb822f30b96 
## values     : -36.43087, 5311.972  (min, max)

Now, let’s reproject the SpatialPolygonsDataFrame representing the municipalities of our department:

(rep_munic = spTransform(munic_sp,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : 537094.2, 837149.8, 514819.8, 780834.7  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 2
## names       : MPIO_CCDGO, MPIO_CNMBR 
## min values  :      15001,    ALMEIDA 
## max values  :      15897,  ZETAQUIRÁ

It is plotting time:

plot(rep_elev, main="Reprojected digital elevation model")
plot(rep_munic, lwd=0.5, add=TRUE)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

To avoid a headache, let’s save our DEM:

# Write the reprojected DEM to disk 
#uncomment and run the following line
#once you run the line, comment it again
writeRaster(rep_elev, filename="../../datos/rep_boyaca_dem.tif", datatype='INT4S', overwrite=TRUE)

5. Basic statistics of elevation data

A quick exploration of the DEM statistics:

# histograma
hist(rep_elev)

# works for large files
promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion  <- cellStats(rep_elev, 'sd')

# create unidimensional vector
metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)

# creacion de data frame con estadisticas de elevacion [metros]
# el paréntesis externo sirve para "imprimir" el contenido de un objeto
(df_estadisticas <- data.frame(metricas, valores))

Interpretation of DEM statistics is helpful in many applications.

6. Obtention of geomorphometric variables

Before proceeding, make sure you have read the geomorphometry chapter written by Victor Olaya for his libro libre en sistemas de informacion geografica.

First, compute slope, aspect, and hillshade:

slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)

Plot elevation. Note the color palette used here.

plot(rep_elev,main="DEM for Boyaca municipalities [meters]", col=terrain.colors(25,alpha=0.8))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

Plot slope. Note the color palette used here.

plot(slope,main="Slope for Boyaca municipalities [degrees]", col=topo.colors(6,alpha=0.6))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

Plot aspect. Note the color palette used here.

plot(aspect,main="Aspect for several municipalities [degrees]", col=rainbow(10,alpha=0.7))
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

In case you are interested, read about color palettes in R.

A combined plot:

# plot DEM w/ hillshade
plot(hill,
        col=grey(1:100/100),  # create a color ramp of grey colors for hillshade
        legend=FALSE,         # no legend, we don't care about the grey of the hillshade
        main="DEM for Boyaca municipalities [m]",
        axes=FALSE)           # makes for a cleaner plot, if the coordinates aren't necessary

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.6), add=TRUE) # color method
         # sets how transparent the object will be (0=transparent, 1=not transparent)
plot(rep_munic, add=TRUE, lwd=0.5)
text(coordinates(rep_munic), labels=as.character(rep_munic$MPIO_CNMBR), cex=0.2)

7. Mapping elevation data with rayshader

The library rayshader is an open source package for producing 2D and 3D data visualizations in R. rayshader uses elevation data in a base R matrix and a combination of raytracing, spherical texture mapping, overlays, and ambient occlusion to generate beautiful topographic 2D and 3D maps. In addition to maps, rayshader also allows the user to translate ggplot2 objects into beautiful 3D data visualizations.

#install.packages("rayshader")

library(rayshader)

#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)

#And we can add a raytraced layer from that sun direction as well:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_map()

It may be helpful to attempt a visualization for a smaller region:

lista <- list("PESCA", "CUÌTIVA", "TOTA", "IZA", "AQUITANIA")
(algunos <- munic %>% filter(MPIO_CNMBR %in% unlist(lista) ))

# this chunk only run interactively
#algunos %>% mapview::mapview(zcol = "MPIO_CNMBR", legend=TRUE, selfcontained = FALSE)

algunos_sp <- as(algunos, "Spatial")

elev_crop = crop(elevation, algunos_sp)
plot(elev_crop, main="Cropped digital elevation model")
plot(algunos_sp, add=TRUE)
text(coordinates(algunos_sp), labels=as.character(algunos_sp$MPIO_CNMBR), cex=0.5)

# Project Raster
# to have control,  we create one raster object
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elev_crop, crs=spatialref)
# Adjust the cell size 
res(pr3) <- 120
# now project
rep_elev <- projectRaster(elev_crop, pr3)

(rep_algunos = spTransform(algunos_sp,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 5 
## extent      : 705258.5, 760124.8, 574692, 624808.9  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 2
## names       : MPIO_CCDGO, MPIO_CNMBR 
## min values  :      15047,  AQUITANIA 
## max values  :      15822,       TOTA

#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)

#detect_water and add_water adds a water layer to the map:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  add_shadow(ambient_shade(elmat), 0) %>%
  plot_map()

Final notes: This notebook illustrated a few DEM visualization options. You can try your own ideas.

The rayshader documentation is here.

A guide to RGL visualization is at _http://www.sthda.com/english/wiki/a-complete-guide-to-3d-visualization-device-system-in-r-r-software-and-data-visualization_.

In case of running into trouble, please review the documentation. Look for tutoriales and/or examples, write your own code, and you will learn a lot.

Good luck!

9. Reproducibility

sessionInfo()

## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rayshader_0.19.2    mapview_2.9.0       forcats_0.5.1      
##  [4] stringr_1.4.0       dplyr_1.0.4         purrr_0.3.4        
##  [7] readr_1.4.0         tidyr_1.1.2         tibble_3.0.6       
## [10] ggplot2_3.3.3       tidyverse_1.3.0     sf_0.9-7           
## [13] elevatr_0.3.4       rgdal_1.5-23        rgl_0.105.22       
## [16] rasterVis_0.49      latticeExtra_0.6-29 lattice_0.20-41    
## [19] raster_3.4-5        sp_1.4-5           
## 
## loaded via a namespace (and not attached):
##  [1] fs_1.5.0                satellite_1.0.2         lubridate_1.7.9.2      
##  [4] doParallel_1.0.16       progress_1.2.2          webshot_0.5.2          
##  [7] RColorBrewer_1.1-2      httr_1.4.2              tools_4.0.4            
## [10] backports_1.2.1         bslib_0.2.4             R6_2.5.0               
## [13] KernSmooth_2.23-18      DBI_1.1.1               colorspace_2.0-0       
## [16] manipulateWidget_0.10.1 withr_2.4.1             prettyunits_1.1.1      
## [19] tidyselect_1.1.0        leaflet_2.0.4.1         curl_4.3               
## [22] compiler_4.0.4          leafem_0.1.3            cli_2.3.0              
## [25] rvest_0.3.6             xml2_1.3.2              sass_0.3.1             
## [28] scales_1.1.1            classInt_0.4-3          hexbin_1.28.2          
## [31] digest_0.6.27           rmarkdown_2.7           base64enc_0.1-3        
## [34] jpeg_0.1-8.1            pkgconfig_2.0.3         htmltools_0.5.1.1      
## [37] highr_0.8               dbplyr_2.1.0            fastmap_1.1.0          
## [40] htmlwidgets_1.5.3       rlang_0.4.10            readxl_1.3.1           
## [43] rstudioapi_0.13         shiny_1.6.0             jquerylib_0.1.3        
## [46] generics_0.1.0          zoo_1.8-8               jsonlite_1.7.2         
## [49] crosstalk_1.1.1         rayimage_0.5.1          magrittr_2.0.1         
## [52] Rcpp_1.0.6              munsell_0.5.0           lifecycle_1.0.0        
## [55] stringi_1.5.3           yaml_2.2.1              grid_4.0.4             
## [58] parallel_4.0.4          promises_1.2.0.1        crayon_1.4.1           
## [61] miniUI_0.1.1.1          haven_2.3.1             hms_1.0.0              
## [64] knitr_1.31              pillar_1.4.7            stats4_4.0.4           
## [67] codetools_0.2-18        reprex_1.0.0            glue_1.4.2             
## [70] evaluate_0.14           modelr_0.1.8            foreach_1.5.1          
## [73] png_0.1-7               vctrs_0.3.6             httpuv_1.5.5           
## [76] cellranger_1.1.0        gtable_0.3.0            assertthat_0.2.1       
## [79] xfun_0.21               mime_0.10               xtable_1.8-4           
## [82] broom_0.7.5             e1071_1.7-4             later_1.1.0.1          
## [85] class_7.3-18            viridisLite_0.3.0       iterators_1.0.13       
## [88] units_0.6-7             ellipsis_0.3.1