Datos de elevación municipal en R

1. INTRODUCCIÓN

En el presente cuaderno de R se explicará cómo ilustrar distintas funcionalidades que permiten obtener, procesar y visualizar modelos digitales de elevación (DEM) en R. Los datos de elevación son comúnmente empleados para visualizaciones, hidrología y modelado ecológico, entre otras aplicaciones. Estos datos en R están disponibles a través de funciones en muchos paquetes o requiere acceso local a los datos debido a que no ha tenido una sola interfaz. Lo anterior ya no es necesario, ya que existe una variedad de API que brinda acceso programático a los datos de elevación, es el caso del paquete elevatr que estandariza el acceso a estos datos desde las API web y el cual se explorará en este cuaderno.

Para realizar la demostración trabajaremos con los datos del municipio de Rionegro del departamento de Antioquía. Antes de empezar es necesario instalar algunas librerías.

#Decomentar el siguiente código para instalar las librerías
#install.packages("rgdal")
#install.packages("raster")
#install.packages("elevatr")
#install.packages("rasterVis")
#install.packages("rgl")

library(rasterVis)

library(raster)

library(rgdal)

library(elevatr)

2. Obtención de datos de elevación ráster

Existen diferentes fuentes de modelos digitales de elevación, como Shuttle Radar Topography Mission (SRTM), USGS National Elevation Dataset (NED), Global DEM (GDEM), entre otras. Mapzen combinó varias de estas fuentes para crear un producto de elevación que utiliza los mejores datos disponibles para una región determinada a un nivel de zoom determinado. Estos datos están actualmente accesibles a través de Terrain Tiles en Amazon Web Services (AWS).

Usando get_elev_raster () para acceder a Terrain Tiles en AWS.

La entrada para get_elev_raster () es un data.frame con ubicaciones x (longitud) e y (latitud), un objeto sp o ráster y src debe establecerse en “mapzen” (este es el valor predeterminado), devolviendo un RasterLayer de los mosaicos que se superponen al cuadro delimitador de la entrada.

No existe diferencia en el uso de datos de entrada sp y ráster. El marco de datos requiere un prj. Se mostrarán ejemplos usando un SpatialPolygonsDataFrame. El nivel de zoom (z) será de 8 para evitar problemas de reinicio del programa al usar valores más altos.

Para empezar, debemos cargar un shapefile de los municipios del departamento de Antioquía obtenidos del geoportal del DANE.

Revisemos el contenido de la carpeta:

list.files("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"    

A través de una función del paquete ráster leamos el shapefile

(mpio <-  shapefile("ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 125 
extent      : -77.12783, -73.88128, 5.418558, 8.873974  (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  :         05,      05001,  ABEJORRAL,                                1541,   15.83597275,      2017,  ANTIOQUIA, 0.172962481259, 0.0012928296672 
max values  :         05,      05895,   ZARAGOZA, Ordenanza 9 de Diciembre 15 de 1983, 2959.36315089,      2017,  ANTIOQUIA,  6.61177822771,   0.23894871119 

mpio

class       : SpatialPolygonsDataFrame 
features    : 125 
extent      : -77.12783, -73.88128, 5.418558, 8.873974  (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  :         05,      05001,  ABEJORRAL,                                1541,   15.83597275,      2017,  ANTIOQUIA, 0.172962481259, 0.0012928296672 
max values  :         05,      05895,   ZARAGOZA, Ordenanza 9 de Diciembre 15 de 1983, 2959.36315089,      2017,  ANTIOQUIA,  6.61177822771,   0.23894871119 

Observemos las primeras columnas

head(mpio)

Ahora, seleccionamos un municipio, en este caso Rionegro

rionegro <- mpio[mpio$MPIO_CNMBR=="RIONEGRO",]
plot(rionegro, main="Rionegro", axes=TRUE)
plot(mpio, add=TRUE)

invisible(text(coordinates(mpio), labels=as.character(mpio$MPIO_CNMBR), cex=0.9))

Para descargar los datos de elevación descomente y ejecute el siguiente código.

#elevation <- get_elev_raster(rionegro, z = 8)

Revisemos que hay en el objeto: elevation

elevation

class      : RasterLayer 
dimensions : 1035, 1033, 1069155  (nrow, ncol, ncell)
resolution : 0.00275, 0.00273  (x, y)
extent     : -75.95125, -73.1105, 4.201768, 7.027318  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -430.9265, 5255.14  (min, max)

A tener en cuenta: * dimensions: “tamaño” del archivo en píxeles, * nrow, ncol: número de filas y columnas en los datos, * ncells: número total de píxeles (celdas) que componen el ráster, * Resolution: tamaño de cada píxel (grados, para este caso) (1 grado decimal representa ~111,11 km), * extent: extensión espacial del ráster (mismas unidades de coordenadas que el sistema de referencia de coordenadas del ráster), * crs: cadena del sistema de referencia de coordenadas para el ráster (para este caso es WGS 84).

plot(elevation, main="DEM de Rionegro (metros)")
plot(rionegro, add=TRUE)

3. Ajustar datos de elevación según la extensión del área de estudio

Mediante el siguiente código se podra cubrir solamente la extensión del municipio seleccionado, en este caso Rionegro

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

elev_crop = crop(elevation, rionegro)
plot(elev_crop, main="DEM de Rionegro ajustado")
plot(rionegro, add=TRUE)

Revisemos el nuevo objeto

elev_crop

class      : RasterLayer 
dimensions : 64, 60, 3840  (nrow, ncol, ncell)
resolution : 0.00275, 0.00273  (x, y)
extent     : -75.4865, -75.3215, 6.058168, 6.232888  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 2070.517, 2692.006  (min, max)

4. Reproyectar los datos de elevación

Cuando trabajamos con DEM es aconsejable utilizar coordenadas de mapa en lugar de coordenadas geográficas, debido a que en las coordenadas geográficas las unidades de dimensiones horizontales son grados decimales, pero la unidad de dimensión vertical son metros.

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

pr3 <- projectExtent(elev_crop, spatialref)
res(pr3) <- 100
rep_elev <- projectRaster(elev_crop, pr3)

Observemos que hay en el objeto rep_elev

rep_elev

class      : RasterLayer 
dimensions : 194, 183, 35502  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 844006.4, 862306.4, 1161800, 1181200  (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     : 2070.828, 2695.198  (min, max)

Ahora, volvamos a proyectar el SpatialPolygonsDataFrame que representa el municipio de Rionegro:

(rep_rionegro = spTransform(rionegro,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 844111.8, 862244.7, 1161773, 1181278  (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       :         05,      05615,   RIONEGRO,       1783, 195.9401959,      2017,  ANTIOQUIA, 0.890721733569, 0.0159994361869 

Visualización del DEM reproyectado

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

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

Realizar una exploración basica de los datos de elevación de Rionegro como su distribución, media, desviaciones.

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

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

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)

Visualización del DEM

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

Visualización de la pendiente

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

Visualización de la orientación

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

Visualización combinada del sombreado y el DEM

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

plot(rep_rionegro, add=TRUE)

7. Mapeo de datos de elevación con rayshader

Para producir visualizaciones de datos 2D y 3D en R se emplea la librería rayshader. 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 mapas topográficos 2D y 3D . Además de los mapas, rayshader también permite al usuario traducir objetos ggplot2 en visualizaciones de datos en 3D.

#install.packages("rayshader")

library(rayshader)

elmat = raster_to_matrix(rep_elev)

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

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

Calculating Surface Normal [------------------------------------] ETA:  0s
Calculating Surface Normal [====================================] ETA:  0s
                                                                          

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

Calculating Surface Normal [------------------------------------] ETA:  0s
Calculating Surface Normal [====================================] ETA:  0s
                                                                          

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

Calculating Surface Normal [------------------------------------] ETA:  0s
Calculating Surface Normal [====================================] ETA:  0s
                                                                          

Raytracing [----------------------------------------------------] ETA:  0s
Raytracing [====================================================] ETA:  0s
                                                                          

8. Otra forma de visualización

Antes de inciar instale la librería jpeg

#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("9pvbHjN.jpg")
out = getv(map, aspect, slope)

no non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Infno non-missing arguments to min; returning Infno non-missing arguments to max; returning -Inf

plotRGB(out)