title: “Angelica Martinez” output: Datos de elevacion en R. 26 de marzo —

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

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

Una vez instalados se deben cargar.

library(raster)
library(rgl)
library(rgdal)
library(lattice)
library(latticeExtra)
library(sp)
library(raster)

library(rasterVis)

library(elevatr)

3. Get Raster Elevation Data

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("ADMINISTRATIVO")

Ahora leeremosel archivo con una utilidad del paquete raster.

(munic <-  shapefile("ADMINISTRATIVO/MGN_MPIO_POLITICO.shp"))

Visualicemos los atributos de munic

head(munic)

Se seleccionara lacapital 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.5))

Ahora se trabajara con el paquete elevation

elevation <- get_elev_raster(valle, z = 8)

#elevation <-  raster("./dem/elev_z8.tif")

Se visualizara el contenido de elevation

elevation

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

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

Ahora se reproyectara el SpatialPolygonsDataFrame de la capital del Cesar:

(rep_valle = spTransform(valle,spatialref))

Esta esla 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.

# histograma
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')

# create unidimensional vector
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))

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)

elmat = raster_to_matrix(rep_elev)

Mapear texturas

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