###01.04.2020 ##Andrea Valentina Galindo Rojas

#Introducción a elevatr ###Los datos de elevación se utilizan para una amplia gama de aplicaciones, entre ellas, por ejemplo, la visualización, la hidrología y la modelización ecológica. El acceso a estos datos en R no ha tenido una sola interfaz, está disponible a través de funciones en muchos paquetes, o requiere un 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 fue escrito para estandarizar el acceso a los datos de elevación de las APIs de la web. Para acceder a los datos de elevación rasterizados (por ejemplo, un DEM) el paquete elevatr utiliza los Azulejos de Terreno de Amazon Web Services.

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

## Warning: package 'rasterVis' was built under R version 3.6.3

## Loading required package: raster

## Warning: package 'raster' was built under R version 3.6.3

## Loading required package: sp

## Loading required package: lattice

## Loading required package: latticeExtra

## Warning: package 'latticeExtra' was built under R version 3.6.3

library(raster)
library(rgl)

## Warning: package 'rgl' was built under R version 3.6.3

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.6.3

## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: D:/Mis documentos/R/win-library/3.6/rgdal/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: D:/Mis documentos/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1

library(elevatr)

## Warning: package 'elevatr' was built under R version 3.6.3

Obtención de datos de elevación Raster

###Hay varias fuentes para los modelos de elevación digital como la Shuttle Radar Topography Mission (SRTM), el USGS National Elevation Dataset (NED), Global DEM (GDEM), y otras. Cada uno de estos DEM tiene pros y contras para su uso. Antes de su cierre en enero de 2018, Mapzen combinó varias de estas fuentes para crear un producto de elevación de síntesis que utiliza los mejores datos de elevación disponibles para una región determinada a un nivel de zoom determinado. Aunque cerrados, estos datos compilados por Mapzen son actualmente accesibles a través de los Terrain Tiles on Amazon Web Services (AWS). La entrada para get_elev_raster() es un data. frame con las ubicaciones x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto sp, o cualquier objeto raster y devuelve un RasterLayer de los azulejos que se superponen al cuadro delimitador de la entrada. Si se recuperan varios azulejos, la salida resultante es una Capa de Raster fusionada

##Uso de get_elev_raster() para acceder a los azulejos del terreno en AWS.

En este cuaderno, cargaremos el archivo de forma descargado del geoportal del DANE. Revisemos el contenido de la carpeta:

list.files("./RStudio_Works/52_NARIÑO/ADMINISTRATIVO")

## character(0)

###Ahora, leamos el archivo de forma usando una función provista por el paquete ráster:

(munic<-shapefile("D:/Mis documentos/RStudio_Works/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 +ellps=WGS84 +towgs84=0,0,0 
## 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

###¿Cuáles son los atributos del objeto municipal?

head(munic)

##   DPTO_CCDGO MPIO_CCDGO   MPIO_CNMBR                         MPIO_CRSLC
## 0         52      52683     SANDONÃ\201               Ordenanza 33 de 1868
## 1         52      52685 SAN BERNARDO Ordenanza 23 Noviembre 26  de 1992
## 2         52      52687  SAN LORENZO                               1886
## 3         52      52560      POTOSÃ\215                               1571
## 4         52      52693    SAN PABLO                               1773
## 5         52      52565  PROVIDENCIA  Ordenanza 34 Noviembre 27 de 1992
##   MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng  Shape_Area
## 0  101.61251      2017    NARIÃ‘O  0.4851210 0.008228253
## 1   65.16451      2017    NARIÃ‘O  0.4702503 0.005281813
## 2  249.14359      2017    NARIÃ‘O  0.7816641 0.020186623
## 3  388.29816      2017    NARIÃ‘O  1.2816483 0.031439880
## 4  111.85964      2017    NARIÃ‘O  0.5022925 0.009068012
## 5   39.94762      2017    NARIÃ‘O  0.3518080 0.003233792

###Seleccionemos solo la ciudad capital de este departamento.

pasto <- munic[munic$MPIO_CNMBR=="PASTO",]
plot(pasto, main="PASTO", axes=TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels=as.character(munic$MPIO_CNMBR), cex=0.5))

elevation <- get_elev_raster(pasto, z = 8)

## 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 +ellps=WGS84 +towgs84=0,0,0

###Ahora, verifiquemos qué hay dentro del objeto de elevación:

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 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : -402.9963, 5807.503  (min, max)

##Observe 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 (imagine una hoja de cálculo o una matriz). ###Ncells: el número total de píxeles o celdas que componen el ráster. ###Resolución: el tamaño de cada píxel (en grados decinales en este caso). Recordemos que 1 grado decimal representa aproximadamente 111,11 km. ###Extensión: la extensión espacial de la trama. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas 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 dato de WGS 84.

plot(elevation, main="This the downloaded DEM [meters]")
plot(pasto, add=TRUE)

Recorte los datos de elevación para que coincidan con la extensión del área de estudio

###Tenga en cuenta que el DEM cubre una extensión mayor que la que necesitamos. Esto se debe a la disposición espacial de los mosaicos de elevación en AWS.

###De todos modos, es una buena idea guardar el DEM para el futuro.

writeRaster(elevation, filename="D:/Mis documentos/RStudio_Works/52_NARIÑO/elev_z8.tif", datatype='INT4S', overwrite=TRUE)

###Ahora, recortemos los datos de elevación correspondientes a Pasto:

elev_crop = crop(elevation, pasto)
plot(elev_crop, main="Cropped Digital elevation model")
plot(pasto, add=TRUE)

###Veamos el nuevo 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 +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 1169.861, 4147.482  (min, max)

#Reproyectar los datos de elevación ###Cuando se trabaja con DEM, siempre es una buena idea usar 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 es metros. Volvemos a proyectar los datos de elevación.

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 reproyectar los datos de elevación de las coordenadas geográficas WGS84 en la zona MAGNA Colombia Bogotá.

pr3 <- projectExtent(elev_crop, spatialref)

res(pr3) <- 100

rep_elev <- projectRaster(elev_crop, pr3)

¿Qué obtenemos?

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.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 
## source     : memory
## names      : layer 
## values     : 1161.735, 4142.733  (min, max)

###Ahora, reproyectemos el SpatialPolygonsDataFrame que representa la capital de nuestro departamento:

(rep_pasto = spTransform(pasto,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 631818.2, 671614, 582736.2, 641676.8  (xmin, xmax, ymin, ymax)
## crs         : +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 
## 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(rep_elev, main="Reprojected Digital elevation model")
plot(rep_pasto, add=TRUE)

Guardamos el DEM:

writeRaster(rep_elev, filename="D:/Mis documentos/RStudio_Works/52_NARIÑO/rep_pasto_elev.tif", datatype='INT4S', overwrite=TRUE)

#Estadísticas básicas de datos de elevación ###Una exploración rápida de las estadísticas DEM:

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

##   metricas  valores
## 1     mean 2937.657
## 2      min 1161.735
## 3      max 4142.733
## 4      std  414.363

#Obtención de variables geomorfométricas.

###Primero, calcule la pendiente, el aspecto 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)

Parcela de elevación. Tenga en cuenta la paleta de colores utilizada aquí.

plot(rep_elev,main="DEM for Pasto [meters]", col=terrain.colors(25,alpha=0.7))

###Parcela pendiente. Tenga en cuenta la paleta de colores utilizada aquí.

plot(slope,main="Slope for Pasto [degrees]", col=topo.colors(25,alpha=0.7))

###Aspecto de la trama. Tenga en cuenta la paleta de colores utilizada aquí.

plot(aspect,main="Aspect for Pasto [degrees]", col=rainbow(25,alpha=0.7))

###Una trama combinada:

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

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

#Mapeo de datos de elevación con sombreador de rayos ###La biblioteca rayshader es un paquete de código abierto para producir visualizaciones de datos 2D y 3D en R. Rayshader utiliza datos de elevación en una matriz base R y una combinación de trazado de rayos, mapeo de texturas esféricas, superposiciones y oclusión ambiental para generar hermosos mapas topográficos 2D y 3D . Además de los mapas, rayshader también permite al usuario traducir objetos ggplot2 en hermosas visualizaciones de datos 3D.

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

## Warning: package 'rayshader' was built under R version 3.6.3

#Convert the DEM into a matrix:
elmat = raster_to_matrix(rep_elev)

#We use another one of rayshader's built-in textures:
elmat %>%
  sphere_shade(texture = "imhof2") %>%
  plot_map()

#detect_water and add_water adds a water layer to the map:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  plot_map()

#And we can add a raytraced layer from that sun direction as well:
elmat %>%
  sphere_shade(texture = "desert") %>%
  add_water(detect_water(elmat), color = "desert") %>%
  add_shadow(ray_shade(elmat), 0.5) %>%
  plot_map()

#Otra forma de visualización. ##Un experto en R sugiere otra forma de visualización aquí.

###Vamos a intentarlo:

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

# Spherical environment mapping hillshade function
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("D:/Mis documentos/RStudio_Works/pasto.jpg")
# Do the mapping
out = getv(map, aspect, slope)
# Plot pretty mountains
plotRGB(out,main="Montañas de Pasto")

DATOS DE ELEVACIÓN EN R

Obtención de datos de elevación Raster

En este cuaderno, cargaremos el archivo de forma descargado del geoportal del DANE. Revisemos el contenido de la carpeta:

Recorte los datos de elevación para que coincidan con la extensión del área de estudio