1. Introducción.

El presente R_notenook, tiene el objetivo de prersentar algunas funcionalidades para obtener, pocesar y visualizar datos de elevación digital DEM del municipio de Villavicencio en el departamento de Meta en Colombia. Para ello, es necesario realizar la instalación de algunos cuadernos de R, que más adelante en este cuaderno, permitirán realizar el análisis de los datos DEM:

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

Llamando a las librerías:

library(rasterVis)
## Loading required package: raster
## Loading required package: sp
## Loading required package: lattice
## Loading required package: latticeExtra
library(raster)
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/EQUIPO/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/EQUIPO/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(rgl)
library(elevatr)

2. Obtención de datos de elevación-Ráster.

Es necesario cargar el shapefile de los municipio del departamento, usando una funcion del paquete ratser, así:

(munic <-  shapefile("C:/Users/EQUIPO/Documents/2020-2/Geomatica/Trabajo_1/50_META/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))
## class       : SpatialPolygonsDataFrame 
## features    : 29 
## extent      : -74.89921, -71.07753, 1.604238, 4.899101  (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  :         50,      50001,     ACACÍAS,                                 1840,  117.34108343,      2017,       META, 0.43640580036, 0.00955216818621 
## max values  :         50,      50711, VISTAHERMOSA, Ordenanza 63 de Noviembre 21 de 1965, 17247.6864753,      2017,       META, 8.63318071739,    1.40216640329

Revisando los atributos del objeto “munic”:

head(munic)
##   DPTO_CCDGO MPIO_CCDGO        MPIO_CNMBR                           MPIO_CRSLC
## 0         50      50001     VILLAVICENCIO                                 1840
## 1         50      50006          ACACÃ\215AS                 Ordenanza 23 de 1960
## 2         50      50110  BARRANCA DE UPIA   Ordenanza 5  de Octubre 16 de 1990
## 3         50      50124          CABUYARO   Decreto 47 de Septiembre 3 de 1912
## 4         50      50150 CASTILLA LA NUEVA      Ordenanza 08 de Julio 7 de 1961
## 5         50      50223          CUBARRAL Ordenanza 23 de Noviembre 28 de 1960
##   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
## 0  1285.9309      2017       META   2.039774 0.10471811
## 1  1123.7379      2017       META   2.074207 0.09150745
## 2   403.8211      2017       META   1.312018 0.03289441
## 3   913.6705      2017       META   2.008848 0.07440326
## 4   514.6870      2017       META   1.349329 0.04189962
## 5  1157.8655      2017       META   2.075091 0.09426999

Ahora, filtarndo los datos para seleccionar únicamente los referentes al municipio de Villavicencio y posteriormente graficando:

villavicencio<- (munic[munic$MPIO_CNMBR=="VILLAVICENCIO",])
plot(villavicencio, main="Villavicencio", axes=TRUE)
plot(munic, add=TRUE)
 invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

2.1. Usando get_elev_raster() para obtener los DEM de Villavicencio en AWS.

A través de Amazon Web Service, se puede obtener los datos de elevación de Villacicencio. por lo tanto, se usa la función get_elev_raster(), así:

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

Revisando los datos descargados recientemente guardados en elevation:

elevation
## class      : RasterLayer 
## dimensions : 1546, 1033, 1597018  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.545, -71.70425, 1.393646, 5.629686  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : -1.570384, 4150.449  (min, max)

Graficando el objeto elevation:

plot(elevation, main="Mapa de elevación]")
plot(villavicencio, add=TRUE)

En la imágen se observa que el area subrayada es la represenatción de los datos de elevación de Villavicencio, sin embargo, para ser más precisos, en el siguiente numeral se filtarn los datos para adecuarlos al area de estudio. ## 3. Filtrando los datos para adecuarlos al area de estudio. A contnuación se filtran los daos para que al graficar, únicamente se muestren los dats referentes a Villavicencio:

elev_crop = crop(elevation, villavicencio)
plot(elev_crop, main= "DEM ajutsado")
plot(villavicencio, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.8))

Revisando las características del objeto “elev_crop”

  elev_crop
## class      : RasterLayer 
## dimensions : 129, 216, 27864  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -73.7695, -73.1755, 3.936366, 4.289826  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : layer 
## values     : 203.7742, 3577.458  (min, max)

Cuando se trabaja con datos DEM, es mejor idea usar map coordinates en vez de geographic coordinates. Eso es debido a que, en geographic coordinates, las unidades horizontales son grados decimales, pero, las dimensiones verticales estan en metros. Ahora bien, para comenzar a trabajar, se usará la proyección MAGNA Colombia, Bogotá, para ello, es necesario obtener la definición para ese marco de referencia espacial en PROJ.4, a trve’s del siguiente comando:

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, se pueden reproyectar los datos de elevación del sistema WGS 84 en MAGNA Colombia Bogota zone.

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)

Llamando al objeto rep_elev:

rep_elev
## class      : RasterLayer 
## dimensions : 391, 660, 258060  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1034192, 1100192, 927079.8, 966179.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 
## source     : memory
## names      : layer 
## values     : 201.3911, 3585.355  (min, max)

Ahora, proyectando el SpatialPolygonsDataFrame que representa a Villavicencio, Meta:

(rep_villavicencio = spTransform(villavicencio,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      : 1034342, 1100224, 926940.8, 965985.1  (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       :         50,      50001, VILLAVICENCIO,       1840, 1285.93087074,      2017,       META, 2.03977418027, 0.104718109413

Graficando:

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

4. Estadística de los datos de elevación.

Realizando un histograma del objeto *rep_elev

hist(rep_elev)

Realizando calculos para entender los datos:

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

Creando vector unidimendional

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

Creando data frame con estadísticas de elevación:

(df_estadisticas <- data.frame(metricas, valores))
##   metricas   valores
## 1     mean  449.7664
## 2      min  201.3911
## 3      max 3585.3554
## 4      std  394.3298

5. Obtención de variables geomorfométricas.

Inicialmente, es necesario obtener la pendiente, la orientación y el sombrado de la colina:

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

Graficando:

plot(rep_elev, main= "DEM para Villavicencio [metros]", col=terrain.colors(36, alpha = 1))
plot(rep_villavicencio, add=TRUE)

Tambien, garficando la pendiente:

 plot(slope,main="Slope for Villavicencio [degrees]", col=topo.colors(20,alpha=0.9))
plot(rep_villavicencio, add=TRUE)

De igual forma, grafiando la orientación:

plot(aspect,main="Aspect for Villavicencio [degrees]", col=rainbow(25,alpha=0.9))
plot(rep_villavicencio, add=TRUE)

La misma gráfica pero con la sombra generada por las colinas:

plot(hill,
        col=grey(1:100/100),  
        legend=FALSE,        
        main="DEM para Villavicencio",
        axes=FALSE)           

plot(rep_elev, 
        axes=FALSE,
        col=terrain.colors(12, alpha=0.35), add=TRUE)
plot(rep_villavicencio, add=TRUE)

6. Graficando datos de elevación con Rayshader.

Rayshader es una herramienta que permite realizar modelos con mayor nivel de detalle, acontinuación se presenta una serie de ejemplos al respecto, para ello, es necesario instalar y llamar a la librería:

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

Convirtendo DEM en una matriz:

elmat <- raster_to_matrix(rep_elev)

Grfaicando con Rayshader

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

Usando la funcion detect_water y add_water, para representar el agua en el mapa:

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

De igual forma, se puede agregar una capa que muestre los rayos del solsobre la superficie:

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

7. Otra forma de visualización.

Para concluir con este cuaderno, se van a mostrar las bonadades del paquete "jpeg, para ello, es necsario instalar y llamar a la misma librería:

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

Usando una función de mapeo de la librería:

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
}

Graficando:

map=readJPEG("C:/Users/EQUIPO/Documents/2020-2/Geomatica/mapas/9pvbHjN.jpg")
out = getv(map, aspect, slope)
plotRGB(out)

Se observa que el municipio de Vilavicencio, presenta únicamente alguna colinas hacia el nor-occidente, lo que, teniendo en cuenta la geografía colombiana, serían las últimas estribaciones de la cordillera orinental antes de salir a los llanos orientales.