1. Why this notebook?

This is an R Notebook created using R Studio Cloud. It illustrates several functionalities to obtain, process and visualize digital elevation models in R. It aims to help Geomatica Basica students at Universidad Nacional with their home activities during this time of social distancing. Everybody should publish a similar notebook, adapted to cover any municipality of your department (the one you like most), not later than on 2nd April 2020.

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.

# uncomment the **install** lines only once

# install.packages("rgdal")
# install.packages("raster")
# install.packages("elevatr")
#install.packages("rasterVis")
#install.packages("rgl")

# otherwise you will waste time installing the packages again and again
library(rasterVis)
## Loading required package: raster
## Loading required package: sp
## Loading required package: lattice
## Loading required package: latticeExtra
library(raster)
library(rgl)
library(rgdal)
## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.2, released 2017/09/15
##  Path to GDAL shared files: /usr/share/gdal/2.2
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.4-1
library(elevatr)

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 prj. I will show examples using a SpatialPolygonsDataFrame. The zoom level (z) defaults to 9 (a trade off between resolution and time for download), but I could not get anything using this zoom. I had to use a lower zoom equal to 8. Otherwise, the RStudio Cloud notebook would crash again and again. Very often, I had to “reload” the page and start it again.

First, we need to load a shapefile representing municipalities of our department. In this notebook, we will load the shapefile downloaded from the DANE’s geoportal as requested in the virtual class on 19 March 2020. This shapefile was previously uploaded to RStudio Cloud under a folder named antioquia.

Let’s review the content of the folder:

list.files("/cloud/project/antioquia/ADMINISTRATIVO")
##  [1] "MGN_DPTO_POLITICO.cpg"     "MGN_DPTO_POLITICO.dbf"    
##  [3] "MGN_DPTO_POLITICO.prj"     "MGN_DPTO_POLITICO.sbn"    
##  [5] "MGN_DPTO_POLITICO.sbx"     "MGN_DPTO_POLITICO.shp"    
##  [7] "MGN_DPTO_POLITICO.shp.xml" "MGN_DPTO_POLITICO.shx"    
##  [9] "MGN_MPIO_POLITICO.cpg"     "MGN_MPIO_POLITICO.dbf"    
## [11] "MGN_MPIO_POLITICO.prj"     "MGN_MPIO_POLITICO.sbn"    
## [13] "MGN_MPIO_POLITICO.sbx"     "MGN_MPIO_POLITICO.shp"    
## [15] "MGN_MPIO_POLITICO.shp.xml" "MGN_MPIO_POLITICO.shx"

Now, let’s read the shapefile using a function provided by the raster package:

# SpatialPolygonsDataFrame example
(munic <-  shapefile("/cloud/project/antioquia/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 125 
## extent      : -77.12783, -73.88128, 5.418558, 8.873974  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                          MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,     Shape_Leng,      Shape_Area 
## min values  :         05,      05001,  ABEJORRAL,                                1541,   15.83597275,      2017,  ANTIOQUIA, 0.172962481259, 0.0012928296672 
## max values  :         05,      05895,   ZARAGOZA, Ordenanza 9 de Diciembre 15 de 1983, 2959.36315089,      2017,  ANTIOQUIA,  6.61177822771,   0.23894871119

What are the attributes of the munic object:

head(munic)

Let’s select only the capital city of this department.

## make sure you review the documentation to  understand the following lines of code
## https://www.neonscience.org/dc-shapefile-attributes-r
medallo <- munic[munic$MPIO_CNMBR=="MEDELLÍN",]
plot(medallo, main="Medellín", axes=TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

Uncomment the following line to download the elevation data:

#elevation <- get_elev_raster(medallo, z = 8)

As I have already downloaded and saved the elevation data, I need to read it using the raster function provided by the raster package. You should not run the chunk.

## Do not run this
elevation <-  raster("./dem/elev_z8.tif")

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

elevation
## class      : RasterLayer 
## dimensions : 1035, 1033, 1069155  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00273  (x, y)
## extent     : -75.95125, -73.1105, 4.201768, 7.027318  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /cloud/project/dem/elev_z8.tif 
## names      : elev_z8 
## values     : -431, 5255  (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 decinal degrees in this case). Let’s remind that 1 decimal degree represents about 111.11 km.
  • 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.
plot(elevation, main="This the downloaded DEM [meters]")
plot(medallo, add=TRUE)

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 Medellín:

# 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, medallo)
plot(elev_crop, main="Cropped Digital elevation model")
plot(medallo, add=TRUE)

Let’s check the new object:

elev_crop
## class      : RasterLayer 
## dimensions : 77, 90, 6930  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00273  (x, y)
## extent     : -75.72025, -75.47275, 6.164638, 6.374848  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : elev_z8 
## values     : 1232, 3115  (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.

We can go to epsg.io and look for the MAGNA Colombia Bogota zone projection. We need to get the definition for this spatial reference in PROJ.4 format (the one used for the sp and raster libraries. Let’s copy the PROJ.4 text and save it.

spatialref <- "+proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"

Now, we can reproject the elevation data from WGS84 geographic coordinates into MAGNA Colombia Bogota zone.

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

What we got?

rep_elev
## class      : RasterLayer 
## dimensions : 233, 275, 64075  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 818156.7, 845656.7, 1173680, 1196980  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : memory
## names      : elev_z8 
## values     : 1214.608, 3109.653  (min, max)

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

(rep_medallo = spTransform(medallo,spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 818324.2, 845631.6, 1173620, 1196843  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.596200416666666 +lon_0=-74.07750791666666 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +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 
## value       :         05,      05001,   MEDELLÍN,       1965, 374.82795927,      2017,  ANTIOQUIA, 1.0327835241, 0.0306072272822

It is plotting time:

plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_medallo, add=TRUE)

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

# Write the reprojected DEM to disk 
#uncomment and run the following line
#writeRaster(rep_elev, filename="./dem/rep_medallo_elev.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))

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 Medellin [meters]", col=terrain.colors(25,alpha=0.7))

Plot slope. Note the color palette used here.

plot(slope,main="Slope for Medellin [degrees]", col=topo.colors(25,alpha=0.7))

Plot aspect. Note the color palette used here.

plot(aspect,main="Aspect for Medellin [degrees]", col=rainbow(25,alpha=0.7))

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 Medellin",
        axes=FALSE)           # makes for a cleaner plot, if the coordinates aren't necessary

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.35), add=TRUE) # color method

         # sets how transparent the object will be (0=transparent, 1=not transparent)

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)
#We use another one of rayshader's built-in textures:
elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

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

#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) %>%
  plot_map()

8. Another way of visualization

An R expert suggest another way of visualization here.

Let’s try it:

install.packages("jpeg")
## Installing package into '/home/rstudio-user/R/x86_64-pc-linux-gnu-library/3.6'
## (as 'lib' is unspecified)
library(jpeg)
# Spherical environment mapping hillshade function
getv=function(i,a,s){
  ct = dim(i)[1:2]/2
  sx = values(s)/90 * ct[1]
  sy = values(s)/90 * ct[2]
  a = values(a) * 0.01745
  px = floor(ct[1] + sx * -sin(a))
  py = floor(ct[2] + sy * cos(a))
  
  
  template = brick(s,s,s)
  values(template)=NA
  
  cellr = px + py * ct[1]*2
  cellg = px + py * ct[1]*2 + (ct[1]*2*ct[2]*2)
  cellb = px + py * ct[1]*2 + 2*(ct[1]*2*ct[2]*2)
  
  template[[1]] = i[cellr]
  template[[2]] = i[cellg]
  template[[3]] = i[cellb]
  
  template = template * 256
  
  template
}
# Load an environment map image
# An example can be found at http://i.imgur.com/9pvbHjN.jpg
map=readJPEG("./dem/9pvbHjN.jpg")

# Do the mapping
out = getv(map, aspect, slope)

# Plot pretty mountains
plotRGB(out, main = "Supposedly pretty mountains in Medellin")

Now, it is time for you to start coding.

In case of running into trouble, do not blame neither R nor RStudio. Please review the documentation. Look for examples and you will learn a lot.

Good luck!