Datos de elevación en R

Introducción a elevatr

Los datos de elevación se utilizan para una amplia gama de aplicaciones, que incluyen, por ejemplo, visualización, hidrología y modelado ecológico. Obtener acceso a estos datos en R no ha tenido una sola interfaz, está disponible a través de funciones en muchos paquetes o requiere 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 se escribió para estandarizar el acceso a los datos de elevación desde las API web.

Para acceder a los datos de elevación ráster (por ejemplo, un DEM), el paquete elevatr utiliza los mosaicos de terreno de Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.

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

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/Brian/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/Brian/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

Existen varias fuentes para modelos digitales de elevación como la Misión de Topografía por Radar Shuttle (SRTM), el Conjunto de datos de elevación nacional (NED) de USGS, el DEM global (GDEM) y otros. Cada uno de estos DEM tiene ventajas y desventajas 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 dada a un nivel de zoom dado. Aunque cerrado, estos datos compilados por Mapzen son actualmente accesibles a través de Terrain Tiles en Amazon Web Services (AWS).

La entrada para get_elev_raster () es un data.frame con ubicaciones x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto sp o cualquier objeto ráster y devuelve un RasterLayer de los mosaicos que se superponen al cuadro delimitador de la entrada Si se recuperan múltiples mosaicos, la salida resultante es una capa ráster fusionada.

**Usando get_elev_raster () para acceder a los Terrain Tiles en AWS.**

Como se mencionó, un dataframe con columnas x e y, un objeto sp o un objeto ráster debe ser la entrada y el src debe establecerse en “mapzen” (este es el valor predeterminado). No hay diferencia en el uso de los tipos de datos de entrada sp y ráster. El dataframe requiere un prj. El nivel de zoom (z) predeterminado es 9 (una compensación entre la resolución y el tiempo de descarga), pero no se pudo obtener nada con este zoom. Se tiene que usar un zoom inferior igual a 8.

Primero, necesitamos cargar un archivo shape que represente a los municipios de nuestro departamento. En este cuaderno, cargaremos el shapefile descargado del geoportal del DANE según lo solicitado en la clase virtual el 19 de marzo de 2020. Revisemos el contenido de la carpeta:

list.files("C:/Users/Brian/Desktop/Daniela/Boyacá")

##  [1] "19_03_2020.qgz"     "Boyacá.cpg"         "Boyacá.dbf"        
##  [4] "Boyacá.prj"         "boyacá.qgz"         "Boyacá.shp"        
##  [7] "Boyacá.shx"         "Departamento.gpkg"  "imagen.jpg"        
## [10] "mun_3116.gpkg"      "municipios.cpg"     "municipios.dbf"    
## [13] "municipios.prj"     "municipios.shp"     "municipios.shx"    
## [16] "n.jpg"              "rep_Ruana_elev.grd" "rep_Ruana_elev.gri"
## [19] "rep_Ruana_elev.tif" "Thumbs.db"

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

(munic = shapefile("C:/Users/Brian/Desktop/Daniela/Boyacá/municipios"))

## class       : SpatialPolygonsDataFrame 
## features    : 123 
## extent      : -74.6743, -71.98309, 4.642001, 7.0275  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 5
## names       :  NAME_1, ID_2,     NAME_2,                     TYPE_2,                  ENGTYPE_2 
## min values  : BoyacÃ¡,  199,    Almeida, Corregimiento Departamento, Corregimiento Departamento 
## max values  : BoyacÃ¡,  321, ZetaquirÃ¡,                  Municipio,               Municipality

¿Cuáles son los atributos del objeto munic?

head(munic)

##    NAME_1 ID_2    NAME_2    TYPE_2    ENGTYPE_2
## 0 BoyacÃ¡  199   Almeida Municipio Municipality
## 1 BoyacÃ¡  200 Aquitania Municipio Municipality
## 2 BoyacÃ¡  201  Arcabuco Municipio Municipality
## 3 BoyacÃ¡  202    BelÃ©n      <NA>         <NA>
## 4 BoyacÃ¡  203    Berbeo Municipio Municipality
## 5 BoyacÃ¡  204 Beteitiva Municipio Municipality

Seleccionemos solo la ciudad capital de este departamento.

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

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

Como ya he descargado y guardado los datos de elevación, necesito leerlos usando la función ráster proporcionada por el paquete ráster. No deberías correr el trozo.

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

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

elevation

## class      : RasterLayer 
## dimensions : 1034, 1033, 1068122  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.545, -71.70425, 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     : -1.570384, 4150.449  (min, max)

plot(elevation, main="Este es el DEM descargado [metros]")
plot(Ruana, 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 = "C:/Users/Brian/Documents/R/win-library/3.6/elevatr", dataType='INT4S', overwrite=TRUE)

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

## class      : RasterLayer 
## dimensions : 1034, 1033, 1068122  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -74.545, -71.70425, 2.796526, 5.629686  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : C:/Users/Brian/Documents/R/win-library/3.6/elevatr.grd 
## names      : layer 
## values     : -2, 4150  (min, max)

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

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

Veamos el nuevo objeto:

elev_crop

## class      : RasterLayer 
## dimensions : 70, 59, 4130  (nrow, ncol, ncell)
## resolution : 0.00275, 0.00274  (x, y)
## extent     : -73.49175, -73.3295, 5.421446, 5.613246  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : layer 
## values     : 2163.53, 3258.219  (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. Vuelva a proyectar los datos de elevación.

Podemos ir a epsg.io y buscar la proyección de la zona MAGNA Colombia Bogotá. Necesitamos obtener la definición de esta referencia espacial en formato PROJ.4 (la que se usa para las bibliotecas sp y raster. Copiemos el texto PROJ.4 y guárdelo.

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)

rep_elev

## class      : RasterLayer 
## dimensions : 212, 180, 38160  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1064897, 1082897, 1091321, 1112521  (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     : 2164.928, 3259.058  (min, max)

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

(rep_Ruana = spTransform(Ruana, spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 1 
## extent      : 1064887, 1082842, 1091256, 1112532  (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   : 5
## names       :  NAME_1, ID_2, NAME_2,    TYPE_2,    ENGTYPE_2 
## value       : BoyacÃ¡,  312,  Tunja, Municipio, Municipality

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

Para evitar un dolor de cabeza, guardemos nuestra DEM:

writeRaster(rep_elev, filename = "C:/Users/Brian/Desktop/Daniela/Boyacá/rep_Ruana_elev", dataType='INT4S', overwrite=TRUE)

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

## class      : RasterLayer 
## dimensions : 212, 180, 38160  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 1064897, 1082897, 1091321, 1112521  (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/Brian/Desktop/Daniela/Boyacá/rep_Ruana_elev.grd 
## names      : layer 
## values     : 2165, 3259  (min, max)

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

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

##   metricas   valores
## 1     mean 2845.1281
## 2      min 2164.9279
## 3      max 3259.0580
## 4       sd  194.8716

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)

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

plot(rep_elev, main="DEM para Tunja [metros]", col=terrain.colors(20, alpha = 1))

plot(slope, main="Pendiente para Tunja [grados]", col=topo.colors(25, alpha = 1))

plot(aspect, main="Aspecto para Tunja [grados]", col=rainbow(25, alpha = 1))

Un plot combinado:

plot(hill, col=grey(1:100/100),legend=FALSE, main="DEM para Tunja", 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)

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 visualización.

Un experto en R sugiere otra forma de visualización aquí. Vamos a intentarlo:

#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/Brian/Pictures/9pvbHjN.jpg")

out = getv(map, aspect, slope)

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

Datos de elevación en R

Marly Daniela Quiroga Buitrago

26/3/2020

Introducción a elevatr

Para acceder a los datos de elevación ráster (por ejemplo, un DEM), el paquete elevatr utiliza los mosaicos de terreno de Amazon Web Services. Exploraremos esta funcionalidad en este cuaderno.

Obtener datos de elevación de ráster

**Usando get_elev_raster () para acceder a los Terrain Tiles en AWS.**

Primero, necesitamos cargar un archivo shape que represente a los municipios de nuestro departamento. En este cuaderno, cargaremos el shapefile descargado del geoportal del DANE según lo solicitado en la clase virtual el 19 de marzo de 2020. Revisemos el contenido de la carpeta:

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

¿Cuáles son los atributos del objeto munic?

Seleccionemos solo la ciudad capital de este departamento.

Como ya he descargado y guardado los datos de elevación, necesito leerlos usando la función ráster proporcionada por el paquete ráster. No deberías correr el trozo.

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

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.

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

Veamos el nuevo objeto:

Reproyectar los datos de elevación

Podemos ir a epsg.io y buscar la proyección de la zona MAGNA Colombia Bogotá. Necesitamos obtener la definición de esta referencia espacial en formato PROJ.4 (la que se usa para las bibliotecas sp y raster. Copiemos el texto PROJ.4 y guárdelo.

Ahora, podemos reproyectar los datos de elevación de las coordenadas geográficas WGS84 en la zona MAGNA Colombia Bogotá.

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

Para evitar un dolor de cabeza, guardemos nuestra DEM:

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

Una exploración rápida de las estadísticas DEM:

Obtención de variables geomorfométricas.

Primero, calcule la pendiente, el aspecto y el sombreado:

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

Un plot combinado:

Mapeo de datos de elevación con sombreador de rayos

Otra forma de visualización.

Un experto en R sugiere otra forma de visualización aquí. Vamos a intentarlo: