Elevation data in R

By Valentina Valle Velasco

April 2nd, 2020

This is an R Notebook with the aim of illustrates a lot of functionalities to obtain, process and visualize Digital Elevation Models (DEM) in R.

*Introduction to elevatr.

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

library(rasterVis)

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.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/vale_/Documents/R/win-library/3.6/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/vale_/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.4-1 

library(elevatr)

*Get Raster Elevation Data

Let’s review the content of the folder:

list.files("c:/Users/vale_/Desktop/UNAL/6to semestre/GB/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"    
 [2] "MGN_DPTO_POLITICO.dbf"    
 [3] "MGN_DPTO_POLITICO.prj"    
 [4] "MGN_DPTO_POLITICO.sbn"    
 [5] "MGN_DPTO_POLITICO.sbx"    
 [6] "MGN_DPTO_POLITICO.shp"    
 [7] "MGN_DPTO_POLITICO.shp.xml"
 [8] "MGN_DPTO_POLITICO.shx"    
 [9] "MGN_MPIO_POLITICO.cpg"    
[10] "MGN_MPIO_POLITICO.dbf"    
[11] "MGN_MPIO_POLITICO.prj"    
[12] "MGN_MPIO_POLITICO.sbn"    
[13] "MGN_MPIO_POLITICO.sbx"    
[14] "MGN_MPIO_POLITICO.shp"    
[15] "MGN_MPIO_POLITICO.shp.xml"
[16] "MGN_MPIO_POLITICO.shx"    

Read the shapefile using a function provided by the raster package:

(munic <- shapefile("c:/Users/vale_/Desktop/UNAL/6to semestre/GB/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 30 
extent      : -74.9466, -73.54184, 8.936489, 11.34891  (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  :         47,      47001,     ALGARROBO,                 1525,  109.48370634,      2017,  MAGDALENA, 0.544326259962, 0.00903317812539 
max values  :         47,      47980, ZONA BANANERA, Ordenanza 74 de 1912, 2347.13929515,      2017,  MAGDALENA,  3.19741434448,   0.194233330475 

What are the attributes of the Munic object:

head(munic)

Now, we are going to select only the capital city of this department.

la_perla <- munic[munic$MPIO_CNMBR=="SANTA MARTA",]    
plot(la_perla, main="Santa Marta", axes=TRUE)             
plot(munic, add=TRUE)  

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

elevation <- get_elev_raster(la_perla, z = 8)

Downloading DEMs [========>------------------]  33% eta:  2s
Downloading DEMs [=============>-------------]  50% eta:  2s
Downloading DEMs [=================>---------]  67% eta:  1s
Downloading DEMs [=====================>-----]  83% eta:  1s
Downloading DEMs [===========================] 100% eta:  0s
Merging DEMs
Reprojecting DEM to original projection
ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -InfNote: Elevation units are in meters.
Note: The coordinate reference system is:
 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

What will be inside?

elevation

class      : RasterLayer 
dimensions : 1546, 1033, 1597018  (nrow, ncol, ncell)
resolution : 0.00275, 0.0027  (x, y)
extent     : -74.545, -71.70425, 8.393864, 12.56806  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -3936.986, 5539.776  (min, max)

plot(elevation, main="This the downloaded DEM [meters]")                                    
plot(la_perla, add=TRUE)

*Crop the elevation data to match the study area extent

writeRaster(elevation, filename="C:/Users/vale_/Documents/R/win-library/3.6/elevatr", datatype ='INT4S', overwrite=TRUE)

class      : RasterLayer 
dimensions : 1546, 1033, 1597018  (nrow, ncol, ncell)
resolution : 0.00275, 0.0027  (x, y)
extent     : -74.545, -71.70425, 8.393864, 12.56806  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/vale_/Documents/R/win-library/3.6/elevatr.grd 
names      : layer 
values     : -3937, 5540  (min, max)

elev_crop = crop(elevation, la_perla)        
plot(elev_crop, main="Cropped Digital Elevation Model")                           
plot(la_perla, add=TRUE)

We are going to checke the new object:

elev_crop

class      : RasterLayer 
dimensions : 194, 246, 47724  (nrow, ncol, ncell)
resolution : 0.00275, 0.0027  (x, y)
extent     : -74.2425, -73.566, 10.82386, 11.34766  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -237.3386, 5539.776  (min, max)

*Reproject the elevation data
Save the PROJ.4

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"

pr3 <- projectExtent(elev_crop, spatialref)   

res(pr3) <- 100

rep_elev <- projectRaster(elev_crop, pr3)

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

rep_elev

class      : RasterLayer 
dimensions : 194, 246, 47724  (nrow, ncol, ncell)
resolution : 300.7187, 298.9142  (x, y)
extent     : 981957.8, 1055935, 1688750, 1746740  (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     : -236.208, 5535.2  (min, max)

Now, let’s reproject the SpatialPolygonsDateFrame representing the capital of our departmen:

rep_la_perla = spTransform(la_perla,spatialref)

Now, plot:

plot(rep_elev, main="Reprojected Digital Elevation Model")                            
plot(rep_la_perla, add=TRUE)

To avoid problems, save our DEM:

writeRaster(rep_elev, filename = "C:/Users/vale_/Desktop/UNAL/6to semestre/GB/rep_la_perla_elev.tif", datatype='INT4S', overwrite=TRUE)

*Basic statistics of elevation data

Now, we are going to explore of the DEM statistics:

hist(rep_elev)

promedio <- cellStats(rep_elev, 'mean')      
minimo <- cellStats(rep_elev, 'min')         
maximo <- cellStats(rep_elev, 'max')         
desviacion <- cellStats(rep_elev, 'sd')

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

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

*Obtention of geomorphometric variables

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.

plot(rep_elev,main="DEM for Santa Marta [meters]", col=terrain.colors(25,alpha=0.7))

Plot slope, with other color palette

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

Plot aspect. It’s with other color palette.

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

A combined plot:

plot(hill, col=grey(1:100/100), legend=FALSE, main="DEM for Santa Marta", axes=FALSE)      
plot(rep_elev, axes=FALSE, col=terrain.colors(12, alpha=0.35), add=TRUE)

*Mapping elevation data with rayshader

Rayshader is useful to create amazing 2D and 3D maps.

# install.packages("rayshader")

library(rayshader)

Convert the DEM into a matrix:

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 246x194."

elmat %>% sphere_shade(texture = "imhof2") %>% plot_map()

Add a water layer to the map, with the next code:

elmat %>% sphere_shade(texture = "desert") %>% add_water(detect_water(elmat), color = "desert") %>% plot_map()

elmat %>% sphere_shade(texture = "desert") %>% add_water(detect_water(elmat), color = "desert") %>% add_shadow(ray_shade(elmat), 0.5) %>% plot_map()

*Another way of visualization

#install.packages("jpeg")

library(jpeg)

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, do the mapping and plot pretty mountains

map=readJPEG("C:/Users/vale_/Desktop/UNAL/6to semestre/GB/9pvbHjN.jpg")
out = getv(map, aspect, slope)              

plotRGB(out, main = "Supposedly pretty mountains in Santa Marta")