1. Introducción

Por medio de este cuaderno se mostrará como obtener, procesar y visualizar modelos digitales de elevación (DEM) en R a partir de los datos de Cundinamarca

Los modelos digitales de elevación son representaciones visuales y matemáticas de los valores de altura con respecto al nivel medio del mar. Estos permiten caracterizan las formas del relive.

El paquete de datos “elevatr” se desarrolló para estandarizar el acceso a los datos de elevación desde la API web. Este utiliza mosiacos de terreno de Amazon Web Services

Primero se deben isntalar y cargar las siguientes librerías para que el paquete “elevatr” pueda funcionar correctamente.

#install.packages("rgdal")
#install.packages("elevatr")
#install.packages("lattice")
#install.packages("latticeExtra")
#install.packages("rasterVis")
#install.packages("rgl")
library(rgdal)
## Loading required package: sp
## 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/JUAN PABLO/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/JUAN PABLO/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)
library(lattice)
library(latticeExtra)
library(rasterVis)
## Loading required package: raster
library(rgl)
library(raster)
library(sp)

2. Obtener datos de elevación tipo ráster

Existen diferentes fuentes para generar un DEM como huttle Radar Topography Mission (SRTM), USGS National Elevation Dataset (NED), Global DEM (GDEM), entre otros.

Actualmente Mapzen combinó varias de las fuentes mencionadas para crear un producto de elevación que utiliza los mejores datos disponibles. Se puede obtener acceso a dichos datos a través de Terrain Tiles en Amazon Web Services (AWS)

Usar la función "get_elev_raster() para acceder a los datos de Tarrain Tiles en AWS

Un data frame corresponde a un conjunto de datos con colummnas e hileras. Los objetos tipo sp o raster deben ser los objetos de entrada y el src debe establecerse en “mapzen” (este es el valor predetermiado)

No hay diferencia en el uso de datos tipo sp o ráster. El marco de datos requiere una extensión de archivo tipo .prj. Se mostrarán unos ejemplos usando “SpatailPolygonsDaraFrame”.

El nivel de zoom (z), predeterminado es 9, pero dicha referencia presenta problemas. Por lo tanto se usará el zoom más bajo igual a 8.

Inicialmente, se debe cargar un shapfile que represente la división adminsitrativa de Cundinamarca. Se cargará el archivo tipo shape del geoportal del DANE, correspondiente a los datos del departamento del año 2017

(munic<- shapefile("C:/Users/JUAN PABLO/Documents/PROGRAMAS GEOMATICA/25_CUNDINAMARCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 116 
## extent      : -74.89063, -73.05256, 3.730129, 5.837258  (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  :         25,      25001, AGUA DE DIOS,                              1537,   43.05740684,      2017, CUNDINAMARCA, 0.292826133295, 0.00351103615261 
## max values  :         25,      25899,    ZIPAQUIRA, Ordenanza 79 Noviembre 29 de 1963, 1196.74949014,      2017, CUNDINAMARCA,  1.99993001495,  0.0975020051483
munic
## class       : SpatialPolygonsDataFrame 
## features    : 116 
## extent      : -74.89063, -73.05256, 3.730129, 5.837258  (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  :         25,      25001, AGUA DE DIOS,                              1537,   43.05740684,      2017, CUNDINAMARCA, 0.292826133295, 0.00351103615261 
## max values  :         25,      25899,    ZIPAQUIRA, Ordenanza 79 Noviembre 29 de 1963, 1196.74949014,      2017, CUNDINAMARCA,  1.99993001495,  0.0975020051483

Para ver los atributos de este archivo se usa la función “head”

head(munic)
##   DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR                       MPIO_CRSLC MPIO_NAREA
## 0         25      25483    NARIÑO                             1899   55.16263
## 1         25      25513      PACHO                             1604  402.62678
## 2         25      25506    VENECIA Decreto 727 Septiembre 5 de 1951  122.20332
## 3         25      25491    NOCAIMA                             1735   68.93823
## 4         25      25489    NIMAIMA    Ordenanza 30 de Julio de 1904   59.49982
## 5         25      25488       NILO                             1899  224.70968
##   MPIO_NANO   DPTO_CNMBR Shape_Leng  Shape_Area
## 0      2017 CUNDINAMARCA  0.5219117 0.004493674
## 1      2017 CUNDINAMARCA  1.2021630 0.032839624
## 2      2017 CUNDINAMARCA  0.4873683 0.009951825
## 3      2017 CUNDINAMARCA  0.3987992 0.005621814
## 4      2017 CUNDINAMARCA  0.4990524 0.004852617
## 5      2017 CUNDINAMARCA  0.8442993 0.018305852

Con este cuadro de código se puede observar un solo municipio. Para este caso será Pacho

pacho<-(munic[munic$MPIO_CNMBR=="PACHO",])
plot(pacho, main="Pacho", axes=TRUE)
plot(munic, add=TRUE) 
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.6))

Este código funciona para descargar los datos de elevación. Lo más recomendable es dejarlo como comentario para que cada vez que se lea el código completo este no se cargue de nuevo

elevation<-get_elev_raster(pacho, 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
## 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

Observar las características del objeto “elevacion”

elevation
## class      : RasterLayer 
## dimensions : 1034, 1033, 1068122  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.545, -71.70425, 2.796526, 5.629686  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -1.570384, 4150.449  (min, max)

A partir de esta información se pueden observar diferentes parámetros: El tamaño del archivo en pixeles La resolución o tamaño del pixel (recordar que un grado decimal corresponde a aproximadamente 111.11 km) El número de filas y columnas El número de celdas La extensión del ráster(este valor estará en las mimsas unidades de coordendas que el sistema de refrencia de coordenadas del ráster) El sistema coordenado de referencia, en este caso el datum de WGS84

Ahora se va a visualizar el DEM

plot(elevation, main="DEM de Pacho (metros)")
plot(pacho, add=TRUE)

3. Recortar datos de elevación para que coincidan con el tamaño del área de estudio

Como se pudo observar, el DEM realizado abarca mucho mayor área que la que es de interés. El siguiente cuadro de código tiene las instrucciones para recortar los datos de elevación correspondientes a Pacho

elev_crop = crop(elevation, pacho)
plot(elev_crop, main="DEM recortado")
plot(pacho, add=TRUE)

Obsevar las propiedades del objeto “elev_crop”

elev_crop
## class      : RasterLayer 
## dimensions : 94, 96, 9024  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.303, -74.039, 5.037846, 5.295406  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1013.265, 3738.211  (min, max)

4. Reproyectar los datos de elevación

Cuando se realizan DEM’s es recomendable usar coordenadas de mapa en lugar de coordendas geográficas. Esto se debe a que en las coordenadas geográficas las unidades de dimensiones horizonatles son grados decimales pero la unidad de dimesión vertical es en metros.

Para encontrar la proyección se puede buscar en epsg.io la proyección de la zona MAGNA Colombia Bogotá. Se necesita la información de dicha referencia en formato PROJ.4 (el usado para las bibliotecas sp y ráster)

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"

Ya al haber cargado la referencia se pueden reproyectar los datos de elevación en el sistema de coordenadas geográficas WGS84 en la Zonaa MAGNA Colombia Bogotá

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
res(pr3) <- 100
rep_elev <- projectRaster(elev_crop, pr3)

Observar las propiedades del objeto “rep_elev”

rep_elev
## class      : RasterLayer 
## dimensions : 285, 293, 83505  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 974994.6, 1004295, 1048824, 1077324  (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     : 1010.015, 3740.806  (min, max)

Ahora se puede realizar la reproyección de los dotados de elevación els sistema de coordenadas geográficas WGS84 en la zona MAGNA Colombia Bogotá

(rep_pacho = spTransform(pacho,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      : 975079.9, 1004237, 1048792, 1077338  (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       :         25,      25513,      PACHO,       1604, 402.62677731,      2017, CUNDINAMARCA, 1.20216295441, 0.0328396243456

Visualizar el DEM reproyectado.

plot(rep_elev, main="Reproyección del DEM")
plot(rep_pacho, add=TRUE)

5. Estádistica básica aplicada a datos de elevación

Este histograma permite tener una primera visión de la heterogeneidad de los datos de elevación en Pacho

hist(rep_elev, main="Histograma de los datos de elevación")

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))
##   metricas   valores
## 1     mean 2179.6584
## 2      min 1010.0148
## 3      max 3740.8058
## 4      std  664.3849

6. Obtención de variables geomórficas

Primero, calcular 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 debe vizualizar el DEM, en este caso se usó la paleta de colores “terrain.colors”

plot(rep_elev,main="DEM de Pacho [metros]", col=terrain.colors(25,alpha=0.7))

Visualización de la pendiente con la paleta de color “topo.colors”

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

Visualización de la orientación en grados usando la paleta de color “rainbow”

plot(aspect,main="Orientación de Pacho [grados]", col=rainbow(25,alpha=0.7))

Visualización de dos variables combinadas: elevación y sombreado

plot(hill,
        col=grey(1:100/100),  
        legend=FALSE,         
        main="DEM para Pacho",
        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” permite producir visualizaciones de datos 2D y 3D en R. Este paquete 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 que permiten generar mapas topográficos 2D y 3D. Además de los mapas, permite transfromar un ggplot2 en visualizaciones 3D

#install.packages("rayshader")
library(rayshader)
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()

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

A partir de este código se puede crear un ráster de sombreado de un DEM sugerido por 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/JUAN PABLO/Desktop/9pvbHjN.jpg")
out = getv(map, aspect, slope)
plotRGB(out, main = "Representación de la topografía de Pacho, Cundinamarca")

Esta visualización permite percibir que Pacho es una zona montañosa, más hacia la zona noroccidental.