Datos de elevacion en R para el departamento de Huila

1. ¿Por qué este cuaderno?

Este es un R Notebook creado con RStudio. Ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación en R. Tiene como objetivo ayudar a los estudiantes de Geomatica Básica de la Universidad Nacional con sus actividades en el hogar durante este tiempo de distanciamiento social. Todo el mundo debería publicar un cuaderno similar, adaptado para cubrir cualquier municipio de su departamento (el que más le guste).

Comenzamos introduciendo estos comandos para mejorar el rendimiento del código y mejorar su publicación:

knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

2. Introducción a elevatr

Los datos de elevación se utilizan para una amplia gama de aplicaciones, que incluyen, por ejemplo, visualización, hidrología y modelado ecológico. Obtener acceso a estos datos en R no ha tenido una sola interfaz, está disponible a través de funciones en muchos paquetes o requiere acceso local a los datos. Esto ya no es necesario, ya que existe una variedad de API que brindan acceso programático a los datos de elevación. El paquete elevatr se escribió para estandarizar el acceso a los datos de elevación desde las API web.

Para acceder a los datos de elevación ráster (por ejemplo, un DEM), el paquete elevatr utiliza los mosaicos de terreno de Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.

# descomenta las líneas ** install ** solo una vez

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

# de lo contrario perderá tiempo instalando los paquetes una y otra vez

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.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/Sergio/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/Sergio/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

library(elevatr)

3. Obtener datos de elevación de ráster

Hay varias fuentes de modelos digitales de elevación, como Shuttle Radar Topography Mission (SRTM), USGS National Elevation Dataset (NED), Global DEM (GDEM) y otras. Cada uno de estos DEM tiene ventajas y desventajas para su uso. Antes de su cierre en enero de 2018, Mapzen combinó varias de estas fuentes para crear un producto de elevación de síntesis que utiliza los mejores datos de elevación disponibles para una región determinada a un nivel de zoom determinado. Aunque cerrados, estos datos compilados por Mapzen son actualmente accesibles a través de Terrain Tiles en Amazon Web Services (AWS).

La entrada para get_elev_raster () es un data.frame con ubicaciones x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto sp o cualquier objeto ráster y devuelve un RasterLayer de los mosaicos que se superponen al cuadro delimitador de la entrada. Si se recuperan varios mosaicos, la salida resultante es una capa ráster combinada.

Usando get_elev_raster () para acceder a Terrain Tiles en AWS. Como se mencionó, un marco de datos con columnas xey, un objeto sp o un objeto raster debe ser la entrada y src debe establecerse en “mapzen” (este es el valor predeterminado).

No hay diferencia en el uso de los tipos de datos de entrada sp y ráster. El marco de datos requiere un prj. Mostraré ejemplos usando un SpatialPolygonsDataFrame. El nivel de zoom (z) predeterminado es 9 (una compensación entre la resolución y el tiempo de descarga), pero no pude obtener nada con este zoom. Tuve que usar un zoom más bajo igual a 8. De lo contrario, el portátil RStudio Cloud se bloquearía una y otra vez. Muy a menudo, tuve que “recargar” la página y comenzar de nuevo.

Primero, necesitamos cargar un shapefile que represente a los municipios de nuestro departamento. En este cuaderno, cargaremos el shapefile descargado del geoportal del DANE según lo solicitado en la clase virtual del 19 de marzo de 2020. Este shapefile fue previamente subido a RStudio Cloud en una carpeta denominada antioquia.

Revisemos el contenido de la carpeta:

list.files("C:/Users/Sergio/Documents/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"

Ahora, leamos el shapefile usando una función proporcionada por el paquete raster

# Ejemplo de SpatialPolygonsDataFrame
(munpios <-  shapefile("C:/Users/Sergio/Documents/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

## class       : SpatialPolygonsDataFrame 
## features    : 37 
## extent      : -76.62466, -74.41303, 1.552125, 3.843208  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +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  :         41,      41001,    ÍQUIRA,                                    1538,   80.44174325,      2017,      HUILA, 0.517464409731, 0.00653203068962 
## max values  :         41,      41885,   YAGUARÁ, Ordenanza 6 del 26 de Noviembre de 1965, 1584.96991931,      2017,      HUILA,  3.07318677018,   0.128978175097

Cuáles son los atributos del objeto munpios:

head(munpios)

Seleccionemos solo la ciudad capital de este departamento.

## asegúrese de revisar la documentación para comprender las siguientes líneas de código
## https://www.neonscience.org/dc-shapefile-attributes-r
neiva <- munpios[munpios$MPIO_CNMBR=="NEIVA",]
plot(neiva, main="Neiva", axes=TRUE)
plot(munpios, add=TRUE)
invisible(text(coordinates(munpios), labels=as.character(munpios$MPIO_CNMBR), cex=0.8))

Descomente la siguiente línea para descargar los datos de elevación:

elevation <- get_elev_raster(neiva, z = 8)

## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output

## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output

## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition

## Merging DEMs

## Reprojecting DEM to original projection

## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  +proj=longlat +datum=WGS84 +no_defs

Como ya descargué y guardé los datos de elevación, necesito leerlos usando la función ráster proporcionada por el paquete ráster. No debe ejecutar el fragmento.

## No ejecute esto
#elevation <- raster("./dem/elev_z8.tif")

Ahora, vamos a revisar qué hay dentro del objeto de elevación:

elevation

## class      : RasterLayer 
## dimensions : 1548, 1033, 1599084  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -75.95125, -73.1105, -0.01287686, 4.228643  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -184.4623, 4500.608  (min, max)

Observe algunas cosas sobre este DEM:

• dimensiones: el “tamaño” del archivo en píxeles
• nrow, ncol: el número de filas y columnas en los datos (imagine una hoja de cálculo o una matriz).
• ncells: el número total de píxeles o celdas que componen el ráster.
• resolución: el tamaño de cada píxel (en grados decimales en este caso). Recordemos que 1 grado decimal representa aproximadamente 111,11 km.
• extensión: la extensión espacial del ráster. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del ráster.
• crs: la cadena del sistema de referencia de coordenadas para el ráster. Este ráster está en coordenadas geográficas con un datum de WGS 84.

plot(elevation, main= "Este es el DEM descargado [metros]")
plot(neiva, add = TRUE)

4. Recorta los datos de elevación para que coincidan con la extensión del área de estudio.

Note que el DEM cubre una extensión mayor que la que necesitamos. Esto se debe a la disposición espacial de los mosaicos de elevación en AWS.

De todos modos, es una buena idea guardar el DEM para el futuro.

# Escriba el RasterLayer en el disco (consulte la documentación de tipos de datos para otros formatos)
# Revise el significado del tipo de datos en el siguiente enlace
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/dataType
# El nombre del archivo incluye un prefijo relacionado con el nivel de zoom (igual a 8 en este caso)
# Descomente y ejecute la siguiente línea:
writeRaster(elevation, filename="C:/Users/Sergio/Documents", datatype='INT4S', overwrite=TRUE)

## class      : RasterLayer 
## dimensions : 1548, 1033, 1599084  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -75.95125, -73.1105, -0.01287686, 4.228643  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : C:/Users/Sergio/Documents.grd 
## names      : layer 
## values     : -184, 4501  (min, max)

Ahora, recortemos los datos de elevación correspondientes a Neiva:

# revise la documentación ráster para comprender cómo funciona la función de recorte
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/crop
elev_crop = crop(elevation, neiva)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(neiva, add=TRUE)

Revisemos el nuevo objeto:

elev_crop

## class      : RasterLayer 
## dimensions : 204, 246, 50184  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -75.62125, -74.94475, 2.724383, 3.283343  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 347.5017, 3157.259  (min, max)

5. Vuelva a proyectar los datos de elevación

Al trabajar con DEM, siempre es una buena idea utilizar coordenadas de mapa en lugar de coordenadas geográficas. Esto se debe al hecho de que, en coordenadas geográficas, las unidades de dimensiones horizontales son grados decimales PERO la unidad de dimensión vertical son metros. Reproyectemos los datos de elevación.

Podemos ir a epsg.io y buscar la proyección de la zona MAGNA Colombia Bogotá. Necesitamos obtener la definición de esta referencia espacial en formato PROJ.4 (el que se usa para las bibliotecas sp y ráster. Copiemos el texto PROJ.4 y guárdelo.

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" 

Ahora, podemos reproyectar los datos de elevación de las coordenadas geográficas WGS84 en la zona MAGNA Colombia Bogotá.

# Proyecto Ráster
# para tener el control, creamos un objeto raster
# usando projectExtent (no se transfieren valores)
pr3 <- projectExtent(elev_crop, spatialref)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

# Ajustar el tamaño de la celda
res(pr3) <- 100
# Ahora Proyecto
rep_elev <- projectRaster(elev_crop, pr3)

¿Qué tenemos?

rep_elev

## class      : RasterLayer 
## dimensions : 619, 753, 466107  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 828323.5, 903623.5, 793057.2, 854957.2  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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     : 349.1514, 3141.076  (min, max)

Ahora, volvamos a proyectar el SpatialPolygonsDataFrame que representa la capital de nuestro departamento:

(rep_neiva = spTransform(neiva,spatialref))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 828277.2, 903590.1, 792961.5, 854990  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +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       :         41,      41001,      NEIVA,       1612, 1269.82755999,      2017,      HUILA, 3.07318677018, 0.103253352109

Es tiempo de trazar:

plot(rep_elev, main="Modelo de elevación digital reproyectado")
plot(rep_neiva, add=TRUE)

Para evitar un dolor de cabeza, guardemos nuestro DEM:

# Escriba el DEM reproyectado en el disco
# descomenta y ejecuta la siguiente línea
writeRaster (rep_elev, filename = "C:/Users/Sergio/Documents", datatype = 'INT4S', overwrite = TRUE)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## class      : RasterLayer 
## dimensions : 619, 753, 466107  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 828323.5, 903623.5, 793057.2, 854957.2  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : C:/Users/Sergio/Documents.grd 
## names      : layer 
## values     : 349, 3141  (min, max)

5. Estadísticas básicas de datos de elevación

Una exploración rápida de las estadísticas de DEM:

# histograma
hist(rep_elev)

# funciona para archivos grandes
promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion  <- cellStats(rep_elev, 'sd')

# crear vector unidimensional
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. Obtención de variables geomorfométricas

Primero, calcule la pendiente, el aspecto y el sombreado:

slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)

Elevación de la parcela. Tenga en cuenta la paleta de colores utilizada aquí.

plot(rep_elev,main="DEM para Huila [metros]", col=terrain.colors(25,alpha=0.8))

Pendiente de la parcela. Tenga en cuenta la paleta de colores utilizada aquí.

plot(slope,main="Pendiente para Huila [grados]", col=topo.colors(25,alpha=0.7))

Aspecto de la parcela. Tenga en cuenta la paleta de colores utilizada aquí.

plot(aspect,main="Aspecto para Huila [grados]", col=rainbow(25,alpha=0.7))

En caso de que esté interesado, lea sobre las paletas de colores en R.

Una trama combinada:

# trazar DEM con sombreado
plot(hill,
        col=grey(1:100/100), # crea una rampa de color de colores grises para sombreado  
        legend=FALSE,  # sin leyenda, no nos importa el gris de la sombra       
        main="DEM para Huila",
        axes=FALSE) # hace una trama más limpia, si las coordenadas no son necesarias

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.35), add=TRUE) # método de color

# establece qué tan transparente será el objeto (0 = transparente, 1 = no transparente)

7. Mapeo de datos de elevación con rayshader

La biblioteca rayshader es un paquete de código abierto para producir visualizaciones de datos 2D y 3D en R. rayshader usa datos de elevación en una matriz R base y una combinación de trazado de rayos, mapeo de textura esférica, superposiciones y oclusión ambiental para generar hermosos mapas topográficos 2D y 3D . Además de los mapas, rayshader también permite al usuario traducir objetos ggplot2 en hermosas visualizaciones de datos en 3D.

# install.packages ("rayshader")

library(rayshader)

#Convierta el DEM en una matriz:
elmat = raster_to_matrix(rep_elev)

# Usamos otra de las texturas integradas de rayshader:
elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

#detect_water y add_water agregan una capa de agua al mapa:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_map()

#Y también podemos agregar una capa con trazado de rayos desde esa dirección del sol:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  plot_map()

8. Otra forma de visualización

Un experto en R sugiere aquí otra forma de visualización.

Vamos a intentarlo:

# install.packages("jpeg")

library(jpeg)

# Función de sombreado de mapeo de entorno esférico
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
}

# Cargar una imagen de mapa del entorno
# Se puede encontrar un ejemplo en http://i.imgur.com/9pvbHjN.jpg
map=readJPEG("C:/Users/Sergio/Documents/9pvbHjN.jpg")

# Hacer el mapeo
out = getv(map, aspect, slope)

# Traza montañas
plotRGB(out, main = "Supposedly pretty mountains in Huila")