DATOS DE ELEVACIÓN EN R

1. INTRODUCCIÓN

El siguiente cuaderno tiene como objetivo ilustrar algunas funcionalidad para obtener, procesar y visualizar Modelos Digitales de Elevación en R.

Para comenzar, los datos de elevación son usados en un amplio rango de aplicaciones, incluyendo, por ejemplo, hidrología y modelamiento ecológico. Para obtener los datos necesarios a través de R se requería de muchos paquetes, sin embargo, esto ya no es necesario debido a que existe una variedad de API que brinda acceso a los datos de elevación. El paquete “elevatr” se escribió para estandarizar el acceso a los datos de elevacón desde las API web.

Para acceder a los datos de elevación, el paquete “elevatr” usa Amazon Web Services Terrain Tiles.

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:/Program Files/R/R-4.0.2/library/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:/Program Files/R/R-4.0.2/library/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)

2. CONSEGUIR DATOS DE ELEVACIÓN

Hay muchos recursos para obtener Modelos de Elevación Digital como el Shuttle Radar Topography Mission (STRM), el USGS National Elevation Dataset (NED), Global DEM (GDEM) y otros. Cada uno de esos DEMs tiene sus ventajas y desventajas.

Antes de su cierre en 2018, Mapzen combinó varias de estos recursos para crear un producto de elevación de síntesis que utiliza los mejores datos de elevación disponibles para una región dada a un nivel de zoom determinado. Esta compilación de datos está dispoible a través de Terrain Tiles en Amazon Web Services (AWS).

La entrada para get_elev_raster() es un data.frame con x(longitud) y y (latitud), cualquier objeto sp o raster se convierte en un RasterLayer de mosaicos que se suponerponen al cuadrado delimitador de la entrada. Si se utilizan varios mosaicos, la salida resultantes es una capa Raster combinada.

3. USANGO GET_ELEV_RASTER() PARA ACCEDER A TERRAIN TILES EN AWS

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

No hay diferencia en usar datos de tipo sp o raster. El data frame requiere un prj. Inicialmente, se necesita subir un shapefile que represente los municipios del departamento, en este caso, Cesar.

Ahora se procede a cargar el shapefile usando la función “shapefile” de la librería ráster

(munic <- shapefile("/Users/Marcela/Documents/UNAL/2020-2/GB/Cesar/MGN_MPIO_POLITICO.shp"))

## class       : SpatialPolygonsDataFrame 
## features    : 25 
## extent      : -74.13916, -72.88575, 7.67435, 10.86767  (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  :         20,      20001,  AGUACHICA,                               1915,   73.00612196,      2017,      CESAR, 0.43487341949, 0.00599288307584 
## max values  :         20,      20787, VALLEDUPAR, Ordenanza 8 de Noviembre 3 de 1971, 4185.79157206,      2017,      CESAR, 4.87780921906,   0.345529320156

Para observar los datos que contiene el archivo se emplea la función head

head(munic)

##   DPTO_CCDGO MPIO_CCDGO  MPIO_CNMBR                               MPIO_CRSLC
## 0         20      20011   AGUACHICA                     Ordenanza 40 de 1914
## 1         20      20032      ASTREA     Ordenanza 13 de Noviembre 26 de 1984
## 2         20      20060    BOSCONIA       Ordenanza 1 de Noviembre 6 de 1979
## 3         20      20175 CHIMICHAGUA                   Ordenanza 0054 de 1892
## 4         20      20178 CHIRIGUANÃ\201                     Ordenanza 04 de 1888
## 5         20      20228   CURUMANÃ\215 Ordenanza 36 del 16 de Noviembre de 1965
##   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
## 0   877.6508      2017      CESAR   1.967938 0.07202222
## 1   637.9028      2017      CESAR   1.335874 0.05252647
## 2   586.2015      2017      CESAR   1.274771 0.04833089
## 3  1376.4446      2017      CESAR   2.473946 0.11324235
## 4  1114.4804      2017      CESAR   2.487206 0.09173969
## 5   915.1712      2017      CESAR   1.606253 0.07529250

En el presente cuaderno se trabajará con la capital del municipio, es decir, Valledupar.

Valle <- munic[munic$MPIO_CNMBR=="VALLEDUPAR",]
plot(Valle, main = "Valledupar", axes= TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.7))

A continuación se descargaran los datos de elevación de la zona de interes.

elevation <- get_elev_raster(Valle, 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

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

Para observar los datos que se acaban de descargar

elevation

## class      : RasterLayer 
## dimensions : 1033, 1544, 1594952  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00271  (x, y)
## extent     : -74.545, -70.299, 8.392522, 11.19195  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -801.2662, 5532.962  (min, max)

La salida anterior permite observar que el sistema de coordenadas de los datos es datum WGS 84. Para graficar los datos se usa la función plot

plot(elevation,main= "Mapa de Elevación")
plot(Valle, add=TRUE)

4. AJUSTAR LOS DATOS A EL ÁREA DE ESTUDIO

Los datos graficados anteriormente presentan un problema debido a que comprenden un área mayor a la de interes, esto se debe a la disposición espacial del arreglo de elevación en AWS.

elev_crop = crop(elevation, Valle)
plot(elev_crop, main= "Modelo de Elevación Digital Ajustado")
plot(Valle, add=TRUE)

Visualizar las características del nuevo objeto

elev_crop

## class      : RasterLayer 
## dimensions : 387, 277, 107199  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00271  (x, y)
## extent     : -73.83825, -73.0765, 9.817982, 10.86675  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 46.89514, 5532.962  (min, max)

Debido al ajuste, el nuevo DEM tiene dimensiones más pequeñas que el anterior, menor cantidad de celdas pero conserva el sistema de referencias WGS 84.

5. REPROJECCIÓN DE LOS DATOS DE ELEVACIÓN

Al trabajar con DEMs resulta más conveniente utilizar coordenadas de un mapa que coordenadas geográficas. Lo anterior se debe a que en las coordenadas geográficas, las unidades horizontales son grados decimales, mientras que las verticales están en metros.

Para comenzar a reproyectar se usará MAGNA Colombia, Bogotá, por lo que se necesita definir esta referencia espacial en formato PROJ.4

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 es posible reroyectar los datos de elevación del sistema de coordenadas WGS 84 a MAGNA Colombia, Bogotá.

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)

Para visualizar lo obtenido solo se debe llamar al objeto

rep_elev

## class      : RasterLayer 
## dimensions : 1162, 837, 972594  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1026160, 1109860, 1577475, 1693675  (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     : 46.24819, 5521.017  (min, max)

A continuación se reproyectaran el SpatialPolygonsDataFrame que representa a Valledupar

(rep_Valle = spTransform(Valle, 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      : 1026104, 1109555, 1577583, 1693669  (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       :         20,      20001, VALLEDUPAR,       1915, 4185.79157206,      2017,      CESAR, 4.87780921906, 0.345529320156

Es momento de graficar

plot(rep_elev, main = "Modelo de Elevación Digital Reproyectado")
plot(rep_Valle, add=TRUE)

6. ESTADÍSTICA BÁSICA DE DATOS ELEVACIÓN

Para una rápida exploración de las estadísticas para DEM

hist(rep_elev)

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

Para crear un vector unidimensional

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

Ahora se puede crear un data frame con las estadísticas de elavación, cabe resaltar que las unidades están en metros

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

##   metricas    valores
## 1     mean  945.31222
## 2      min   46.24819
## 3      max 5521.01701
## 4       sd 1182.98377

7. OBTENER VARIABLES GEOFORMOMÉTRICAS

Primero se calcula la pendiente, aspecto y sombreado

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 Valledupar (metros)", col=terrain.colors(25, alpha = 0.8))

También es posible graficar la pendiente

plot(slope, main="Pendiente de Valledupar (grados)", col=topo.colors(25, alpha = 0.8))

Asimismo, es posible graficar el aspecto

plot(aspect, main= "Aspecto para Valledupar (grados)", col=rainbow(25, alpha= 0.8))

Y por último, la sombra

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

8. MAPEANDO DATOS DE ELEVACIÓN CON RAYSHADER

La librería rayshader se utiliza para producir visualizaciones en 2D y 3D en R. Esta librería usa datos de elevación en una matrix R de base y una combinación de trazado de rayos, mapeo de texturas esfércias, superposiciones y oclusión ambiental para generar mapas topográficas en 2D y 3D. Además de los mapas, rayshader permite al usuario traducir objetos ggplot2 en visualizaciones en 3D.

library(rayshader)

Inicialmente se debe convertir el DEM a una matrix

elmat = raster_to_matrix(rep_elev)

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

Detectar agua y agregar una capa de agua al mapa

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

También se puede agregar una capa con trazado de los rayos del sol

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

9. OTRA FORMA DE VISUALIZACIÓN

Se comenzará instalando 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("./9pvbHjN.jpg")
out = getv(map, aspect, slope)
plotRGB(out, main = "Supuestamente montañas bonitas de valledupar")