Our fifth notebook: Elevation data in R

1. Introduction

This is Our fifth notebook: Elevation data in R, written using R Markdown Notebook. It illustrates several functionalities to obtain, process and visualize digital elevation models in R

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.

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. This shapefile is part of the marco geostadistico. It is located within a local folder named 70_SUCRE.

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

# SpatialPolygonsDataFrame example
munic <-  read_sf("/Users/Acer/Documents/Materias U/Geomatica Basica/70_SUCRE/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"); select(munic, MPIO_CCDGO, MPIO_CNMBR,)

Let’s download the elevation data for Sucre.

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]]]]

he 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:

# Created a file .tif in the specifed path
# write to geotiff file (depends on rgdal)
if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("/Users/Acer/Documents/Materias U/Geomatica Basica/sucre_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.

elevation

## class      : RasterLayer 
## dimensions : 4077, 2581, 10522737  (nrow, ncol, ncell)
## resolution : 0.000681119, 0.000681119  (x, y)
## extent     : -75.9375, -74.17953, 7.71089, 10.48781  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Acer/AppData/Local/Temp/RtmpqSlkOC/file199465a343e4.tif 
## names      : file199465a343e4 
## 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.

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 : 2742, 1722, 4721724  (nrow, ncol, ncell)
## resolution : 0.000681119, 0.000681119  (x, y)
## extent     : -75.70592, -74.53303, 8.277581, 10.14521  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : file199465a343e4 
## values     : -952, 822  (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 : 1148, 718, 824264  (nrow, ncol, ncell)
## resolution : 180, 180  (x, y)
## extent     : 422259.9, 551499.9, 914911.3, 1121551  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file199465a343e4 
## values     : -341.1644, 804.363  (min, max)

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

(rep_munic = spTransform(munic_sp,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 26 
## extent      : 422484.7, 551436.6, 915025, 1121539  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                 MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,      Shape_Area 
## min values  :         70,      70001, BUENAVISTA,                       1799,   56.55589823,      2017,      SUCRE, 0.426884228585, 0.0046519224476 
## max values  :         70,      70823, TOLÚ VIEJO, Ordenanzas 27 y 42 de 1938, 1494.16295608,      2017,      SUCRE,   2.9849722936,   0.12277744607

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="/Users/Acer/Documents/Materias U/Geomatica Basica/rep_boyaca_dem.tif", datatype='INT4S', overwrite=TRUE)

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

7. 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 Sucre 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 Sucre 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 Sucre 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)

8. 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("COVEÑAS", "PALMITO", "SAN ONOFRE", "SANTIAGO DE TOLÚ", "SINCELEJO")
(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      : 422484.7, 463739.2, 1019785, 1121539  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                  MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,      Shape_Area 
## min values  :         70,      70001,    COVEÑAS,                        1840,   56.55589823,      2017,      SUCRE, 0.426884228585, 0.0046519224476 
## max values  :         70,      70820,  SINCELEJO, Ley 47 de Agosto 30 de 1966, 1035.07793303,      2017,      SUCRE,  2.08052417308, 0.0852548821666

#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-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18363)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.1252  LC_CTYPE=Spanish_Colombia.1252   
## [3] LC_MONETARY=Spanish_Colombia.1252 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.1252    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] rayshader_0.24.5    mapview_2.9.0       forcats_0.5.1      
##  [4] stringr_1.4.0       dplyr_1.0.5         purrr_0.3.4        
##  [7] readr_1.4.0         tidyr_1.1.3         tibble_3.1.0       
## [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.50.1    latticeExtra_0.6-29 lattice_0.20-41    
## [19] terra_1.1-4         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.10       
##  [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             utf8_1.1.4             
## [13] R6_2.5.0                KernSmooth_2.23-18      DBI_1.1.1              
## [16] colorspace_2.0-0        manipulateWidget_0.10.1 withr_2.4.1            
## [19] prettyunits_1.1.1       tidyselect_1.1.0        leaflet_2.0.4.1        
## [22] curl_4.3                compiler_4.0.4          leafem_0.1.3           
## [25] cli_2.3.1               rvest_1.0.0             xml2_1.3.2             
## [28] sass_0.3.1              scales_1.1.1            classInt_0.4-3         
## [31] hexbin_1.28.2           digest_0.6.27           rmarkdown_2.7          
## [34] base64enc_0.1-3         jpeg_0.1-8.1            pkgconfig_2.0.3        
## [37] htmltools_0.5.1.1       highr_0.8               dbplyr_2.1.0           
## [40] fastmap_1.1.0           htmlwidgets_1.5.3       rlang_0.4.10           
## [43] readxl_1.3.1            rstudioapi_0.13         shiny_1.6.0            
## [46] jquerylib_0.1.3         generics_0.1.0          zoo_1.8-9              
## [49] jsonlite_1.7.2          crosstalk_1.1.1         rayimage_0.5.1         
## [52] magrittr_2.0.1          Rcpp_1.0.6              munsell_0.5.0          
## [55] fansi_0.4.2             lifecycle_1.0.0         stringi_1.5.3          
## [58] yaml_2.2.1              grid_4.0.4              parallel_4.0.4         
## [61] promises_1.2.0.1        crayon_1.4.1            miniUI_0.1.1.1         
## [64] haven_2.3.1             hms_1.0.0               knitr_1.31             
## [67] pillar_1.5.1            stats4_4.0.4            codetools_0.2-18       
## [70] reprex_1.0.0            glue_1.4.2              evaluate_0.14          
## [73] modelr_0.1.8            foreach_1.5.1           png_0.1-7              
## [76] vctrs_0.3.6             httpuv_1.5.5            cellranger_1.1.0       
## [79] gtable_0.3.0            assertthat_0.2.1        xfun_0.21              
## [82] mime_0.10               xtable_1.8-4            broom_0.7.5            
## [85] e1071_1.7-4             later_1.1.0.1           class_7.3-18           
## [88] viridisLite_0.3.0       iterators_1.0.13        units_0.7-0            
## [91] ellipsis_0.3.1