Datos de elevacion en R

1. Introducción

Mediante este trabajo se quiere ilustrar varias funciones en donde se puede visualizar los modelos digitales de elevación (DEM) en R. Para esto se trabajara con el Paquete elevatr el cual estandarizo el accesso a los datos de elevación desde las API web. En este caso para dar ejemplo de estas herramientas el trabajo tendra como foco el departamento de Nariño.

El primer paso es instalar los paquetes necesarios:

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

Lo siguiente es llamar a las librerias:

library(sp)
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/USER/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/USER/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(rasterVis)

## Loading required package: lattice

## Loading required package: latticeExtra

library(elevatr)

2. Obtencion daros de elevación

Se debe tener en cuenta que hay varias fuentes de donde obtener los modelos de elevación digital. Con la entrada para getelevraster() para acceder a terrain tiles en AWS el cual es un data.frame con ubicaciones x para longitus e Y para latitud, un objeto sp o un objeto raster necesita ser la entrada y src debe establecerse en “mapzen” (este es un valor predeterminado).

El primer paso es cargar el archivo .shp del departamento de Nariño el cual se descargo del geoportal del DANE. Para esto se usa la funcion *list.files(), podemos observar los archivos que hay dentro del directorio de trabajo:

list.files("C:/Users/USER/Desktop/Geomatica/52_NARIÑO/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"

(mnici <- shapefile("/Users/USER/Desktop/Geomatica/52_NARIÑO/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

## class       : SpatialPolygonsDataFrame 
## features    : 64 
## extent      : -79.01021, -76.83368, 0.3613481, 2.683898  (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  :         52,      52001, ALBÁN (San José),                      1540,   25.31281457,      2017,    NARIÑO, 0.240922725368, 0.00205017524724 
## max values  :         52,      52885,         YACUANQUER, Ordenanzas 7 y 11 de 1871, 3611.35178489,      2017,    NARIÑO,  6.45337364952,   0.291654453146

Lo siguiente es seleccionar la capital del departamento de Nariño, el cual es San Juan de Pasto:

(pasto<-mnici[mnici$MPIO_CNMBR=="PASTO",])

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : -77.38396, -77.0266, 0.8215102, 1.353557  (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 
## value       :         52,      52001,      PASTO,       1540, 1099.41143169,      2017,    NARIÑO, 1.86958100889, 0.089064952915

plot(pasto, main="Pasto", axes=TRUE)
plot(mnici, add=TRUE)
invisible(text(coordinates(mnici), labels= as.character(mnici$MPIO_CNMBR), cex=0.8))

Lo siguiente es descargar los datos de elevación:

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

Lo siguiente es ver lo que esta dentro del objeto elevation y es importante entender los componentes que se muestran:

elevation

## class      : RasterLayer 
## dimensions : 1033, 1544, 1594952  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00275  (x, y)
## extent     : -78.76375, -74.51775, -1.420891, 1.419859  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -402.9963, 5807.503  (min, max)

library(scales)

plot(elevation, main="Mapa elevación Pasto")
plot(pasto, add=TRUE)

4. Enfocar datos de acuerdo al area de estudio

Como se observo, el DEM esta cubriendo un area mayor a la extension que queremos trabajar. Esto se debe a la disposición espacial de los mosaicos de elevación en AWS. Es por ello que se busca filtrar los datos para obtener la grafica de lo que queremos:

elev_crop = crop(elevation, pasto)
plot(elev_crop, main= "DEM pasto ajustado")
plot(pasto, add=TRUE)

Y se comprueban las caracteristicas del objeto:

elev_crop

## class      : RasterLayer 
## dimensions : 194, 130, 25220  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00275  (x, y)
## extent     : -77.38325, -77.02575, 0.8203588, 1.353859  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 1169.861, 4147.482  (min, max)

5. Reproyeccion de 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 coordenadas geográficas, las unidades de dimensiones horizontales son grados decimales pero la unidad de dimensión vertical son 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"

Con esto se puede proyectar los datos con coordenadas geográficas WGS84 en la zona 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)

se observa lo que tiene el objeto:

rep_elev

## class      : RasterLayer 
## dimensions : 591, 399, 235809  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 631838.3, 671738.3, 582620.9, 641720.9  (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     : 1161.735, 4142.733  (min, max)

Lo siguiente es proyectar el SpatialPolygonsDataFrame para Pasto, Nariño:

(rep_pasto = spTransform(pasto, 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      : 631818.2, 671614, 582736.2, 641676.8  (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       :         52,      52001,      PASTO,       1540, 1099.41143169,      2017,    NARIÑO, 1.86958100889, 0.089064952915

Y se realiza el plot:

plot(rep_elev, main = "DEM Pasto Reproyectado")
plot(rep_pasto, add=TRUE)

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

Se realiza un histograma buscando hacer una exploracion rápida de las estadisticas del DEM

hist(rep_elev)

Es posible calcular las estadisticas del objeto Raster que se esta trabajando:

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

Se crean dos vectores unidimensionales donde se almacenaran los datos calculados:

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

Por ultimo se crea el data.frame con las estadisticas de elevacion:

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

7. Obtencion de variables geomorfométricas

Lo primero es calcular la pendiente, la orientación 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)

Se realiza el plot de la elevacion:

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

Ahora se realiza el plot de la pendiente:

plot(slope, main="Pendiente de Pasto [Grados]", col=topo.colors(25, alpha = 0.7))
plot(rep_pasto, add=TRUE)

Tambien se realiza el Plot de la orientación:

plot(aspect, main="Orientación para Pasto [Grados]", col=rainbow(25, alpha = 0.7))
plot(rep_pasto, add=TRUE)

Por ultimo se plotea las sombras:

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

8. Mapeo de datos de elevacion con Rayshader

Mediante el uso de la biblioteca rayshader se puede producir bisualizaciones de datos 2D y 3D en R. Con esto es posiblre realizar modelos que presentan un mayor nivel de detalle. Para el uso de esta herramienta se debe instalar el paquete y realizar el siguiente procedimietno:

#install.packages("rayshader")

library(rayshader)

El primer paso en convertir el DEM a una matriz:

elmat = raster_to_matrix(rep_elev)

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

Para representar el agua en el mapa se puede trabajar con detect_water y add_water:

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

Y tambien es posible agragar una capa del trazado de los rayos solares:

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

9. Otra manera de Visualizar

Aun asi si se quiere hacer una visualizacion distinta, es posible usar el paquete jpeg, para ello como siempre es necesario instalar y llamar la libreria:

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

Y se usa la funcion de sombreado de mapeo de entorno esférico:

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/USER/Desktop/Geomatica/9pvbHjN.jpg")
out = getv(map, aspect, slope)
plotRGB(out, main = "Montañas de Pasto")