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 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 22 September 2020.

Note that this notebook complements the information provided in an early notebook described last week in class.

A few tips for writing your own notebook:

It may be that you find helpful the following information:

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.

# 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)
library(raster)
library(rgl)
library(rgdal)
library(elevatr)

It is important to know what is the working directory:

getwd()
## [1] "/Users/ivanlizarazo/Documents/ivan/GB/notebooks"

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 SpatialPolygonsDataFrame. 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 = 12. Previously, using RStudio Cloud, 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 RStudio Cloud page and start it again. Using RStudio on my own laptop, I did not have many problems.

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 review the content of the administrativo folder:

list.files("../datos/shapefiles/15_boyaca/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("../datos/shapefiles/15_boyaca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : -74.66496, -71.94885, 4.655196, 7.055557  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 2
## names       : MPIO_CCDGO, MPIO_CNMBR 
## min values  :      15001,    ALMEIDA 
## max values  :      15897,  ZETAQUIRÁ

What are the attributes of the munic object:

head(munic)

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

tunja <- munic[munic$MPIO_CNMBR=="TUNJA",]

Plotting time:

## make sure you review the documentation to  understand the following lines of code
## https://www.neonscience.org/dc-shapefile-attributes-r
plot(tunja,  axes=TRUE, col="green")
plot(munic, add=TRUE)
text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR))

Uncomment the following line to runk the code and download the elevation data. Once the code has run, comment the chunk.

#elevation <- get_elev_raster(tunja, z = 12)

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/tunja_dem_z12.tif"), format="GTiff", #overwrite=TRUE)
#}

In my case, running the code above produced a .TIF file:

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/tunja_dem_z12.tif")

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

elevation
## class      : RasterLayer 
## dimensions : 1545, 1543, 2383935  (nrow, ncol, ncell)
## resolution : 0.000172, 0.000171  (x, y)
## extent     : -73.47742, -73.21203, 5.352646, 5.616841  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : /Users/ivanlizarazo/Documents/ivan/GB/datos/tunja_dem_z12.tif 
## names      : tunja_dem_z12 
## values     : 2026.88, 3543.023  (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 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.
plot(elevation, main="This the downloaded DEM for Tunja [meters]")
plot(tunja, add=TRUE)
text(coordinates(tunja), labels=as.character(tunja$MPIO_CNMBR))

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, tunja)
plot(elev_crop, main="Cropped Digital elevation model for Tunja")
plot(tunja, add=TRUE)
text(coordinates(tunja), labels=as.character(tunja$MPIO_CNMBR))

Let’s check the new object:

elev_crop
## class      : RasterLayer 
## dimensions : 805, 834, 671370  (nrow, ncol, ncell)
## resolution : 0.000172, 0.000171  (x, y)
## extent     : -73.44732, -73.30387, 5.446354, 5.584009  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : tunja_dem_z12 
## values     : 2231.158, 3262.087  (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) <- 20
# now project
rep_elev <- projectRaster(elev_crop, pr3)

What we got?

rep_elev
## class      : RasterLayer 
## dimensions : 762, 796, 606552  (nrow, ncol, ncell)
## resolution : 20, 20  (x, y)
## extent     : 1069823, 1085743, 1094051, 1109291  (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      : tunja_dem_z12 
## values     : 2228.549, 3261.334  (min, max)

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

(rep_tunja = spTransform(tunja,spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1069842, 1085712, 1094054, 1109295  (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   : 2
## names       : MPIO_CCDGO, MPIO_CNMBR 
## value       :      15001,      TUNJA

It is plotting time:

plot(rep_elev, main="Reprojected Digital elevation model for Tunja")
plot(rep_tunja, add=TRUE)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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_tunja_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))

Do not forget to describe the DEM statistics.

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 Tunja [meters]", col=terrain.colors(25,alpha=0.8))
plot(rep_tunja, add=TRUE, lwd=2)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

Plot slope. Note the color palette used here.

plot(slope,main="Slope for Tunja [degrees]", col=topo.colors(25,alpha=0.7))
plot(rep_tunja, add=TRUE, lwd=2)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

Plot aspect. Note the color palette used here.

plot(aspect,main="Aspect for Tunja [degrees]", col=rainbow(25,alpha=0.7))
plot(rep_tunja, add=TRUE, lwd=2)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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 Tunja",
        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)
plot(rep_tunja, add=TRUE, lwd=2)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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")
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
}

Before going on, make sure to download to your working directory this jpg file:

# Load an environment map image
# An example can be found at http://i.imgur.com/9pvbHjN.jpg
map=readJPEG("../datos/9pvbHjN.jpg")

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

# Plot pretty mountains
plotRGB(out, main = "Supposedly pretty mountains in Tunja")
polygons(rep_tunja)
## class       : SpatialPolygons 
## features    : 1 
## extent      : 1069842, 1085712, 1094054, 1109295  (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
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

Now, it is time for you to start coding. Note that this notebook is just an example. You can put into practice your own DEM visualization 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, do not blame neither R nor RStudio. Please review the documentation. Look for examples and you will learn a lot.

Good luck!

8. Reproducibility

sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] jpeg_0.1-8.1        rayshader_0.19.2    elevatr_0.2.0      
##  [4] rgdal_1.4-8         rgl_0.100.54        rasterVis_0.47     
##  [7] latticeExtra_0.6-29 lattice_0.20-41     raster_3.3-13      
## [10] sp_1.4-2           
## 
## loaded via a namespace (and not attached):
##  [1] progress_1.2.2          zoo_1.8-8               xfun_0.17              
##  [4] colorspace_1.4-1        vctrs_0.3.4             miniUI_0.1.1.1         
##  [7] htmltools_0.5.0         viridisLite_0.3.0       yaml_2.2.1             
## [10] rayimage_0.3.1          rlang_0.4.7             manipulateWidget_0.10.1
## [13] hexbin_1.28.1           later_1.1.0.1           RColorBrewer_1.1-2     
## [16] lifecycle_0.2.0         foreach_1.5.0           stringr_1.4.0          
## [19] munsell_0.5.0           htmlwidgets_1.5.1       codetools_0.2-16       
## [22] evaluate_0.14           knitr_1.29              fastmap_1.0.1          
## [25] doParallel_1.0.15       httpuv_1.5.4            crosstalk_1.1.0.1      
## [28] parallel_3.6.3          Rcpp_1.0.5              xtable_1.8-4           
## [31] scales_1.1.1            promises_1.1.1          webshot_0.5.2          
## [34] jsonlite_1.7.1          mime_0.9                hms_0.5.3              
## [37] png_0.1-7               digest_0.6.25           stringi_1.5.3          
## [40] shiny_1.5.0             grid_3.6.3              tools_3.6.3            
## [43] magrittr_1.5            crayon_1.3.4            pkgconfig_2.0.3        
## [46] prettyunits_1.1.1       rmarkdown_2.3.5         iterators_1.0.12       
## [49] R6_2.4.1                compiler_3.6.3