##Angelica Martinez## Datos de elevacion en R —

En este notebook se ilustraran diferentes funciones para obtener, procesar y visualizar modelos de elevacion digital.

1. Instalamoslos paquetes necesarios

Para lograrlos se deben instalar los siguientes paquetes.

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

2. Se adjuntan los paquetes Una vez instalados se deben cargar.

library(raster)

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

## Loading required package: sp

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

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: C:/Users/Usuario/Documents/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: C:/Users/Usuario/Documents/R/win-library/3.6/rgdal/proj
##  Linking to sp version: 1.4-1

library(lattice)
library(latticeExtra)

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

library(sp)
library(raster)

library(rasterVis)

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

library(elevatr)

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

3. Obtener datos de elevación raster

Se trabajara con el archivo del Censo del Dane, que fue descargado con antelación. Se listara el contenido del archivo llamado “ADMINISTRATIVO”

list.files("20_CESAR/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"

Ahora leeremos el archivo con una utilidad del paquete raster.

(munic <-  shapefile("20_CESAR/ADMINISTRATIVO/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 +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  :         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

Visualicemos los atributos de munic

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

Se seleccionara la capital del departamento de Cesar “Valledupar” y realizaremos un ploteo

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

Ahora se trabajara con el paquete elevation

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

Se visualizara el contenido de elevation

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

Es de resaltar que el sistema de coordenadas de refencia es el WGS84

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

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

#writeRaster(elevation, cesar1="./dem/elev_z8.tif", datatype='INT4S', overwrite=TRUE)

Cortaremos los datos correspondientes a Valledupar, y su abreviatura sera " valle"

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

Se chequea el nuevo elemento.

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

5. Reproyectar los datos de elevación

En la pagina epsg se pueden obtener los datos de MAGNA Colombia Bogota zone projection. se necesita la definicion de PROJEC 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" 

Podemos reproyectarlos datos de elevacion de WGS84 geographic coordinates a MAGNA Colombia Bogota zone.

pr3 <- projectExtent(elev_crop, spatialref)
 
res(pr3) <- 100

rep_elev <- projectRaster(elev_crop, pr3)

Lo que obtendriamos sera:

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.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     : 46.24819, 5521.017  (min, max)

Ahora se reproyectara el SpatialPolygonsDataFrame de la capital del Cesar:

(rep_valle = spTransform(valle,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1026104, 1109555, 1577583, 1693669  (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       :         20,      20001, VALLEDUPAR,       1915, 4185.79157206,      2017,      CESAR, 4.87780921906, 0.345529320156

Esta es la imagen que obtenemos al plotear

plot(rep_elev, main="Reprojected Digital elevation model")
plot(rep_valle, add=TRUE)

Con elsiguiente codigo guardemos el cambio.

#writeRaster(rep_elev, cesar1="./dem/rep_valle_elev.tif", datatype='INT4S', overwrite=TRUE)

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

Haremos un histograma de las elevaciones , para revisar algo de estadistica.

hist(rep_elev)

# para archivos grandes sepueden trabajar los siguientes codigos 
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)

Obtendremos estos datos

# creacion de data frame con estadisticas de elevacion [metros]
# el paréntesis externo sirve para "imprimir" el contenido de un objeto
(df_estadisticas <- data.frame(metricas, valores))

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

7. Estadísticas básicas de datos de dimensiones

Antes de introducirse en este tema se deben tener conocimientos previos del tema de pendientes

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 for Valledupar [meters]", col=terrain.colors(25,alpha=0.7))

En nuestros ploters , la pendiente se observa por los cambios de tonos y de colores

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

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

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

8. Mapeo de datos de elevación con rayshader

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

elmat = raster_to_matrix(rep_elev)

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

Detectar agua

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

Podemos agregar una capa de trazado de rayos desde esa dirección 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 visualizacion

#install.packages(jpeg)

library(jpeg)

Función 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
}

Ahora con una imagen jpeg

map=readJPEG("./9pvbHjN.jpg")


out = getv(map,aspect, slope)


plotRGB(out, main = "El Valle")