Datos de elevación en R

Luisa Fernanda Carrión Ramírez

02/04/2020

Este cuaderno creado en R tiene como objetivo ilustrar algunas funcionalidades contenidas en este programa que son útiles para el procesamiento y visualización de unos modelos digitales de elevación. Para esta práctica se seleccionó el municipio de Yopal.

title: “introducción a elevart” Los datos de elevación pueden ser utilizados de diferentes maneras por su amplia gama de aplicaciones. Es por ello que anteriormente se requerían instalar varios paquetes aunque, en la actualidad se utilizan API’s que simplifican el proceso. El paquete elevatr fue programado con el propósito de estandarizar el acceso de los datos de elevación desde las API web. El paquete elevatr emplea mosaicos de terrenos de Amazon Web Services para acceder a los datos de elevación de los ráster. Para iniciar, se deben descargar los paquetes rgdal, raster, elevatr, rasterVis y rgl.

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

Luego, se corren los paquetes usando la función library()

library(rgdal)
library(raster)
library(elevatr)
library(rgl)
library(rasterVis)

En adición, se requieren los paquetes sp, lattice y latticeExtra

library(sp)
library(latticeExtra)
library(lattice)

Obtenga datos de elevación de ráster

Hay varias fuentes para los modelos de elevación digital y cada uno de estos DEM tiene ventajas y desventajas para 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 con un nivel de zoom determinado. Aunque cerró estos fueron compilados por Mapzen y pueden ser consultados en Amazon Web Services (AWS).

En este cuaderno, cargaremos se trabajará con un shapefile descargado del geoportal del DANE ,que contiene datos correspondientes al departamento de Casanare que se cargará con el siguiente código.

list.files("C:/Users/LUISA CARRION/Documents/Geomática/85_CASANARE/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"     "MGN_DPTO_POLITICO.dbf"     "MGN_DPTO_POLITICO.prj"    
 [4] "MGN_DPTO_POLITICO.sbn"     "MGN_DPTO_POLITICO.sbx"     "MGN_DPTO_POLITICO.shp"    
 [7] "MGN_DPTO_POLITICO.shp.xml" "MGN_DPTO_POLITICO.shx"     "MGN_MPIO_POLITICO.cpg"    
[10] "MGN_MPIO_POLITICO.dbf"     "MGN_MPIO_POLITICO.prj"     "MGN_MPIO_POLITICO.sbn"    
[13] "MGN_MPIO_POLITICO.sbx"     "MGN_MPIO_POLITICO.shp"     "MGN_MPIO_POLITICO.shp.xml"
[16] "MGN_MPIO_POLITICO.shx"    

Luego se lee el shapefile:

(munic <- shapefile("C:/Users/LUISA CARRION/Documents/Geomática/85_CASANARE/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 19 
extent      : -73.07777, -69.83591, 4.287476, 6.346111  (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  :         85,      85001,    AGUAZUL,                                  1748,  181.22675497,      2017,   CASANARE, 0.674577402276, 0.0147760191324 
max values  :         85,      85440,      YOPAL, Ordenanza 0038 del 8 de Julio de 1942, 12115.0422213,      2017,   CASANARE,  8.57006487338,  0.985967805522 

Por medio de la función head() vemos los archivos

head(munic)

Ahora, se selecciona únicamente la capital del departamento

yopal <- munic[munic$MPIO_CNMBR=="YOPAL",]
plot(yopal, main="Yopal", axes=TRUE)
plot(munic, add=TRUE)

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

Procedemos a crear el objeto elevacion.

elevacion <- get_elev_raster(yopal, z=8)

Mediante el siguiente comando se puede ver lo que compone el objeto elevacion:

elevacion

class      : RasterLayer 
dimensions : 1034, 1033, 1068122  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -73.13875, -70.298, 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     : -42.2496, 3841.442  (min, max)

En el objeto elevacion se tienen las dimensiones, número de filas y columnas en los datos, el número total de píxeles, la resolución cada píxel y la extensión espacial de la trama. Además, también se puede ver que el ráster esta en coordenadas geográficas WGS 84.

Para trazar el DEM se escribe el siguiente código:

plot(elevacion, main="DEM descargado (metros)")
plot(yopal, add=TRUE)

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

Como se pudo observar en la aterior gráfica, el área correpsondieente a Yopal es muy pequeña, por ello, debemos cortar el área que necesitamos. Sin embargo, es mejor guardar la capa sin cortar:

writeRaster(rep_elev, filename = "C:/Users/LUISA CARRION/Documents/Geomática/85_CASANARE/elevacion", dataType='INT4S', overwrite=TRUE)

argument "datatype" misspelled as "dataType"

class      : RasterLayer 
dimensions : 739, 716, 529124  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1165579, 1237179, 1034149, 1108049  (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/LUISA CARRION/Documents/Geomática/85_CASANARE/elevacion.grd 
names      : layer 
values     : 161, 2421  (min, max)

Por medio del siguiente código se recorta el departamento

elev_crop = crop(elevacion, yopal)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(yopal, add=TRUE)

Ahora, visualizamos el nuevo objeto

elev_crop

class      : RasterLayer 
dimensions : 243, 234, 56862  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -72.58325, -71.93975, 4.903586, 5.569406  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 160.9496, 2428.565  (min, max)

Reproyectar los datos de elevación

Cuando se trabaja con DEM, es preferible utilizar las coordenadas del lugar en donde se esta trabajando. Es por ello quese visitará (epsg.io) y se buscará la proyección “MAGNA- SIRGAS Colombia Bogota zone”. Necesitamos obtener la definición de esta referencia espacial en formato PROJ.4 ( Pues esta es la que se usa para las bibliotecas sp y raster ). Copiemamos el texto PROJ.4 y lo guardamos.

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"

pr3 <- projectExtent(elev_crop, spatialref)

res(pr3) <- 100

rep_elev <- projectRaster(elev_crop, pr3)

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

Posteriormente, se reproyectarán las coordenadas:

rep_elev

class      : RasterLayer 
dimensions : 739, 716, 529124  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1165579, 1237179, 1034149, 1108049  (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     : 160.9715, 2420.886  (min, max)

Ahora, debemos reproyectar el poligono de atributos de Yopal:

(rep_yopal = spTransform(yopal,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 1165549, 1237211, 1034289, 1107987  (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       :         85,      85001,      YOPAL, Ordenanza 0038 del 8 de Julio de 1942, 2482.9042923,      2017,   CASANARE, 3.15923255859, 0.202336772314 

Luego, se traza el objeto obtenido

plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_yopal, add=TRUE)

Ahora, se guarda el DEM

writeRaster(rep_elev, filename="C:/Users/LUISA CARRION/Documents/Geomática/85_CASANARE/elevacion.tif", datatype='INT4S', overwrite=TRUE)

Estadísticas básicas de los datos de elevación

Se exploran algunas estadísticas contenidas en el DEM empezando por el histograma:

hist(rep_elev)

Luego, miramos la elevación por metros:

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

Obtención de Variables geomorfomótricas

En esta sección se calcula 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)

Luego, se trazan la pendieente, aspecto y sombreado:

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

Grafiamos la pendiente:

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

Graficamos el aspecto:

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

Podemos cambiar la paleta de colores, en este caso, escogí “Viridis”

install.packages("viridis")

library(viridis)

Graficamos e aspecto con “Viridis”

plot(aspect,main="Aspecto de Yopal [Grados]", col=viridis(25,alpha=0.7))

Se puede hacer un plot combinado:

plot(hill,
        col=grey(1:100/100),  
        legend=FALSE,         
        main="DEM de Yopal",
        axes=FALSE)          

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.35), add=TRUE)

Mapeo de datos de elevación con rayshader.

La librería rayshader es un paquete que permite visualizar los datos en formato 2D y 3D. Instalamos y cargamos la libreri:

install.packages("rayshader")

Ahora, ase debe convertir el DEM en matríz:

library(rayshader)

Ahora, se le asigna una textura RGB a un sombreado:

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 716x739."

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

Con “detect_water” se pueden detectar cuerpos de agua.

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

Se puede agregar una capa de trazado de rayos:

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

Otra forma de visualización

Para llevarla a cabo, se necesita instalar el paquete “jpeg”

install.packages("jpeg")

library(jpeg)

Luego, se ejecuta el siguiente código:

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
}

Por ultimo, se escribe:

map=readJPEG("C:/Users/LUISA CARRION/Documents/Geomática/85_CASANARE/elevacion.jgp")

out = getv(map, aspect, slope)

plotRGB(out,main="Montañas de Yopal")