Datos de elevacion en R

Angie Alejandra Juyo Buitrago

2 de abril de 2020

1.Introducción

Este es un cuaderno R donde se Ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación en R. Es realizado para el curso de geomatica básica de la Universidad Nacional de Colombia, Sede Bogotá

2.El paquete elevatr

El paquete elevatr fue escrito para estandarizar el acceso a los datos de elevación desde las APIs (interfaz de programaci?n de aplicaciones) de la web.Pues , de forma anterior el acceso a estos datos en R no tenia una sola interfaz, pues se precisaba el uso de funciones en muchos paquetes, o requieria el acceso local a los datos.

Los datos de elevación se utilizan para una amplia gama de aplicaciones, entre ellas, por ejemplo, la visualización, la hidrología y la modelización ecológica.

Para acceder a los datos de elevación rasterizados (por ejemplo, un DEM) el paquete elevatr utiliza Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.

Primero , instalamos el paquete elevatr junto a los otros paquetes que van a ser utilizados en el desarrollo de este cuaderno

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

Luego , procedemos a cargar las librerias necesarias

library(rgdal)

rgdal: version: 1.4-8, (SVN revision 845)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
 Path to GDAL shared files: C:/Users/LUISA CARRION/Documents/R/win-library/3.6/rgdal/gdal
 GDAL binary built with GEOS: TRUE 
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: C:/Users/LUISA CARRION/Documents/R/win-library/3.6/rgdal/proj
 Linking to sp version: 1.4-1 

library(raster)
library(elevatr)
library(rasterVis)

Loading required package: lattice
Loading required package: latticeExtra

library(rgl)

3. Obtener los datos de elevación de Raster

Hay varias fuentes para los modelos de elevación digital como la Shuttle Radar Topography Mission (SRTM),el USGS National Elevation Dataset (NED), Global DEM (GDEM), y otras.

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 que utiliza los mejores datos de elevaci?n disponibles para una región determinada a un nivel de zoom determinado. Aunque cerrado, estos datos compilados por Mapzen son actualmente accesibles a través de Amazon Web Services (AWS).

Para comenzar con la exploración de estas funciones,primero, necesitamos cargar un archivo de forma que represente los municipios de nuestro departamento.

En este cuaderno, cargaremos un shapefile descargado del geoportal del DANE ,que contiene datos correspondientes al departamento de Vichada

Primero, podemos ver que contiene la carpeta donde fueron descargados los datos

list.files("C:/Users/LUISA CARRION/Documents/99_VICHADA/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 leemos nuestro shapefile con ayuda de una de las funciones del paquete ráster

(munic <-  shapefile("C:/Users/LUISA CARRION/Documents/99_VICHADA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

class       : SpatialPolygonsDataFrame 
features    : 4 
extent      : -71.07793, -67.4098, 2.737109, 6.324317  (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  :         99,      99001,       CUMARIBO,       Decreto 1594 de Ago 5  de 1974, 3898.56891769,      2017,    VICHADA, 3.29670807195, 0.316732435778 
max values  :         99,      99773, SANTA ROSALÍA, Ordenanza 66 de Noviembre 22 de 1996, 65599.7022767,      2017,    VICHADA,  18.794382661,   5.3085802966 

Una vez leidos ,podemos ver , cuales son los atributos de nuestro objeto, con ayuda de la función “head”

head(munic)

Luego, seleccionamos solo la capital de nuestro departamento

Pto_Carreno <- munic[munic$MPIO_CNMBR=="PUERTO CARREÑO",]
plot(Pto_Carreno, main="Puerto Carreño", axes=TRUE)
plot(munic, add=TRUE)

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

3.1Usar get_elev_raster() para acceder a los Terrain Tiles en AWS.

La entrada para get_elev_raster() es un data.frame con ubicaciones x(longitud) e y(latitud) como las dos primeras columnas, cualquier objeto sp o cualquier objeto r?ster y devuelve un RasterLayer de los mosaicos que se superponen al cuadro delimitador de la entrada. Si se recuperan m?ltiples mosaicos, la salida resultante sar? una capa r?ster fusionada.

Como se mencionó, un data frame con columnas xey, un objeto sp o un raster debe ser la entrada y src debe establecerse en “mapzen” (este es el valor predeterminado).

No hay diferencia en el uso de las sp y raster de entrada de tipos de datos. El marco de datos requiere a prj.

El nivel de zoom determina la resolución del raster de salida y su valor predeterminado es 9 (Pues tiene una buena resolución y no es mucho el tiempo de descarga)

elevation <- get_elev_raster(Pto_Carreno, z = 8)

Luego,podemos ver que está contenido en el objeto “Elevacion”, obtendremos información como dimensiones: el “tamaño” del archivo en píxeles, numero de filas, numero de columnas ,el número total de celdas que componen el raster inclusive podemos obtener información de la referencia de sistema de coordenadas para el ráster

elevation

class      : RasterLayer 
dimensions : 2570, 2567, 6597190  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -71.7325, -64.67325, -0.01443207, 7.027368  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -419.5286, 2829.102  (min, max)

Podemos trazar los datos de elevacion del objeto de la siguiente manera

plot(elevation, main=" DEM [metros]")

plot(Pto_Carreno, add=TRUE)

4.Recortar los datos de elevación para que coincidan con el área de estudio

Como observamos, 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,antes de recortar es una buena idea guardar el DEM para el futuro.

writeRaster(elevation, filename = "Aqui iria la ruta de mi computador", dataType='INT4S', overwrite=TRUE)

argument "datatype" misspelled as "dataType"

class      : RasterLayer 
dimensions : 2570, 2567, 6597190  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -71.7325, -64.67325, -0.01443207, 7.027368  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : C:/Users/LUISA CARRION/Downloads/Aqui iria la ruta de mi computador.grd 
names      : layer 
values     : -420, 2829  (min, max)

Y recortamos nuestros datos correspondientes al departamento de Vichada

elev_crop = crop(elevation, Pto_Carreno)
plot(elev_crop, main="Modelo de elevación recortado")

plot(Pto_Carreno, add=TRUE)

Luego, vemos el nuevo objeto

elev_crop

class      : RasterLayer 
dimensions : 1309, 1334, 1746206  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -71.078, -67.4095, 2.736528, 6.323188  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : -320.8194, 1803.522  (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 las coordenadas geográficas, las unidades de dimensiones horizontales son grados decimales pero la unidad de dimensi?n vertical son los metros.

Podemos ir a (epsg.io) y 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" 

Ahora, podemos reproyectar los datos de elevación de las coordenadas geográficas WGS84 en MAGNA Colombia Bogota zone

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

rep_elev

class      : RasterLayer 
dimensions : 4010, 4111, 16485110  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1332035, 1743135, 794733, 1195733  (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     : -310.4869, 1802.343  (min, max)

Ahora, se va a reproyectar el SpatialPolygonsDataFrame Poligono con atributos que representa a Puerto Carreño:

(rep_Pto_Carreno = spTransform(Pto_Carreno,spatialref))

Después, trazamos el objeto obtenido

plot(rep_elev, main = "DEM reproyectado")

plot(rep_Pto_Carreno, add= TRUE)

Error in plot(rep_Pto_Carreno, add = TRUE) : 
  objeto 'rep_Pto_Carreno' no encontrado

Luego , guardamos nuestro DEM

writeRaster(rep_elev, filename = "Ac? va la ruta de archivo", dataType='INT4S', overwrite=TRUE)

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

Se hace una exploaración básica de las estadísticas del DEM

Primero,realizamos un histograma

hist(rep_elev)

Luego realizamos un dataframe con estadásticas de elevación en 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

Primero, calcule la pendiente, el aspecto y el sombreado

La pendiente puede ser calculada con la función terrain de la libreria Raster, mientras que la función hillshade, calcula la sombra de una montaña o una colina a partir de la pendiente ,el aspecto el ángulo de elevación del sol en grados y “azimuth”" el ángulo dirección del sol en grados

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

Luego, trazamos elevación , pendiente y aspecto , fijandonos en las paletas de color que son usadas para cada uno de los ploteos

Elevación

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

Pendiente

plot(slope,main="Pendiente para Puerto Carreño [Grados]", col=topo.colors(25,alpha=0.7))

Aspecto

plot(aspect,main="Aspecto para Puerto Carreño [Grados]", col=rainbow(25,alpha=0.7))

También ,podemos hacer un plot combinado

plot (hill, col= grey(1:100/100), legend= FALSE, main = "DEM para Puerto Carreño",  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 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 de R y una combinación de trazado de rayos, mapeo de textura esf?rica, superposiciones y oclusi?n ambiental para generar mapas topográficos en 2D y 3D.

Además de los mapas, rayshader también permite al usuario traducir objetos ggplot2 en visualizaciones de datos 3D.

Instalamos y cargamos la libreria

#install.packages("rayshader")
library(rayshader)

Primero , convertimos nuestro DEM , en una matriz

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 4111x4010."

Podemos asignar o incorporar una textura RGB a un sombreado mediante un mapeo esférico con la función “Sphere_Shade”

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

Podemos usar “detect_water” para detectar cuerpos de agua ( Rayshader usa un algoritmo de inundación para este fin ) con la función para agregar un color de agua al mapa

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

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

8.Otra forma de visualización

Un experto en R, sugiere otra forma de visualización aqui

#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("ACA IRIA UNA RUTA DE ARCHIVO") 

out = getv(map, aspect, slope)

plotRGB(out, main = "Montañas en Vichada")