1. Why this notebook?

This notebook illustrates several functionalities to obtain, process and visualize digital elevation models (DEMs) in R. It aims to help Geomatica Basica students at Universidad Nacional to get acquainted with DEMs and geomorphometric variables.

A few tips for writing your own notebook:

Setup

The required packages shoul be installed in advance from the Console.

## SETUP
# INSTALL THESE PACKAGES FROM THE CONSOLE, NOT FROM THIS CHUNK
# paquetes = c("rgdal","raster", "rgeos", terra","elevatr","rasterVis", "rgl", "mapview")
# install.packages(paquetes)

It is advisable to start cleaning the memory:

rm(list=ls())

Now, let’s load the libraries:

library(rasterVis)
library(raster)
#library(terra)
library(rgeos)
library(rgdal)
library(elevatr)
library(sf)
library(tidyverse)
library(mapview)
library(ggplot2)
library(leaflet)
library(tmap)
library(RColorBrewer)

When running a chunk, pay attention to any message indicating error or warning eventual conflicts.

Conflicts may be caused by functions, from differents packages, which have the same names (e.g. raster::extract and tidyr::extract). In order to avoid problems, we need to call such functions using the long way, i.e. package_name::function.

After running a chunk for the first time, it is a good practice to hide messages and warnings using: {r message=FALSE, warning=FALSE}.

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.

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() can be 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 coordinate reference system). The zoom level (z) defaults to 9 (a trade off between resolution and time for download). You could start trying a higher zoom level (e.g. 10).

3. Get raster elevation data for your department

First, we need to load a shapefile or a geopackage representing municipalities of our department. In this notebook, I will use a shapefile downloaded from the DANE’s geoportal.

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

# SpatialPolygonsDataFrame example
munic <-  st_read("./68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `/Users/ivanlizarazo/Documents/GB/68_SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 87 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.52895 ymin: 5.707536 xmax: -72.47706 ymax: 8.14501
## Geodetic CRS:  WGS 84

Let’s check geometry, bounding box, CRS, and attributes of the munic object:

head(munic)

Let’s create a new object with centroids of municipalities:

# The task is to find a center point for each region
# We convert the simple features to spatial polygons
sp.munic = as_Spatial(munic)
# Now we use the *gCentroid* function from the rgeos package
centers <- data.frame(gCentroid(sp.munic, byid = TRUE))
centers$region <- row.names(sp.munic)
centers$munic <- sp.munic$MPIO_CNMBR
centers

Now, let’s download elevation data for our department:

# z denotes the zoom level of the data
# the higher the zoom the higher the spatial resolution
#elevation <- get_elev_raster(munic, z = 12)
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.

What is the elevation object?

elevation
## class      : RasterLayer 
## dimensions : 4085, 3087, 12610395  (nrow, ncol, ncell)
## resolution : 0.0006833192, 0.0006833192  (x, y)
## extent     : -74.53125, -72.42184, 5.615809, 8.407168  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : file69a81303bc7d.tif 
## names      : file69a81303bc7d

Notice a few things about this DEM:

In your notebook, write down the cell size of your raster (i.e. the spatial resolution), in meters.

Note that elevation is a temporary raster layer (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)
# NOTE THAT z10 means Zoom Level 10
# CHANGE IT IF NEEDED
if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("./stder_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}

In my case, elevation data at zoom 12 was about 340.3 MB. When zoom was 10, the file size lowered to about 28.5 MB.

Write down the size of your DEM.

4. Visualize the elevation data

For visualization purposes, a lower resolution speeds the task:

#This chunk is optional
#Use it only when zoom data was very high
elevation2 <- aggregate(elevation, fact=4, fun=mean)
crs(elevation2) <- "+proj=longlat +datum=WGS84"

A good palette is key for an appropiate visualization. Let’s try this one:

pal <- colorNumeric(c("forestgreen","yellow","tan","brown"), values(elevation2),
  na.color = "transparent")

Now, we will use the leaflet library to view the elevation data:

leaflet(munic) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
                          "Km2: ", munic$MPIO_NAREA, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(elevation2, colors = pal, opacity = 0.9)  %>%
  addLegend(pal = pal, values = values(elevation2),
    title = "Elevation data for Santander (mts)")

Clic each region ID to get municipality name and extension. Change the plot options to make more readable your map.

5. Reprojecting 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.

In order to reproject the elevation data, we conduct two steps.

First, we define the target CRS:

#library(sp)
# WGS 84 / UTM zone 18N
spatialref <- CRS(SRS_string="EPSG:32618")

Then, we reproject the DEM using the function projectRaster from the raster package:

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

What we got?

rep_elev
## class      : RasterLayer 
## dimensions : 3869, 2925, 11316825  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 785603.9, 620725.8, 930245.8  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file69a81303bc7d 
## values     : -1413.431, 4502.723  (min, max)

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

(rep_munic = spTransform(sp.munic,spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 87 
## extent      : 552105.4, 778842.2, 630996, 900533.3  (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  :         68,      68001,     AGUADA,                 1539,   19.69444164,      2017,  SANTANDER, 0.270018178002, 0.00160984367921 
## max values  :         68,      68895,   ZAPATOCA, Ordenanza 33 de 1968, 3174.27867176,      2017,  SANTANDER,  4.36038638488,   0.259459176725

To avoid a headache, let’s save our DEM. Make sure to indicate a pathname and a filename which suits your needs (i.e. appropriate for your platform, and your department).

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

6. Basic statistics of elevation data

A quick exploration of the DEM statistics.

This chunk “cleans” elevation lower than 0.0:

rep_elev[rep_elev < 0.0] <- 0.0

Now, let’s create the histogram of elevation values:

# histograma
hist(rep_elev)

We can compute a summary of DEM statistics:

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

Let’s create an unidimensional stats vector:

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

Now, the result:

(df_estadisticas <- data.frame(metricas, valores))

7. Obtention of geomorphometric variables

Before proceeding, make sure you have understood basic concepts about: DEM, slope, and aspect.

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)

Let’s check the first output:

slope
## class      : RasterLayer 
## dimensions : 3869, 2925, 11316825  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 785603.9, 620725.8, 930245.8  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 73.96779  (min, max)

Now, the second output:

aspect
## class      : RasterLayer 
## dimensions : 3869, 2925, 11316825  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 551603.9, 785603.9, 620725.8, 930245.8  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : aspect 
## values     : 0, 360  (min, max)

In order to plot slope, it is a good idea to reduce the size of the RasterLayer:

#This chunk is optional
#Use it only when zoom data was very high
slope2 <- aggregate(slope, fact=4, fun=mean)

It may be convenient to project the elevation dataset:

slope3 <- projectRasterForLeaflet(slope2, "ngb")

Let’s create a color palette:

pal2 <- colorNumeric(c("#537188", "#CBB279", "#E1D4BB", "#EEEEEE"),   values(slope3),na.color = "transparent")

Now, plot time:

leaflet(munic) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
                          "Km2: ", munic$MPIO_NAREA, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(slope3, colors = pal2, opacity = 1.0)  %>%
  addLegend(pal = pal2, values = values(slope3),
    title = "Slope for Santander (%)")

Let’s remind what is the slope2 object:

slope2
## class      : RasterLayer 
## dimensions : 968, 732, 708576  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 551603.9, 785843.9, 620485.8, 930245.8  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0.02368777, 65.15777  (min, max)

Note that there are two libraries to manage raster data in R: - raster, and - terra

The raster library is well know (but, in some tasks, quite slow). The terra library is more recent (but faster).

We will convert a RasterLayer (i.e. a raster object) into a SpatRaster (i.e. a terra object):

(nslope2  <- as(slope2, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 968, 732, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 551603.9, 785843.9, 620485.8, 930245.8  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :       slope 
## min value   :  0.02368777 
## max value   : 65.15776795

Slope classification is a common way to understand relief variability. You need to include here a table with slope intervals, class names, and description, as defined by the Colombian Geographic Institute (IGAC).

Let’s conduct a slope classification based on common ranges. First, we create a reclassification matrix.

m <- c(0, 3, 1,  3, 7, 2, 7,12,3,  12,
       25,  4, 25, 50, 5, 50, 75, 6)
(rclmat <- matrix(m, ncol=3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]    0    3    1
## [2,]    3    7    2
## [3,]    7   12    3
## [4,]   12   25    4
## [5,]   25   50    5
## [6,]   50   75    6

In your notebook, explain the meaning of the above matrix.

Now, the reclassification task:

# for values >= 0 (instead of > 0), do
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)

Let’s check the result:

slope_rec
## class       : SpatRaster 
## dimensions  : 968, 732, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 551603.9, 785843.9, 620485.8, 930245.8  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        : slope 
## min value   :     1 
## max value   :     6

8. Static visualization

We will use tmap for creating and saving a static slope classification map.

slope.map <- tm_shape(slope_rec) +
  tm_raster(alpha = 0.6, style= "cat", 
            title="Slope") +  
  # Borders only
  tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.3, alpha=0.7) +
  tm_scale_bar(text.size = 0.5) +
  tm_compass(position = c("right", "top"), size = 2) +
  tm_graticules(alpha = 0.2) +
  tm_layout(outer.margins = 0.01, legend.position = c("left", "top"), title= 'Slope classes', title.position = c('right', 'top')) +
  tm_credits("Data source: Mapzen & DANE", size=0.3)

Let’s view the output:

slope.map

More tips to plot raster and vector objects can be found at [this link] (https://r-tmap.github.io/tmap-book/). Make sure to try additional options to those included in this notebook.

Now, backup time:

## Note we are using letter paper size
tmap_save(slope.map, "./stder_slope.png", width = 8.5, height =11, dpi=600)
## Map saved to /Users/ivanlizarazo/Documents/GB/stder_slope.png
## Resolution: 5100 by 6600 pixels
## Size: 8.5 by 11 inches (600 dpi)

Make sure to check the output in your working directory

9. Interactive visualization

Now, we want to create an interactive slope map.

Let’s convert the sf object into a SpatVector object:

(nmunic  <- as(rep_munic, "SpatVector"))
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 87, 9  (geometries, attributes)
##  extent      : 552105.4, 778842.2, 630996, 900533.3  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       : DPTO_CCDGO MPIO_CCDGO  MPIO_CNMBR      MPIO_CRSLC MPIO_NAREA
##  type        :      <chr>      <chr>       <chr>           <chr>      <num>
##  values      :         68      68001 BUCARAMANGA            1623      152.9
##                        68      68013      AGUADA            1775      75.23
##                        68      68020     ALBANIA Ordenanza 33 d~      166.2
##  MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
##      <int>      <chr>      <num>      <num>
##       2017  SANTANDER     0.6923    0.01251
##       2017  SANTANDER     0.4758   0.006146
##       2017  SANTANDER     0.8761    0.01357

Let’s convert the raster object into a terra object:

(nslope  <- as(slope, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 3869, 2925, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 551603.9, 785603.9, 620725.8, 930245.8  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 73.96779

Now, let’s calculate slope categories aggregated by municipality:

## S4 method for signature 'SpatRaster,SpatVector'
munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)

Note that the result is a terra object:

munic.slope
## class       : SpatRaster 
## dimensions  : 3869, 2925, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 551603.9, 785603.9, 620725.8, 930245.8  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :     slope 
## min value   :  1.152868 
## max value   : 26.008008

Now, let’s view the map:

tmap_mode("view")
  tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(munic.slope) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Mean slope")+
tm_layout(legend.outside = TRUE)

10. Bibliography

If you reuse this notebook’s code please cite this work as:

Lizarazo, I., 2023. Elevation data processing and analysis in R. Available at https://rpubs.com/ials2un/dem_analysis

Reproducibility

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/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] RColorBrewer_1.1-3 tmap_3.3-3         leaflet_2.1.1      mapview_2.11.0    
##  [5] forcats_0.5.2      stringr_1.5.0      dplyr_1.0.10       purrr_0.3.5       
##  [9] readr_2.1.3        tidyr_1.2.1        tibble_3.1.8       ggplot2_3.4.0     
## [13] tidyverse_1.3.2    sf_1.0-9           elevatr_0.4.2      rgdal_1.6-6       
## [17] rgeos_0.5-9        raster_3.6-20      sp_1.5-1           rasterVis_0.51.5  
## [21] lattice_0.20-45   
## 
## loaded via a namespace (and not attached):
##   [1] googledrive_2.0.0       leafem_0.2.0            colorspace_2.0-3       
##   [4] deldir_1.0-6            ellipsis_0.3.2          class_7.3-20           
##   [7] satellite_1.0.4         markdown_1.4            base64enc_0.1-3        
##  [10] fs_1.5.2                dichromat_2.0-0.1       rstudioapi_0.14        
##  [13] proxy_0.4-27            farver_2.1.1            hexbin_1.28.3          
##  [16] fansi_1.0.3             lubridate_1.9.0         xml2_1.3.3             
##  [19] codetools_0.2-18        cachem_1.0.6            knitr_1.42             
##  [22] jsonlite_1.8.4          tmaptools_3.1-1         broom_1.0.2            
##  [25] dbplyr_2.2.1            png_0.1-8               compiler_4.2.2         
##  [28] httr_1.4.4              backports_1.4.1         assertthat_0.2.1       
##  [31] fastmap_1.1.0           gargle_1.2.1            cli_3.4.1              
##  [34] s2_1.1.1                prettyunits_1.1.1       leaflet.providers_1.9.0
##  [37] htmltools_0.5.4         tools_4.2.2             gtable_0.3.1           
##  [40] glue_1.6.2              wk_0.7.1                Rcpp_1.0.10            
##  [43] cellranger_1.1.0        jquerylib_0.1.4         vctrs_0.5.1            
##  [46] progressr_0.10.1        leafsync_0.1.0          crosstalk_1.2.0        
##  [49] lwgeom_0.2-8            xfun_0.35               rvest_1.0.2            
##  [52] mime_0.12               timechange_0.1.1        lifecycle_1.0.3        
##  [55] XML_3.99-0.13           googlesheets4_1.0.1     terra_1.7-29           
##  [58] zoo_1.8-10              scales_1.2.1            hms_1.1.2              
##  [61] slippymath_0.3.1        parallel_4.2.2          curl_4.3.3             
##  [64] yaml_2.3.7              sass_0.4.4              latticeExtra_0.6-30    
##  [67] stringi_1.7.8           highr_0.10              e1071_1.7-12           
##  [70] rlang_1.0.6             pkgconfig_2.0.3         evaluate_0.19          
##  [73] htmlwidgets_1.6.0       tidyselect_1.2.0        magrittr_2.0.3         
##  [76] R6_2.5.1                generics_0.1.3          DBI_1.1.3              
##  [79] pillar_1.8.1            haven_2.5.1             withr_2.5.0            
##  [82] units_0.8-1             stars_0.5-5             abind_1.4-5            
##  [85] modelr_0.1.10           crayon_1.5.2            interp_1.1-4           
##  [88] KernSmooth_2.23-20      utf8_1.2.2              tzdb_0.3.0             
##  [91] rmarkdown_2.19          progress_1.2.2          jpeg_0.1-10            
##  [94] grid_4.2.2              readxl_1.4.1            reprex_2.0.2           
##  [97] digest_0.6.31           classInt_0.4-8          webshot_0.5.4          
## [100] stats4_4.2.2            munsell_0.5.0           viridisLite_0.4.1      
## [103] bslib_0.4.2