R Notebook

Erika Aparicio

02/04/2020

Datos de elevación en R

1. 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.

library(rasterVis)

Loading required package: raster
Loading required package: sp
Loading required package: lattice
Loading required package: latticeExtra

library(raster)
library(rgl)

Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

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

2. Obtenga 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.

se mencionó, un marco de datos 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 marco de datos requiere un prj. Mostraré ejemplos usando un SpatialPolygonsDataFrame. El nivel de zoom (z) predeterminado es 9 (una compensación entre la resolución y el tiempo de descarga), pero no pude obtener nada con este zoom. Tuve que usar un zoom inferior igual a 8. De lo contrario, el portátil RStudio Cloud se bloquearía una y otra vez. Muy a menudo, tuve que “volver a cargar” la página y comenzarla de nuevo.

Primero, se necesita cargar un archivo shape que represente a los municipios del departamento Meta. En este cuaderno, cargaremos el archivo de forma descargado del geoportal del DANE.

#####Revisemos el contenido de la carpeta:

list.files("C:/Users/oaparicio15/Desktop/QGIS/ColombiaDatos/50_META/ADMINISTRATIVO")

 [1] "MGN_DPTO_POLITICO.cpg"                            
 [2] "MGN_DPTO_POLITICO.dbf"                            
 [3] "MGN_DPTO_POLITICO.prj"                            
 [4] "MGN_DPTO_POLITICO.sbn"                            
 [5] "MGN_DPTO_POLITICO.sbx"                            
 [6] "MGN_DPTO_POLITICO.shp"                            
 [7] "MGN_DPTO_POLITICO.shp.DG_EST118.1768.8604.sr.lock"
 [8] "MGN_DPTO_POLITICO.shp.xml"                        
 [9] "MGN_DPTO_POLITICO.shx"                            
[10] "MGN_MPIO_POLITICO.cpg"                            
[11] "MGN_MPIO_POLITICO.dbf"                            
[12] "MGN_MPIO_POLITICO.prj"                            
[13] "MGN_MPIO_POLITICO.sbn"                            
[14] "MGN_MPIO_POLITICO.sbx"                            
[15] "MGN_MPIO_POLITICO.shp"                            
[16] "MGN_MPIO_POLITICO.shp.DG_EST118.1768.8604.sr.lock"
[17] "MGN_MPIO_POLITICO.shp.xml"                        
[18] "MGN_MPIO_POLITICO.shx"                            

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

¿Cuáles son los atributos del objeto munic?

Seleccionemos solo la ciudad capital del departamento (Meta)

Descomente la siguiente línea para descargar los datos de elevación:

#elevation <- get_elev_raster(vila, z = 8)

elevation <- get_elev_raster(villa, z = 8)

Downloading DEMs [========>------------------]  33% eta:  2s
Downloading DEMs [=============>-------------]  50% eta:  3s
Downloading DEMs [=================>---------]  67% eta:  2s
Downloading DEMs [=====================>-----]  83% eta:  1s
Downloading DEMs [===========================] 100% eta:  0s
Merging DEMs
Reprojecting DEM to original projection
ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -InfNote: 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 deberia correrlo

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

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

elevation

class      : RasterLayer 
dimensions : 1546, 1033, 1597018  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -74.545, -71.70425, 1.393646, 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)

3. 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/oaparicio15/Documents/R/win-library/3.6/elevatr", datatype='INT4S', overwrite=TRUE)

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

Ahora, recortemos los datos de elevación correspondientes a Villavicencio

Veamos el nuevo objeto:

elev_crop

class      : RasterLayer 
dimensions : 129, 216, 27864  (nrow, ncol, ncell)
resolution : 0.00275, 0.00274  (x, y)
extent     : -73.7695, -73.1755, 3.936366, 4.289826  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
source     : memory
names      : layer 
values     : 203.7742, 3577.458  (min, max)

4. 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.

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)

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf

rep_elev

class      : RasterLayer 
dimensions : 391, 660, 258060  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1034192, 1100192, 927079.8, 966179.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 
source     : memory
names      : layer 
values     : 201.3911, 3585.355  (min, max)

Ahora, reproyectemos el SpatialPolygonsDataFrame que representa la capital del departamento del Meta:

(rep_villa = spTransform(villa,spatialref))

class       : SpatialPolygonsDataFrame 
features    : 1 
extent      : 1034342, 1100224, 926940.8, 965985.1  (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       :         50,      50001, VILLAVICENCIO,       1840, 1285.93087074,      2017,       META, 2.03977418027, 0.104718109413 

Está tramando el tiempo:

Para evitar un dolor de cabeza, guardemos nuestra DEM

writeRaster(rep_elev, filename = "C:/Users/oaparicio15/Desktop/Rstudio/rep_villa_elev", dataType='INT4S', overwrite=TRUE)

argument "datatype" misspelled as "dataType"

class      : RasterLayer 
dimensions : 391, 660, 258060  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 1034192, 1100192, 927079.8, 966179.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 
source     : C:/Users/oaparicio15/Desktop/Rstudio/rep_villa_elev.grd 
names      : layer 
values     : 201, 3585  (min, max)

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

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

crear vector unidimensional

creación de marco de datos con estadísticas de elevación [metros]

6. Obtención de variables geomorfométricas.

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

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

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

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

Una trama combinada:

7. 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 visualizaciones de datos 3D.

#install.packages("rayshader")

Convertir el DEM en una matriz:

elmat = raster_to_matrix(rep_elev)

[1] "Dimensions of matrix are: 660x391."

Usamos otra de las texturas incorporadas del sombreador de rayos:

######detect_water y add_water agrega una capa de agua al mapa:

Y también podemos agregar una capa de trazado de rayos desde esa dirección del sol:

8. Otra forma de visualización.

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

Vamos a intentarlo:

#install.packages("jpeg")

Función de sombreado de mapeo de entorno esférico:

out = getv(map, aspect, slope)

ning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Infning昼㹡n argumento finito para min; retornando Infningun argumento finito para max; retornando -Inf