Datos de elevación en R

Introducción a elevatr

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

library(rasterVis)

## Loading required package: raster

## Loading required package: sp

## Loading required package: lattice

## Loading required package: latticeExtra

library(raster)
library(rgl)
library(rgdal)

## 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(elevatr)

Obtener datos de elevación de ráster

list.files("C:/Users/USUARIO/Documents/Geomatica/Tolima")

##  [1] "9pvbHjN.jpg"    "dem.tif"        "depto.cpg"      "depto.dbf"     
##  [5] "depto.prj"      "depto.qpj"      "depto.shp"      "depto.shx"     
##  [9] "elevacion"      "elevacion.grd"  "elevacion.gri"  "municipios.cpg"
## [13] "municipios.dbf" "municipios.prj" "municipios.qpj" "municipios.shp"
## [17] "municipios.shx" "vias.dbf"       "vias.prj"       "vias.qpj"      
## [21] "vias.shp"       "vias.shx"

Ahora, leamos shapefile usando una función provista por el paquete ráster:

(munic <- shapefile("C:/Users/USUARIO/Documents/Geomatica/Tolima/municipios.shp"))

## class       : SpatialPolygonsDataFrame 
## features    : 46 
## extent      : -76.1007, -74.5033, 2.9669, 5.2915  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 11
## names       : ID_0, ISO,   NAME_0, ID_1, NAME_1, ID_2,     NAME_2,    TYPE_2,    ENGTYPE_2, NL_NAME_2,     VARNAME_2 
## min values  :   53, COL, Colombia,   29, Tolima, 1000,  Alpujarra, Municipio, Municipality,        NA,     Mariquita 
## max values  :   53, COL, Colombia,   29, Tolima,  999, Villarrica, Municipio, Municipality,        NA, Villa Hermosa

¿Cuáles son los atributos del objeto munic?

head(munic)

##   ID_0 ISO   NAME_0 ID_1 NAME_1 ID_2      NAME_2    TYPE_2    ENGTYPE_2
## 0   53 COL Colombia   29 Tolima  966   Alpujarra Municipio Municipality
## 1   53 COL Colombia   29 Tolima  967    Alvarado Municipio Municipality
## 2   53 COL Colombia   29 Tolima  968    Ambalema Municipio Municipality
## 3   53 COL Colombia   29 Tolima  969 Anzoátegui Municipio Municipality
## 4   53 COL Colombia   29 Tolima  970      Armero Municipio Municipality
## 5   53 COL Colombia   29 Tolima  971       Ataco Municipio Municipality
##   NL_NAME_2 VARNAME_2
## 0      <NA>      <NA>
## 1      <NA>      <NA>
## 2      <NA>      <NA>
## 3      <NA>      <NA>
## 4      <NA>      <NA>
## 5      <NA>      <NA>

Seleccionemos solo la ciudad capital de este departamento.

lechona <- munic[munic$NAME_2=="Ibagué",]
plot(lechona, main="Ibagué", axes=TRUE)
plot(munic, add=TRUE)
invisible(text(coordinates(munic), labels =        as.character(munic$NAME_2), cex=0.8))

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

elevation

## class      : RasterLayer 
## dimensions : 1034, 1033, 1068122  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -75.95125, -73.1105, 2.796526, 5.629686  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 139.1862, 5253.441  (min, max)

plot(elevation, main="Este es el DEM descargado (m)")
plot(lechona, add=TRUE)

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

elev_crop = crop(elevation, lechona)
plot(elev_crop, main="Modelo de elevación digital recortado")
plot(lechona, add=TRUE)

Veamos el nuevo objeto:

elev_crop

## class      : RasterLayer 
## dimensions : 140, 189, 26460  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -75.54425, -75.0245, 4.262426, 4.646026  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 453.4822, 4417.787  (min, max)

reproyectacion de los datos de elevacion en un nuevo sistema de coordenadas (MAGNA-Colombia)

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 "

Ajuste de las caracterisicas del raster

pr3 <- projectExtent(elev_crop, spatialref)
# Adjust the cell size 
res(pr3) <- 100
# now project
rep_elev <- projectRaster(elev_crop, pr3)

Ahora observamos los atributos del raster recien creado

rep_elev

## class      : RasterLayer 
## dimensions : 425, 578, 245650  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 837153.9, 894953.9, 963178.5, 1005679  (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     : 442.3035, 4498.308  (min, max)

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

(rep_lechona = spTransform(lechona, spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 837006.5, 894877, 963387.1, 1005692  (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   : 11
## names       : ID_0, ISO,   NAME_0, ID_1, NAME_1, ID_2,  NAME_2,    TYPE_2,    ENGTYPE_2, NL_NAME_2, VARNAME_2 
## value       :   53, COL, Colombia,   29, Tolima,  987, Ibagué, Municipio, Municipality,        NA,        NA

plot(rep_elev, main="Modelo de elevación digital reproyectado")
plot(rep_lechona, add=TRUE)

A continuación guardaremos nuestro DEM:

writeRaster(rep_elev, filename = "C:/Users/USUARIO/Documents/Geomatica/Tolima/elevacion", dataType='INT4S', overwrite=TRUE)

## Warning in .datatype(...): argument "datatype" misspelled as "dataType"

## class      : RasterLayer 
## dimensions : 425, 578, 245650  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 837153.9, 894953.9, 963178.5, 1005679  (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     : C:/Users/USUARIO/Documents/Geomatica/Tolima/elevacion.grd 
## names      : layer 
## values     : 442, 4498  (min, max)

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

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', 'sd')
valores = c(promedio, minimo, maximo, desviacion)

(df_estadisticas = data.frame(metricas, valores))

##   metricas   valores
## 1     mean 2090.3218
## 2      min  442.3035
## 3      max 4498.3082
## 4       sd  906.5085

Obtención de variables geomorfométricas.

Los siguientes objetos seran necesarios para las variables geomorfometricas, que son 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)

Plot de elevación el cual nos muestra los datos de elevacion bajo el sistema de coordenadas MAGNA

plot(rep_elev, main="DEM para Ibagué (m)", col=terrain.colors(20, alpha = 1))

plot(slope, main="Pendiente para Ibagué (grados)", col=topo.colors(25, alpha = 1))

plot(aspect, main="Aspecto para Ibagué (grados)", col=rainbow(25, alpha = 1))

A continuación es un Plot combinado:

plot(hill, col=grey(1:100/100),legend=FALSE, main="DEM para Ibagué", 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 Primero debemos instalar rayshader con el siguiente comando:

#install.packages("rayshader")

library(rayshader)

elmat = raster_to_matrix(rep_elev)

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

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

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

### Otra forma de visualizarlo:

#install.packages("jpeg")

library(jpeg)

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/USUARIO/Documents/Geomatica/Tolima/9pvbHjN.jpg")
out = getv(map, aspect, slope)

plotRGB(out, main = "Montañas supuestamente bonitas en Ibagué")

Gracias por la atención prestada.