Juan Pablo Cañon Sosa

18/09/2020

1. Porque Usar este notebook?

Mediante este cuaderno de R usaremos diversas funcionalidades para obtener, procesar y visualizar \(modelos\) \(de\) \(elevacion\) \(digital\) en R. En esta ocasion Trabajaremos con los datos correspondientes a la ciudad de Cali, la cual es capital del departamento del Valle del Cauca.

2. Introducion a \(elevatr\)

Los datos de elevacion son usados para diversas aplicaciones tales como: visualizacion, hidrologia y modelacion ecologica entre otros. Obtener acceso a estos datos en R no ha tenido una sola interfaz, está disponible a través de funciones en muchos paquetes o requiere acceso local a los datos. Esto ya no es necesario ya que ahora existe una variedad de API que proporcionan acceso programático a los datos de elevación. El paquete elevatr se escribió para estandarizar el acceso a datos de elevación de API web. Para acceder a los datos de elevación ráster (por ejemplo, un DEM), el paquete elevatr utiliza los mosaicos de terreno de Amazon Web Services.

Nosotros exploraráremos esta funcionalidad en este cuaderno.

Para ello instalamos las siguientes librerias:

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

Llamamos las librerias anteriormente instaladas:

library(lattice)
library(latticeExtra)
library(rasterVis)
## Loading required package: raster
## Loading required package: sp
library(rgl)
library(rgdal)
## rgdal: version: 1.5-16, (SVN revision 1050)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/FLCS/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/FLCS/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
library(raster)
library(elevatr)

3. Get Datos de Raster de elevacion:

Hay varias fuentes de modelos digitales de elevación, como Shuttle Radar Topography Mission (SRTM), el Conjunto de datos de elevación nacional (NED) del USGS, DEM global (GDEM) y otros. Cada uno de estos DEMs tiene pros y contras para su uso. Antes de su cierre en enero de 2018, Mapzen combinó varias de estas fuentes para crear una síntesis Producto de elevación que utiliza los mejores datos de elevación disponibles para una región determinada a un nivel de zoom determinado. A pesar de que fue cerrado, estos datos compilados por Mapzen son actualmente accesibles a través de Terrain Tiles en Amazon Web Services (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 varios mosaicos, la salida resultante es una capa ráster combinada.

Revisamo el directorio de trabajo usando el comando getwd():

getwd()
## [1] "C:/Users/FLCS/Documents/R/Geomatica Basica/GB/Notebook"

configuramos el directorio de trabajo usando el comando setwd():

setwd("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/datos")

Usando get_elev_raster() para acceder a los mosaicos de terreno en AWS.

Como se mencionó, un marco de datos con columnas x e y, un objeto sp o un objeto ráster deben ser la entrada y el src debe establecerse en “mapzen” (este es el valor predeterminado).

No hay diferencia en el uso de los tipos de datos de entrada sp y ráster. El marco de datos requiere un prj. Mostraré ejemplos utilizando un SpatialPolygonsDataFrame. El nivel de zoom (z) predeterminado es 9 (una compensación entre resolución y tiempo para descargar), pero no pude obtener nada con este zoom.

Primero, necesitamos cargar un shapefile que represente a los municipios de nuestro departamento. En este cuaderno, cargaremos el shapefile descargado del geoportal del DANE.

Usando la funcion list.files(), vemos los archivos que hay dentro de nuestro directorio de trabajo:

list.files("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/datos/76_VALLE_DEL_CAUCA/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, vamos a leer el shapefile usando una funcion \(shapefile\) provista por el paquete \(raster\) y almacenamos los datos en un objeto que llamaremos munic:

# SpatialPolygonsDataFrame ejemplo
(munic <- shapefile("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/datos/76_VALLE_DEL_CAUCA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 42 
## extent      : -77.54977, -75.70724, 3.091239, 5.047394  (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  :         76,      76001,    ALCALÁ,                          1536,   41.86090736,      2017, VALLE DEL CAUCA, 0.453826056161, 0.00340901158098 
## max values  :         76,      76895,     ZARZAL, Ordenanza 9 de Diciembre 1954, 6292.50083741,      2017, VALLE DEL CAUCA,  6.59527269127,   0.510778519244

Vamos a ver los primeros datos del objeto munic:

head(munic)

Instalamos la libreria \(RColorBrewer\), la cual usaremos para cambiar los colores de los plot:

#install.packages("RColorBrewer")

llamamos la libreria:

library(RColorBrewer)

Seleccionamos solo la Ciudad Capital del Departamento, en este caso Cali:

cali <- munic[munic$MPIO_CNMBR=="CALI",]
plot(cali, main="Cali", axes=TRUE, col=brewer.pal(n=3, name = 'Dark2'))
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

para descagrar los datos de elevacion y almacenarlos en un nuevo objeto que llamaremos elevatio usamos:

elevation <- get_elev_raster(cali, z=8)
## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output

## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output

## Warning in sp::proj4string(locations): CRS object has comment, which is lost in
## output
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Merging DEMs
## Reprojecting DEM to original projection
## Note: Elevation units are in meters.
## Note: The coordinate reference system is:
##  +proj=longlat +datum=WGS84 +no_defs

Ahora, vamos a revisar que hay dentro del objeto elevation:

elevation
## class      : RasterLayer 
## dimensions : 1035, 1033, 1069155  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -77.3575, -74.51675, 1.392743, 4.228643  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -598.1536, 5331.562  (min, max)

Observamos algunas cosas sobre este DEM:

-dimensiones: el “tamaño” del archivo en píxeles

-nrow, ncol: el número de filas y columnas en los datos.

-ncells: el número total de píxeles o celdas que componen el ráster.

-resolution: el tamaño de cada píxel (en grados decimales en este caso). Recordemos que 1 grado decimal representa unos 111,11 km.

-extent: la extensión espacial del ráster. Este valor estará en las mismas unidades de coordenadas que la coordenada sistema de referencia del ráster.

-crs: la cadena del sistema de referencia de coordenadas para el ráster. Este ráster está en coordenadas geográficas con un datum de WGS 84.

library(scales)

Ploteamos los datos de elevacion de Cali:

plot(elevation, main="Descarga del DEM de Cali [metros]")
plot(cali, add= TRUE)

4. Recortar los datos de elevacion para que conincidan con el area de estudio

El DEM cubre una extensión mayor que la que necesitamos. Esto se debe a la disposición espacial de la mosaicos de elevación en AWS.

Usamos la funcion \(crop\) con la cual alineamos Areas incluidas a x e y (Ajustar el area de la ciudad que etamos trbajando), creando una extension que llamaremos elev_crop:

elev_crop = crop(elevation, cali)
plot(elev_crop, main="Modelo de elevacion digital de Cali")
plot(cali, add=TRUE)

Vamos a revisar el nuebo objeto elev_crop:

elev_crop
## class      : RasterLayer 
## dimensions : 101, 90, 9090  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -76.70575, -76.45825, 3.272383, 3.549123  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 944.2498, 3979.236  (min, max)

5. Reproyectar los datos de elevacion

Al trabajar con DEM, siempre es una buena idea utilizar coordenadas de mapa en lugar de coordenadas geográficas. Esto se debe al hecho de que, en las coordenadas geográficas, las unidades de las dimensiones horizontales son grados decimales PERO las unidades de dimension vertical son metros. Reproyectemos los datos de elevación. Podemos ir a epsg.io y buscar la proyección de la zona MAGNA Colombia Cali. Necesitamos obtener la definición de esta referencia espacial en formato PROJ.4 (el que se usa 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"

Ahora, podemos reprojectar los datos de elevacion usando la funcion \(projectExtent\), de coordenadas geograficas WGS84 a MAGNA para la zona de Cali Colombia y los almacenamos en un nuevo objeto Raster que llamaremos pr3:

pr3 <- projectExtent(elev_crop, spatialref)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

Ajustamos el tamaño de celda del objeto pr3 usando la funcion \(res\):

res(pr3) <- 100

Ahora proyectamos los datos de elev_crop usando la funcion \(projectRaster\) y almacenamos los datos en un nuevo objeto que llamaremos rep_elev:

rep_elev <- projectRaster(elev_crop, pr3)

Observamos que hay en el objeto rep_elev:

rep_elev
## class      : RasterLayer 
## dimensions : 307, 276, 84732  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 707796.8, 735396.8, -146070.7, -115370.7  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : 943.4611, 3963.625  (min, max)

Ahora, vamos a reproyectar el SpacialPolygonsDataFrame representando la capital del Valle del Cauca, almacenamos la reproyeccion en un nuevo objeto que llamaremos rep_cali:

(rep_cali = spTransform(cali, spatialref))
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved
## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 707694.8, 735285.7, -145948.9, -115539.4  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=4.59620041666667 +lon_0=-74.0775079166667 +k=1 +x_0=1000000 +y_0=0 +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       :         76,      76001,       CALI,       1536, 563.04276445,      2017, VALLE DEL CAUCA, 1.16366974876, 0.0457332109324

Ploteamos los objetos rep_elev y rep_cali usando la funcion \(plot\):

plot(rep_elev, main="Reproyeccion del Modelo Digital de elevacion de Cali")
plot(rep_cali, add=TRUE)

Guardamos el DEM usando la funcion \(writeRaster\), almacenandolo en nuestro directorio de trabajo:

writeRaster(rep_elev,filename = "C:/Users/FLCS/Documents/R/Geomatica Basica/GB/Notebook/rep_cali_elev.tif", datatype='INT4S', overwrite=TRUE)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
##  but +towgs84= values preserved

6. Estadisticas basicas de datos de elevacion

Una exploracion rapida de las estadisticas del DEM. (Histograma) :

hist(x=rep_elev, main="Histograma de Reproyeccion", col ="blue")

A continuacion calculamos las etadisticas (promedio, maximo, minimo y desviacion) del objeto Raster rep_elev, usando la funcion \(cellStats\):

promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion <- cellStats(rep_elev, 'sd')

Creamos dos vectores unidimensionales que llamaremos metricas y valores donde almacenamos los calculos anteriores:

metricas <- c('mean', 'min', 'max', 'sd')
valores <- c(promedio, minimo, maximo, desviacion)

Creamo un data.frame que llamaremos df_estadisticas, con las estadisticas de elevacion en metros usando la funcion \(data.frame\). (usamos los parentesis exteriores para “ver” el contenido del objeto):

(df_estadisticas <- data.frame(metricas, valores))

7. Obtencion de Variables geomorfometricas

Primero usamos las funciones \(terrain\) (con la que calculamos Pendiente y Aspecto), y la funcion \(hillShade\) (para calcular la sombra de de colinas):

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

Ploteamos los datos de elevacion almacenados en el objeto rep_elev:

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

Ploteamos los datos de la pendiente almacenados en el objeto slope:

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

Ploteamos los datos de Aspecto almacenados en el objeto aspect:

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

Ploteamos una combinacion entre Sombra de colina hillShade y el modelo digital de elevacion DEM:

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

7. Cartografiando datos de elevacion con rayshader

La biblioteca rayshader es un paquete de código abierto para producir visualizaciones de datos 2D y 3D en los usos de R. rayshader datos de elevación en una matriz R base y una combinación de trazado de rayos, mapeo de textura esférica, superposiciones y ambiente oclusión para generar hermosos mapas topográficos 2D y 3D. Además de los mapas, rayshader también permite al usuario traduce objetos ggplot2 en hermosas visualizaciones de datos en 3D.

Instalamos la libreria \(rayshader\):

#install.packages("rayshader")

llamamos la libreria:

library(rayshader)

Usamos la funcion raster_to_matrix para convertir el DEM(rep_elev) en matriz, almacenamos la matriz en un objeto que llamaremos elmat:

elmat = raster_to_matrix(rep_elev)

Usamos otra de las texturas integradas de rayshader para crear:

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

Usamos las funciones \(detect\)\(water\) y \(add\)\(water\) para agregar al mapa una capa de agua:

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

Agregamos una capa con la cual podemos ver la sombra trazada por los rayos del sol usando las funciones add_shadow y ray_shade:

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

8. Otro modo de visualizacion

Instalamos la libreria \(jpeg\):

#install.packages("jpeg")

Llamamos la libreria:

library(jpeg)

En este plot usamo una funcion de sombreado de mapeo de entorno espefico:

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
}

Plotemos otro modo en el cual usamos un jpg que hemos descargado previamente:

map=readJPEG("C:/Users/FLCS/Documents/R/Geomatica Basica/GB/Notebook/9pvbHjN.jpg")
out = getv(map, aspect, slope)
plotRGB(out, main = "Supuestamente montañas hermosas en Cali")