1. Objetivo de este cuaderno.

El objetivo de este cuaderno es ilustrar varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación en R. El trabajo se realizará con un municipio de nuestra preferencia, en este caso, uno perteneciente al departamento de Norte de Santander.

2. Introducción a elevator.

Los datos de elevación tienen múltiples aplicaciones, 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 varios paquetes o requiere acceso local a los datos. Esto ya no es necesario, ya que actualmente, existe una variedad de API (interfaz de programación de aplicaciones) que proporcionan 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 (p.e. un DEM), el paquete elevatr utiliza los mosaicos de terreno de Amazon Web Services. En este cuaderno, se explora esta funcionalidad. Se necesitan instalar los siguientes paquetes:

# install.packages("rgdal")
# install.packages("raster")
# install.packages("elevatr")
# install.packages("rasterVis")
# install.packages("rgl")
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.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/user/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/user/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1
library(elevatr)

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

Existen varias fuentes para modelos de elevación digital como la Misión de Topografía por Radar Shuttle (SRTM), el conjunto de Datos de Elevación Nacional (NED) de USGS, el DEM global (GDEM), y otros. Cada uno de estos DEM posee ventajas y desventajas en 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 dada a un nivel de zoom dado. Los datos compilados por Mapzen son 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 múltiples mosaicos, la salida resultante es una capa ráster fusionada.

Usar get_elev_raster() para acceder a los Terrain Tiles en AWS.

Como se mencionó anteriormente, un data frame con columnas x e y, un objeto sp o un objeto ráster debe ser la entrada y el 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 data frame requiere un prj. Se mostrarán algunos 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 se pudo obtener nada con este zoom. Se tuvo que utilizar un zoom inferior o igual a 8. Primero, se necesita cargar un archivo shapefile que represente a los municipios del departamento de Norte de Santander. En este cuaderno, se va a cargar el shapefile descargado del geoportal del DANE, según lo solicitado en la clase virtual del 19 de marzo de 2020. Se procede a revisar el contenido de la carpeta Norte de Santander:

list.files("C:/Users/user/Documents/Norte de Santander/54_NORTE DE SANTANDER/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, se lee el shapefile usando una función que proporciona el paquete ráster:

(munic = shapefile("C:/Users/user/Documents/Norte de Santander/54_NORTE DE SANTANDER/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 40 
## extent      : -73.63379, -72.04761, 6.872201, 9.290847  (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  :         54,      54001,            ABREGO,                                1555,   44.94540102,      2017, NORTE DE SANTANDER, 0.455016349208, 0.00368621836307 
## max values  :         54,      54874, VILLA DEL ROSARIO, Ordenanza 9 de Noviembre 25 de 1948, 2680.07859017,      2017, NORTE DE SANTANDER,  4.08629755008,   0.220098910874

¿Cuáles son los atributos del objeto munic?

head(munic)
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA MPIO_NANO
## 0         54      54001    CÚCUTA       1972  1133.2031      2017
## 1         54      54003     ABREGO       1806  1382.4498      2017
## 2         54      54051  ARBOLEDAS       1835   456.1490      2017
## 3         54      54099  BOCHALEMA       1826   170.2669      2017
## 4         54      54109 BUCARASICA       1838   270.7909      2017
## 5         54      54125    CÃ\201COTA       1630   138.9222      2017
##           DPTO_CNMBR Shape_Leng Shape_Area
## 0 NORTE DE SANTANDER  2.8636253 0.09289252
## 1 NORTE DE SANTANDER  2.2393759 0.11336240
## 2 NORTE DE SANTANDER  1.1611936 0.03736204
## 3 NORTE DE SANTANDER  0.7565329 0.01394455
## 4 NORTE DE SANTANDER  0.8281364 0.02220518
## 5 NORTE DE SANTANDER  0.5561253 0.01136823

Ahora, se selecciona solo la capital de este departamento.

Cucuta <- munic[munic$MPIO_CNMBR== "CÚCUTA",]
plot(Cucuta, main="Cúcuta", axes = TRUE)
plot(munic, add= TRUE)
invisible(text(coordinates(munic), labels= as.character(munic$MPIO_CNMBR), cex= 0.8))

elevation = get_elev_raster(Cucuta,z=8)
## 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 +ellps=WGS84 +towgs84=0,0,0

Como ya se han descargado y guardado los datos de elevación, se necesita leerlos usando la función ráster del paquete raster. No se debe correr el chunk.

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

Ahora, se verifica qué hay dentro del objeto de elevación:

elevation
## class      : RasterLayer 
## dimensions : 1547, 1033, 1598051  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00272  (x, y)
## extent     : -73.13875, -70.298, 5.601438, 9.809278  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : -267.3149, 5239.53  (min, max)

En este DEM se observan sus dimensiones, el número de filas y columnas de los datos, el número de celdas del ráster,la resolución, la extensión y la crs, que representa 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= "DEM descargado [metros]")
plot(Cucuta, add= TRUE)

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

Tenga en cuenta que el DEM cubre una extensión mayor que la que se necesita. 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.

writeRaster(elevation, filename = "C:/Users/user/Documents/Intro_to_R/win-library3.6elevatr", dataType= "INT4S", overwrite = TRUE)
## Warning in .datatype(...): argument "datatype" misspelled as "dataType"
## class      : RasterLayer 
## dimensions : 1547, 1033, 1598051  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00272  (x, y)
## extent     : -73.13875, -70.298, 5.601438, 9.809278  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/user/Documents/Intro_to_R/win-library3.6elevatr.grd 
## names      : layer 
## values     : -267, 5240  (min, max)
library(elevatr)
elev_crop = crop(elevation, Cucuta)
plot(elev_crop, main = "Modelo de elevación digital recortado")
plot(Cucuta, add= TRUE)

Veamos el nuevo objeto:

elev_crop
## class      : RasterLayer 
## dimensions : 261, 95, 24795  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00272  (x, y)
## extent     : -72.608, -72.34675, 7.720318, 8.430238  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 40.66319, 1757.218  (min, max)

5. Reproyectar los datos de elevación.

Cuando se trabaja con DEMs, es buena idea trabajar con coordenadas de mapa en lugar de coordenadas geográficas. Esto se debe a que, en coordenadas geográficas, las unidades de dimensiones horizontales son grados decimales, pero la unidad de dimensión vertical es metros. Se van a reproyectar los datos de elevación. Se puede visitar la página epsg.io y buscar la proyección MAGNA- SIRGAS Colombia Bogota zone. Se necesita obtener la definición de esta referencia espacial en formato PROJ.4, que es la que se usa para las bibliotecas sp y raster; y se puede copiar este texto PROJ.4 de la misma página epsg.io.

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"

En seguida, se pueden reproyectar los datos de elevación desde las coordenadas geográficas WGS84 en MAGNA Colombia Bogota zone

pr3 <- projectExtent(elev_crop, spatialref)
res(pr3) <- 100
rep_elev <- projectRaster(elev_crop, pr3)
rep_elev
## class      : RasterLayer 
## dimensions : 787, 291, 229017  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1161846, 1190946, 1345724, 1424424  (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     : 40.74319, 1769.531  (min, max)

Ahora, se va a reproyectar el SpatialPolygonsDataFrame que representa la capital del departamento de Norte de Santander:

(rep_Cucuta = spTransform(Cucuta, spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1162057, 1190861, 1345912, 1424479  (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       :         54,      54001,    CÚCUTA,       1972, 1133.20308169,      2017, NORTE DE SANTANDER, 2.8636253172, 0.0928925236454

Después, se plotea:

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

Se puede guardar el DEM:

writeRaster(rep_elev, filename = "C:/Users/user/Documents/Intro_to_R/rep_Cucuta_elev", dataType= "INT4S", overwrite= TRUE)
## Warning in .datatype(...): argument "datatype" misspelled as "dataType"
## class      : RasterLayer 
## dimensions : 787, 291, 229017  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1161846, 1190946, 1345724, 1424424  (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     : C:/Users/user/Documents/Intro_to_R/rep_Cucuta_elev.grd 
## names      : layer 
## values     : 41, 1770  (min, max)

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

A continuación, se realiza una exploración rápida de las estadísticas del DEM. Primero se realiza un histograma.

hist(rep_elev)

promedio <- cellStats(rep_elev, "mean")
minimo <- cellStats(rep_elev, "min")
maximo <- cellStats(rep_elev, "max")
desviacion <- cellStats(rep_elev, "sd")
# Crear un vector unidimensional
metricas <- c("mean", "min", "max", "sd")
valores <- c(promedio, minimo, maximo, desviacion)
# Creación del data frame con estadísticas de elevación [metros]
(df_estadisticas <- data.frame(metricas, valores))
##   metricas    valores
## 1     mean  355.03938
## 2      min   40.74319
## 3      max 1769.53113
## 4       sd  347.61460

6. Obtención de variables geomorfométricas.

Es necesario leer el capítulo geomorfometría de Victor Olaya para seguir ejecutando este cuaderno de R.

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

Ahora, se realiza el ploteo de elevación. Tenga en cuenta la paleta de colores que se utiliza aquí.

plot(rep_elev, main= "DEM para Cúcuta [metros]", col= terrain.colors(25, alpha = 0.7))

Luego, se va a plotear la pendiente. Note la paleta de color usada aquí.

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

Ahora, el aspecto:

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

Plot combinado:

plot (hill, 
      col= grey(1:100/100),
      legend= FALSE,
      main = "DEM para Cúcuta",
      axes= FALSE)
plot(rep_elev, axes= FALSE, 
     col= terrain.colors(12, alpha = 0.35), add= TRUE)

7. Mapeo de datos de elevación con rayshader.

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

# install.packages("rayshader")
library(rayshader)
# Convierta el DEM a una matriz:
elmat = raster_to_matrix(rep_elev)
elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

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

# También, se puede agregar una capa de 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.

Se va a intentar otra forma de visualización que sugiere un experto en R:

#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
}
#map=readJPEG("C:/Users/user/Documents/Intro_to_R/9pvbHjN.jpg")
#out = getv(map, aspect, slope)

#plotRGB(out, main = "Montañas supuestamente bonitas en Cúcuta")