Datos de elevación en R

1. Introducción:

En este cuaderno se implementan diferentes métodos de obtener, procesar y visualizar modelos de elevación digital (DEMs) en R, en el departamento de Casanare, en el municipio de Yopal.

Los datos de elevación se implementan en muchas áreas, para muchas operaciones, como la visualización espacial. Para acceder a esta información mediante R, es necesario instalar una serie de librerías, mencionadas a continuación:

###
library(rasterVis)

## Loading required package: raster

## Loading required package: sp

## Loading required package: lattice

## Loading required package: latticeExtra

library(raster)
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/yarya/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/yarya/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(elevatr)
library(lattice)
library(latticeExtra)
##

2. Obtener los datos de elevación raster:

Existen diferentes fuentes para obtener los datos de elevación que pueden ser ejecutados en R, como: the Shuttle Radar Topography Mission (SRTM), the USGS National Elevation Dataset (NED), Global DEM (GDEM), entre otros.Pese a su cierre en enero del 2018, Mapzen combinó muchas de estas fuentes para crear un producto de evaliación sintetizado que utiliza los mejoreas datos de elevación disponibles, para una región dada a un nivel de zoom determinado.

La entrada para “get_elev_raster()” es un data.frame con valores x(longitud) y y(latitud).

El data frame requiere un prj. Los archivos ingresados deben presentar formato .shp, los cuales fueron obtenidos de la base de datos del DANE, específicamente para el departamento de Casanare.

###
munic<-  shapefile("C:/Users/yarya/Desktop/Geomatica/Clase5/DatosCasanare2017/85_CASANARE/ADMINISTRATIVO//MGN_MPIO_POLITICO.shp")
head(munic)

class(munic)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

##

En este caso, se trabajará con la capital del departamento de Casanare, es decir, Yopal

###
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.8))

##

Es necesario descargar los datos de elevación, mediante la función get_elev_raster(municipio, z = zoom deseado). Con esto, es posible generar un mapa de los datos de elevación de la zona seleccionada

###
elevation <- get_elev_raster(yopal, z = 9)

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

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

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

elevation

## class      : RasterLayer 
## dimensions : 1545, 1550, 2394750  (nrow, ncol, ncell)
## resolution : 0.00137, 0.00137  (x, y)
## extent     : -73.13185, -71.00835, 3.506186, 5.622836  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -25.86266, 3845.725  (min, max)

plot(elevation, main="Modelo digital de elevación")
plot(yopal, add=TRUE)

##

3. Recorta los datos de elevación para que coincidan con la extensión del área de estudio:

Para realizar el estudio en la zona especificada (Yopal), se debe restringir la zona, esto mediante la selección de los datos de elevación a emplear, utilizando la función “crop”.

###
elev_crop = crop(elevation, yopal)
plot(elev_crop, main="Modelo digital de elevación (DEM) cortado")
plot(yopal, add=TRUE)

elev_crop

## class      : RasterLayer 
## dimensions : 487, 471, 229377  (nrow, ncol, ncell)
## resolution : 0.00137, 0.00137  (x, y)
## extent     : -72.58385, -71.93858, 4.903586, 5.570776  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 144.3594, 2499.776  (min, max)

##

4. Reproyectar los datos de elevación:

Cuando trabajamos con DEMs, siempre es buena idea usar coordenadas de mapa en lugar de coordenadas geográficas, debido a que, en coordenadas geográficas, las unidades de dimensión hprizontal son grados decimales, pero las unidades de dimensión verticales son metros. Esto para trabajar en ambas dimensiones bajo las mismas unidades.

Esta información s epuede obtener de la página web epsg.io, donde se especificael uso de la proyección MAGNA Colombia Bogota, específicamente en formato PROJ.4, debido a que es el formato implementado en las librerias sp y raster.

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

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

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

## class      : RasterLayer 
## dimensions : 740, 718, 531320  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1165512, 1237312, 1034201, 1108201  (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     : 159.1454, 2491.141  (min, max)

(rep_yopal = spTransform(yopal,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      : 1165549, 1237211, 1034289, 1107987  (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       :         85,      85001,      YOPAL, Ordenanza 0038 del 8 de Julio de 1942, 2482.9042923,      2017,   CASANARE, 3.15923255859, 0.202336772314

##

Ahora podemos realizar un mapa con estas unidades

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

##

5. Estadística básica para los datos de elevación

Es posible realizar diferentes análisis estadísticos de los datos, como la freciencia de los valores de elevación, lo que permite analizar como se presentan estos datos a lo largo del territorio

Por otro lado, es posible analizar los valores medios, máximos y mínimos, lo cual puede permitir una interpretación mucho más clara

###
hist(rep_elev)

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 las variables geomorfométricas

Para el mapeo de las variables como aspecto y pendiente, se requiere computar su determinación, mediante las funciones “slope” y “aspect”. Es recomendable implementar diferentes paletas de colores. El mapa de sombras permite hacer constrastes entre la elevación y las variables determinadas.

###
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)
plot(rep_elev,main="DEM para Yopal [metros]", col=terrain.colors(25,alpha=0.7))

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

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

##

De igual forma, es posible realizar un mapa combinado, el cual contenga el mapa de sombras, y variables como la elevación, la pendiente o el aspecto

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

plot(hill,col=grey(1:100/100),legend=FALSE,main="Pendiente para Yopal",axes=FALSE)          
plot(slope, axes=FALSE,col=topo.colors(12, alpha=0.35), add=TRUE)

plot(hill,col=grey(1:100/100),legend=FALSE,main="Pendiente para Yopal",axes=FALSE)          
plot(aspect, axes=FALSE,col=rainbow(12, alpha=0.35), add=TRUE)

##

6. Mapeo de datos de elevación con rayshader

La librería “rayshader” permite visualizar datos en R de forma bi y tridimensional.

Es necesario instalar la librería “rayshader”, y posteriormente convertir el DEM en una matríz

###
library(rayshader)
elmat = raster_to_matrix(rep_elev) #DEM a matríz
elmat %>% 
  sphere_shade(texture = "imhof2") %>%
plot_map()

##

Es posible añadir capas de agua a la imagen con las funciones detect_water y add_water.

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

##

También es posible añadir sombras para esa dirección del sol, con la fucnión add_shadow, la cual, en el presenta caso, permite dar una idea básica de la diferencia de altitud entre las diferentes zonas de yopal, mediante la proyección de la sombra.

###
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 visuialización

Grant Williamson, experto en R, sugiere esta metodología para visualizar los datos. Para esto, es necesario instalar la librería “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/yarya/Desktop/Geomatica/Clase8/ambiente_azul.jpg")
out = getv(map, aspect, slope)
plotRGB(out, main = "Supposedly pretty mountains in Yopal")

##