Datos de elevacion en R

1.Introducción

En este cuaderno ilustraremos algunas funciones para obtener, procesar y visualizar Modelos Digitales de Elevación en R. Los datos de elevación son usados en un amplio rango de aplicaciones, por ejemplo, hidrología y modelamiento ecológico. se trabajara con el Paquete elevatr el cual estandarizo el accesso a los datos de elevación desde las API web. En este caso para dar ejemplo de estas herramientas el trabajo tendra como foco el departamento de NBoyacá. primero instalar los paquetes necesarios:

install.packages("rgdal")
install.packages("elevatr")
install.packages("lattice")
install.packages("latticeExtra")
install.packages("rasterVis")
install.packages("rgl")

library(rgdal)

library(elevatr)
library(lattice)
library(latticeExtra)
library(rasterVis)
library(rgl)
library(raster)
library(sp)

2. Obtencion daros de elevación

Se debe tener en cuenta que hay varias fuentes de donde obtener los modelos de elevación digital. Con la entrada para getelevraster() para acceder a terrain tiles en AWS el cual es un data.frame con ubicaciones x para longitus e Y para latitud, un objeto sp o un objeto raster necesita ser la entrada y src debe establecerse en “mapzen” (este es un valor predeterminado).

El primer paso es cargar el archivo .shp del departamento de Boyacá el cual se descargo del geoportal del DANE. Para esto se usa la funcion *list.files(), podemos observar los archivos que hay dentro del directorio de trabajo:

list.files("C:/Users/David Perdomo/Desktop/Geomatica/15_BOYACA/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"                            
 [2] "MGN_DPTO_POLITICO.dbf"                            
 [3] "MGN_DPTO_POLITICO.prj"                            
 [4] "MGN_DPTO_POLITICO.sbn"                            
 [5] "MGN_DPTO_POLITICO.sbx"                            
 [6] "MGN_DPTO_POLITICO.shp"                            
 [7] "MGN_DPTO_POLITICO.shp.DG_EST54.1344.12652.sr.lock"
 [8] "MGN_DPTO_POLITICO.shp.xml"                        
 [9] "MGN_DPTO_POLITICO.shx"                            
[10] "MGN_MPIO_POLITICO.cpg"                            
[11] "MGN_MPIO_POLITICO.dbf"                            
[12] "MGN_MPIO_POLITICO.prj"                            
[13] "MGN_MPIO_POLITICO.sbn"                            
[14] "MGN_MPIO_POLITICO.sbx"                            
[15] "MGN_MPIO_POLITICO.shp"                            
[16] "MGN_MPIO_POLITICO.shp.DG_EST54.1344.12652.sr.lock"
[17] "MGN_MPIO_POLITICO.shp.xml"                        
[18] "MGN_MPIO_POLITICO.shx"                            

(munic <- shapefile("C:/Users/David Perdomo/Desktop/Geomatica/15_BOYACA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 123 
extent      : -74.66496, -71.94885, 4.655196, 7.055557  (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  :         15,      15001,    ALMEIDA,                               1537,   25.35271459,      2017,    BOYACÁ, 0.204497978002, 0.00206924119091 
max values  :         15,      15897, ZETAQUIRÁ, Ordenanza 8 de Diciembre 4 de 1965, 1513.60104233,      2017,    BOYACÁ,  2.40391520946,   0.123609862522 

head(munic)

Seleccionemos solo la ciudad capital de este departamento.

tunja <- munic[munic$MPIO_CNMBR=="TUNJA",]

plot(tunja,  axes=TRUE, col="green")
plot(munic, add=TRUE)

text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR))

Descomente la siguiente línea para ejecutar el código y descargar los datos de elevación. Una vez que se haya ejecutado el código, comente el fragmento.

elevation <- get_elev_raster(tunja, z = 12)

elevation

class      : RasterLayer 
dimensions : 1545, 1543, 2383935  (nrow, ncol, ncell)
resolution : 0.000172, 0.000171  (x, y)
extent     : -73.47742, -73.21203, 5.352646, 5.616841  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 2026.88, 3543.023  (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 Tunja [metros]")
plot(tunja, add=TRUE)

text(coordinates(tunja), labels=as.character(tunja$MPIO_CNMBR))

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

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

text(coordinates(tunja), labels=as.character(tunja$MPIO_CNMBR))

elev_crop

class      : RasterLayer 
dimensions : 805, 834, 671370  (nrow, ncol, ncell)
resolution : 0.000172, 0.000171  (x, y)
extent     : -73.44732, -73.30387, 5.446354, 5.584009  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : layer 
values     : 2231.158, 3262.087  (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)
res(pr3) <- 20
rep_elev <- projectRaster(elev_crop, pr3)

rep_elev

class      : RasterLayer 
dimensions : 762, 796, 606552  (nrow, ncol, ncell)
resolution : 20, 20  (x, y)
extent     : 1069823, 1085743, 1094051, 1109291  (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     : 2228.55, 3261.334  (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_tunja = spTransform(tunja,spatialref))

Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
 but +towgs84= values preserved

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 1069842, 1085712, 1094054, 1109295  (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       :         15,      15001,      TUNJA,       1541, 119.68956961,      2017,    BOYACÁ, 0.572374355128, 0.00976630124258 

Visualizar el DEM reproyectado.

plot(rep_elev, main="Proyección del DEM")
plot(rep_tunja, add=TRUE)

text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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

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 Tunja [metros]", col=terrain.colors(25,alpha=0.8))
plot(rep_tunja, add=TRUE, lwd=2)

text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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

plot(slope,main="pendiente para Tunja [grados]", col=topo.colors(25,alpha=0.7))
plot(rep_tunja, add=TRUE, lwd=2)

text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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

plot(aspect,main="orientación de Tunja [grados]", col=rainbow(25,alpha=0.7))
plot(rep_tunja, add=TRUE, lwd=2)

text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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

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

plot(rep_tunja, add=TRUE, lwd=2)
text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

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)

[1] "Dimensions of matrix are: 796x762."

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/David Perdomo/Desktop/Geomatica/9pvbHjN.jpg")
out = getv(map, aspect, slope)

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

plotRGB(out, main = "Montañas en Tunja")
polygons(rep_tunja)

class       : SpatialPolygons 
features    : 1 
extent      : 1069842, 1085712, 1094054, 1109295  (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 

text(coordinates(rep_tunja), labels=as.character(rep_tunja$MPIO_CNMBR))

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Matrix products: default

locale:
[1] LC_COLLATE=Spanish_Colombia.1252 
[2] LC_CTYPE=Spanish_Colombia.1252   
[3] LC_MONETARY=Spanish_Colombia.1252
[4] LC_NUMERIC=C                     
[5] LC_TIME=Spanish_Colombia.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
[1] jpeg_0.1-8.1  raster_3.3-13 sp_1.4-2     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        rstudioapi_0.11   knitr_1.29       
 [4] magrittr_1.5      tidyselect_1.1.0  lattice_0.20-41  
 [7] R6_2.4.1          rlang_0.4.7       fansi_0.4.1      
[10] blob_1.2.1        dplyr_1.0.2       tools_4.0.2      
[13] grid_4.0.2        xfun_0.17         cli_2.0.2        
[16] DBI_1.1.0         dbplyr_1.4.4      htmltools_0.5.0  
[19] ellipsis_0.3.1    assertthat_0.2.1  digest_0.6.25    
[22] tibble_3.0.3      lifecycle_0.2.0   crayon_1.3.4     
[25] purrr_0.3.4       codetools_0.2-16  htmlwidgets_1.5.1
[28] vctrs_0.3.4       glue_1.4.2        compiler_4.0.2   
[31] pillar_1.4.6      generics_0.0.2    pkgconfig_2.0.3