Elevation data in R

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.

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

list.files("/cloud/project/valledelcauca/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:

(munic <-  shapefile("/cloud/project/valledelcauca/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 42 
extent      : -77.54977, -75.70724, 3.091239, 5.047394  (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  :         76,      76001,     ALCALÁ,                          1536,   41.86090736,      2017, VALLE DEL CAUCA, 0.453826056161, 0.00340901158098 
max values  :         76,      76895,     ZARZAL, Ordenanza 9 de Diciembre 1954, 6292.50083741,      2017, VALLE DEL CAUCA,  6.59527269127,   0.510778519244 

What are the attributes of the munic object:

head(munic)

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

cali <- munic[munic$MPIO_CNMBR=="CALI",]
plot(cali, main="", axes=TRUE)
plot(munic, add=TRUE)

invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

elevation <-get_elev_raster(cali, z = 8)

Downloading DEMs [===================>-------]  75% eta:  0s
Downloading DEMs [===========================] 100% eta:  0s
Merging DEMs
Reprojecting DEM to original projection
no non-missing arguments to min; returning Infno non-missing arguments to max; returning -InfNote: Elevation units are in meters.
Note: The coordinate reference system is:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

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

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

elevation

class      : RasterLayer 
dimensions : 1035, 1033, 1069155  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -77.3575, -74.51675, 1.392743, 4.228643  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -598.1536, 5331.562  (min, max)

plot(elevation, main="This the downloaded DEM [meters]")
plot(cali, 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.

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

elev_crop = crop(elevation, cali)
plot(elev_crop, main="Cropped Digital elevation model")
plot(cali, add=TRUE)

Let’s check the new object:

elev_crop

class      : RasterLayer 
dimensions : 101, 90, 9090  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -76.70575, -76.45825, 3.272383, 3.549123  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 944.2498, 3979.236  (min, max)

5. Reproject the elevation data

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)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

rep_elev

class      : RasterLayer 
dimensions : 307, 276, 84732  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 707796.8, 735396.8, 853929.3, 884629.3  (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      : layer 
values     : 943.4611, 3963.625  (min, max)

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

(rep_cali = spTransform(cali,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 707694.8, 735285.7, 854051.1, 884460.6  (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       :         76,      76001,       CALI,       1536, 563.04276445,      2017, VALLE DEL CAUCA, 1.16366974876, 0.0457332109324 

It is plotting time:

plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_cali, add=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 CALI [meters]", col=terrain.colors(25,alpha=0.7))

Plot slope. Note the color palette used here.

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

Plot aspect. Note the color palette used here.

plot(aspect,main="Aspect for CALI [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 CALI",
        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

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.

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

[1] "Dimensions of matrix are: 276x307."

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

Calculating Surface Normal [-----------------------------------] ETA:  0s
Calculating Surface Normal [===================================] ETA:  0s
                                                                         

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

Calculating Surface Normal [-----------------------------------] ETA:  0s
Calculating Surface Normal [===================================] ETA:  0s
                                                                         

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

Calculating Surface Normal [-----------------------------------] ETA:  0s
Calculating Surface Normal [===================================] ETA:  0s
                                                                         

Raytracing [---------------------------------------------------] ETA:  0s
Raytracing [===================================================] ETA:  0s
                                                                         

8. Another way of visualization

An R expert suggest another way of visualization here.

Let’s try it:

# 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("./valledelcauca/9pvbHjN.jpg")

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

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

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