Datos de Elevacion en R

1. Introduccion

En este notebook se utilizara el paquete estadistico “elevatr” como una fuente de acceso de datos de elevacion y se realizaran diferentes mapas de elevacion digital del municipio de Sopo perteneciente al departamento de Cundinamarca, debido a que no se puede representar el mapa de su capital Bogota, ya que esta no le pertenece al departamento de Cundinamarca.

2. Elevacion

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, la elevación de una ubicación geográfica es la altura por encima o por debajo de un punto de referencia fijo, más comúnmente una referencia geoide, un modelo matemático del nivel del mar de la Tierra como una superficie equipotencial, tener acceso a esta informacion en R no ha tenido una sola interfaz, está disponible a través de funciones en muchos paquetes, como por ejemplo el paquete elevatr el cual se escribió para estandarizar el acceso a los datos de elevación desde las API web.

Para esto primero vamos a instalar las siguientes librerias

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

library(rasterVis)

library(raster)
library(rgl)
library(rgdal)

library(elevatr)

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

Un ráster consta de una matriz de celdas (o píxeles) organizadas en filas y columnas (o una cuadrícula) en la que cada celda contiene un valor que representa información, como la temperatura, altura. Los rásteres son fotografías de aéreas digitales, imágenes de satélite, imágenes digitales o incluso mapas escaneados.

Lo primero que realizaremos sera traer los datos del departamento de Cundinamarca descargados de la base de datos del geoportal del DANE y utilizaremos la funcion “list.files()” la cual genera un vector de caracteres de los nombres de archivos o directorios.

list.files("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Cundinamarca/DANE/25_CUNDINAMARCA/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 archivo de forma usando una función provista por el paquete ráster, designando a un objeto llamado “munic”

# Ejemplo de un datafream de un poligono espacial
(munic <-  shapefile("C:/Users/mario/Documents/Ing. Agronomica/9 Semestre/Geomatica/Cundinamarca/DANE/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 +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  :         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

¿Cuáles son los atributos del objeto munic? Para responder esto podemos utilizar la funcion “head()”, para ver el encabezado del objeto.

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

Ahora seleccionemos solamente la ciudad capital del departamento de Cundinamarca, para este caso aunque en Bogota se encuentra el Distrito Capital no le pertenece a Cundinamarca, seleccionaremos al municipio de “SOPO”

Si no recuerda como estan escritos los datos de los municipios de su departamento puede solamente correr como una linea de codigo su objeto “munic” y utizando el signo pesos para traer solo los datos de los municipios de la base de datos, que en este caso sale como “MPIO_CNMBR”.

#munic$MPIO_CNMBR

## make sure you review the documentation to  understand the following lines of code
## https://www.neonscience.org/dc-shapefile-attributes-r
sopo <- munic[munic$MPIO_CNMBR=="SOPO",]
plot(sopo, main="Cundinamarca", axes=TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

Ahora utilizaremos la siguiente línea para descargar los datos de elevación

elevation <- get_elev_raster(sopo, z = 8)

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 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : -1.570384, 4150.449  (min, max)

Tenga en cuenta algunas cosas sobre este DEM:

Dimensiones: Es 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 decinales en este caso). Recordemos que 1 grado decimal representa aproximadamente 111,11 km. Extensión: La extensión espacial de el ráster. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del ráster. crs: Es el sistema de referencia de coordenadas para el ráster. Este ráster está en coordenadas geográficas con un dato de WGS 84.

plot(elevation, main="DEM descargado (metros)")
plot(sopo, add=TRUE)

4. Recorte 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 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.

#writeRaster(elevation, filename="./layer.tif", datatype='INT4S', overwrite=TRUE)

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

# review the raster documentation to understand how the crop function works
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/crop
elev_crop = crop(elevation, sopo)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(sopo, add=TRUE)

elev_crop

## class      : RasterLayer 
## dimensions : 59, 35, 2065  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.0115, -73.91525, 4.813166, 4.974826  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 2546.713, 3213.353  (min, max)

5. Reproyectar los datos de elevación

Cuando se trabaja con DEM, siempre es una buena idea usar 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 es metros. Vuelva a proyectar los datos de elevación.

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

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

rep_elev

## class      : RasterLayer 
## dimensions : 179, 107, 19153  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1007320, 1018020, 1023971, 1041871  (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     : 2546.998, 3206.572  (min, max)

Ahora, reproyectemos el datafream del poligono espacial, el cual representa la capital de nuestro departamento de Cundinamarca

(rep_sopo = spTransform(sopo,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1007370, 1017935, 1023919, 1041782  (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       :         25,      25758,       SOPO,       1653, 111.08048482,      2017, CUNDINAMARCA, 0.59825817844, 0.00905630487218

plot(rep_elev, main="Reproyeccion digital del modelo de Elevacion")
plot(rep_sopo, add=TRUE)

Para evitar inconvenientes, guardemos nuestra DEM

writeRaster(rep_elev, filename="./rep_sopo_elev.tif", datatype='INT4S', overwrite=TRUE)

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

Ahora realizaremos una exploración rápida de las estadísticas DEM:

# histograma
hist(rep_elev,main = "Datos de Elevacion de la reproyeccion digital", xlab = "Elevacion (m.s.n.m)", ylab = "Frecuencia", col = "lightblue", ylim=c(0,12000))

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

##   metricas  valores
## 1     mean 2685.659
## 2      min 2546.998
## 3      max 3206.572
## 4      std  161.765

7. Obtención de variables geomorfométricas.

Antes de continuar, se recomienda haber leído el capítulo de geomorfometría escrito por Victor Olaya para su libro libre en sistemas de información geográfica, se puede encontrar en este enlace https://volaya.github.io/libro-sig/chapters/Geomorfometria.html.

Primero, calcule la pendiente “slope”, el aspecto “aspect” y el sombreado “hill”:

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

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

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

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

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

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

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

Plot combinado

#plot DEM w/ hillshade
plot(hill,
        col=grey(1:100/100),  # crear una rampa de color de colores grises para sombreado
        legend=FALSE,         # sin leyenda, no nos importa el gris de la ladera
        main="DEM de Sopo",
        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) # metodo de color

8. 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 utiliza datos de elevación en una matriz base R y una combinación de trazado de rayos, mapeo de texturas esféricas, 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 3D.

Lo primero que realizaremos es instalar este paquete estadistico “rayshader”

#install.packages("rayshader")

library(rayshader)

Ahora debemos convertir el DEM en una matrix, lo cual realizaremos con la siguiente funcion “raster_to_matrix()”

elmat = raster_to_matrix(rep_elev)

# Utilizaremos una de las texturas incorporadas de la libreria rayshader
elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

#detect_water y add_water agrega 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 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 ()

9. Otra forma de visualización.

#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
}

Ahora debemos descargar la siguiente imagen de ejemplo del entorno del mapa en el siguiente enlace http://i.imgur.com/9pvbHjN.jpg, y utilizando la funcion “readJPEG()” leeremos este archivo.

map=readJPEG("./9pvbHjN.jpg")

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

# Plot de montañas
plotRGB(out, main = "Monatañas de Sopo")