2.Introducción a elevatr
Los datos de elevación se utilizan para una amplia gama de
aplicaciones, incluyendo, por ejemplo, visualización, hidrología y
modelización ecológica. El acceso a estos datos en R no ha tenido una
interfaz única, está disponible a través de funciones en muchos
paquetes, o requiere acceso local a los datos. Esto ya no es necesario,
puesto que ahora existen diversas API que proporcionan acceso
programático a los datos de elevación. El paquete elevatr se creó para
estandarizar el acceso a los datos de elevación desde las API web.
Para acceder a los datos de elevación raster (por ejemplo, un DEM)
el paquete elevatr utiliza los Terrain Tiles de Amazon Web Services.
Exploraremos esta funcionalidad en este cuaderno.
Existen varias fuentes de modelos digitales de elevación, como el
Shuttle Radar Topography Mission (SRTM), el USGS National Elevation
Dataset (NED), el Global DEM (GDEM) y otros. Cada uno de estos DEM tiene
pros y contras 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 determinada en un nivel de zoom determinado. Aunque cerrados,
estos datos compilados por Mapzen son actualmente accesibles a través de
los Terrain Tiles en Amazon Web Services (AWS).
La entrada para get_elev_raster() puede ser un data.frame con
ubicaciones x (longitud) e y (latitud) como las dos primeras columnas,
cualquier objeto sp, o cualquier objeto raster y devuelve un RasterLayer
de los mosaicos que se solapan con el cuadro delimitador de la entrada.
Si se recuperan varios mosaicos, el resultado es una capa raster
combinada.
** Utilizando get_elev_raster() para acceder a los mosaicos del
terreno en AWS**.
Como se ha mencionado, un marco de datos con columnas x e y, un
objeto sp, o un objeto raster necesita ser la entrada y el src necesita
ser establecido a “mapzen” (este es el valor por defecto).
No hay diferencia en el uso de los tipos de datos de entrada sp y
raster.
El marco de datos requiere un CRS (es decir, un sistema de
referencia de coordenadas). El nivel de zoom (z) por defecto es 9 (un
compromiso entre resolución y tiempo de descarga). Puede empezar
probando con un nivel de zoom superior (por ejemplo, 10).
3. Obtenga datos de elevación raster para su departamento
Primero, necesitamos cargar un shapefile o un geopackage que
represente municipios de nuestro departamento. En este cuaderno
utilizaré un shapefile descargado del geoportal del DANE.
Vamos a leer el shapefile utilizando una función proporcionada por
el paquete sf:
munic <- st_read("C:/Users/Dasil/Desktop/Trabajos Universidad/geomatica/DEM/47_MAGDALENA/ADMINISTRATIVO/MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source
## `C:\Users\Dasil\Desktop\Trabajos Universidad\geomatica\DEM\47_MAGDALENA\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 30 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.9466 ymin: 8.936489 xmax: -73.54184 ymax: 11.34891
## Geodetic CRS: WGS 84
## Comprobemos la geometría, el cuadro delimitador, el CRS y los atributos del objeto munic:
head(munic)
## Simple feature collection with 6 features and 9 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.9466 ymin: 9.618984 xmax: -73.82743 ymax: 10.39318
## Geodetic CRS: WGS 84
## DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR
## 1 47 47541 PEDRAZA
## 2 47 47030 ALGARROBO
## 3 47 47058 ARIGUANÍ
## 4 47 47161 CERRO DE SAN ANTONIO
## 5 47 47170 CHIVOLO
## 6 47 47205 CONCORDIA
## MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
## 1 Decreto dptal 1322 del 19 dediciembre de 1908 325.4799 2017 MAGDALENA
## 2 Ordenanza 008 de Junio 24 de 1999 406.3884 2017 MAGDALENA
## 3 Ordenanza 14 Bis de Noviembre 30 de 1966 1131.7442 2017 MAGDALENA
## 4 Ordenanza 3038 de 1912 177.1961 2017 MAGDALENA
## 5 Decreto 107 de Marzo 8 de 1974 536.5437 2017 MAGDALENA
## 6 Ordenanza 007 del 24 de Junio de 1999 109.4837 2017 MAGDALENA
## Shape_Leng Shape_Area geometry
## 1 0.9740441 0.026847311 POLYGON ((-74.8832 10.21502...
## 2 1.2656800 0.033536941 POLYGON ((-74.10861 10.3915...
## 3 1.9928810 0.093278409 POLYGON ((-74.03381 9.99541...
## 4 0.7003465 0.014622408 POLYGON ((-74.85084 10.367,...
## 5 1.3263733 0.044254449 POLYGON ((-74.49589 10.2315...
## 6 0.5443263 0.009033178 POLYGON ((-74.81188 10.2748...
# Creemos un nuevo objeto con los centroides de los municipios
## La tarea es encontrar un punto central para cada región
## Convertimos las características simples en polígonos espaciales
sp.munic = as_Spatial(munic)
# Ahora usamos la función *gCentroid* del paquete rgeos
centers <- data.frame(gCentroid(sp.munic, byid = TRUE))
centers$region <- row.names(sp.munic)
centers$munic <- sp.munic$MPIO_CNMBR
## Ahora, vamos a descargar los datos de elevación de nuestro departamento:
# z indica el nivel de zoom de los datos
# a mayor zoom mayor resolución espacial
#elevación <- get_elev_raster(munic, z = 12)
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
#¿Cuál es el objeto de elevación?
elevation
## class : RasterLayer
## dimensions : 4078, 2589, 10557942 (nrow, ncol, ncell)
## resolution : 0.0006789017, 0.0006789017 (x, y)
## extent : -75.23438, -73.4767, 8.754526, 11.52309 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : file451c58f93eac.tif
## names : file451c58f93eac
Fíjese en algunas cosas sobre este DEM:
Dimensions: el “tamaño” del archivo en píxeles
nrow, ncol: el número de filas y columnas de los datos (imagina una
hoja de cálculo o una matriz).
ncells: el número total de píxeles o celdas que componen el
raster.
Resolución: el tamaño de cada píxel (en grados decimales en este
caso). Recordemos que 1 grado decimal representa aproximadamente 111,11
km en el Ecuador.
Extensión: la extensión espacial de la trama. Este valor estará en
las mismas unidades de coordenadas que el sistema de referencia de
coordenadas del raster.
crs: la cadena del sistema de referencia de coordenadas del ráster.
Este ráster está en coordenadas geográficas con un datum de WGS 84.
Tenga en cuenta que la elevación es una capa ráster temporal (es
decir, sólo existe en memoria). Sin embargo, es posible escribir el MDE
en el directorio local utilizando la función writeRaster de la
biblioteca rgdal:
4. Visualizar los datos de elevación
A efectos de visualización, una resolución más baja agiliza la
tarea:
#Este chunk es opcional
#Usalo solo cuando los datos del zoom sean muy altos
elevation2 <- aggregate(elevation, fact=4, fun=mean)
crs(elevation2) <- "+proj=longlat +datum=WGS84"
Una buena paleta es clave para una visualización adecuada. Probemos
con ésta:
pal <- colorNumeric(c("forestgreen","yellow","tan","brown"), values(elevation2),
na.color = "transparent")
Ahora, utilizaremos la biblioteca leaflet para ver los datos de
elevación:
leaflet(munic) %>% addTiles() %>%
addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
opacity = 0.25, fillOpacity = 0.15,
popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
"Km2: ", munic$MPIO_NAREA, "<br>")) %>%
addLabelOnlyMarkers(data = centers,
lng = ~x, lat = ~y, label = ~region,
labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
addRasterImage(elevation2, colors = pal, opacity = 0.9) %>%
addLegend(pal = pal, values = values(elevation2),
title = "Elevation data for Magdalena (mts)")
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
Haga clic en el ID de cada región para obtener el nombre y la
extensión del municipio. Cambia las opciones de trazado para que el mapa
sea más legible.
5. Reproyección de los datos de elevación
Cuando se trabaja con DEMs, siempre es una buena idea utilizar
coordenadas cartográficas en lugar de coordenadas geográficas. Esto se
debe al hecho de que, en las coordenadas geográficas, las unidades de
dimensión horizontal son los grados decimales PERO la unidad de
dimensión vertical son los metros.
Para reproyectar los datos de elevación, realizamos dos pasos.
En primer lugar, definimos el CRS objetivo:
#library(sp)
# WGS 84 / UTM zone 18N
spatialref <- CRS(SRS_string="EPSG:32618")
A continuación, reproyectamos el DEM utilizando la función
projectRaster del paquete raster:
# Project Raster
# First, we create one raster object
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elevation, crs=spatialref)
# Adjust the cell size to your cell size
res(pr3) <- 80
# now project
rep_elev <- projectRaster(elevation, pr3)
¿Qué tenemos?
rep_elev
## class : RasterLayer
## dimensions : 3832, 2417, 9261944 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 474221.7, 667581.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : file451c58f93eac
## values : -2439.451, 5664.858 (min, max)
Ahora, reproyectemos el SpatialPolygonsDataFrame que representa los
municipios de nuestro departamento.
(rep_munic = spTransform(sp.munic,spatialref))
## class : SpatialPolygonsDataFrame
## features : 30
## extent : 505849.9, 659441.8, 988021.7, 1254730 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +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
## min values : 47, 47001, ALGARROBO, 1525, 109.48370634, 2017, MAGDALENA, 0.544326259962, 0.00903317812539
## max values : 47, 47980, ZONA BANANERA, Ordenanza 74 de 1912, 2347.13929515, 2017, MAGDALENA, 3.19741434448, 0.194233330475
7. Obtención de variables geomorfométricas
Antes de continuar, asegúrese de haber comprendido los conceptos
básicos sobre: DEM, pendiente y aspecto.
En primer lugar, calcula la pendiente, el aspecto y el relieve:
slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)
Comprobemos el primer resultado:
slope
## class : RasterLayer
## dimensions : 3832, 2417, 9261944 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 474221.7, 667581.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : slope
## values : 0, 76.59431 (min, max)
Ahora el siguiente resultado
aspect
## class : RasterLayer
## dimensions : 3832, 2417, 9261944 (nrow, ncol, ncell)
## resolution : 80, 80 (x, y)
## extent : 474221.7, 667581.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : aspect
## values : 0, 360 (min, max)
Para trazar la pendiente, es una buena idea reducir el tamaño del
RasterLayer:
#This chunk is optional
#Use it only when zoom data was very high
slope2 <- aggregate(slope, fact=4, fun=mean)
Puede ser conveniente proyectar el conjunto de datos de
elevación:
slope3 <- projectRasterForLeaflet(slope2, "ngb")
Vamos a crear una paleta de colores:
pal2 <- colorNumeric(c("#537188", "#CBB279", "#E1D4BB", "#EEEEEE"), values(slope3),na.color = "transparent")
Ahora plot time!
leaflet(munic) %>% addTiles() %>%
addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
opacity = 0.25, fillOpacity = 0.15,
popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
"Km2: ", munic$MPIO_NAREA, "<br>")) %>%
addLabelOnlyMarkers(data = centers,
lng = ~x, lat = ~y, label = ~region,
labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
addRasterImage(slope3, colors = pal2, opacity = 1.0) %>%
addLegend(pal = pal2, values = values(slope3),
title = "Slope for Magdalena (%)")
Recordemos qué es el objeto slope2:
slope2
## class : RasterLayer
## dimensions : 958, 605, 579590 (nrow, ncol, ncell)
## resolution : 320, 320 (x, y)
## extent : 474221.7, 667821.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## crs : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## source : memory
## names : slope
## values : 0, 61.90718 (min, max)
Tenga en cuenta que existen dos bibliotecas para gestionar datos
raster en R: - raster, y - terra
La biblioteca raster es bien conocida (pero, en algunas tareas,
bastante lenta). La biblioteca terra es más reciente (pero más
rápida).
Vamos a convertir un RasterLayer (es decir, un objeto raster) en un
SpatRaster (es decir, un objeto terra):
(nslope2 <- as(slope2, "SpatRaster"))
## class : SpatRaster
## dimensions : 958, 605, 1 (nrow, ncol, nlyr)
## resolution : 320, 320 (x, y)
## extent : 474221.7, 667821.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 61.90718
La clasificación de las pendientes es una forma común de entender la
variabilidad del relieve. Es necesario incluir aquí una tabla con los
intervalos de pendiente, los nombres de las clases y la descripción, tal
como lo define el Instituto Geográfico Colombiano (IGAC).
Vamos a realizar una clasificación de pendientes basada en
intervalos comunes. Primero, creamos una matriz de reclasificación.
m <- c(0, 3, 1, 3, 7, 2, 7,12,3, 12,
25, 4, 25, 50, 5, 50, 75, 6)
(rclmat <- matrix(m, ncol=3, byrow = TRUE))
## [,1] [,2] [,3]
## [1,] 0 3 1
## [2,] 3 7 2
## [3,] 7 12 3
## [4,] 12 25 4
## [5,] 25 50 5
## [6,] 50 75 6
Ahora, la tarea de reclasificación:
# for values >= 0 (instead of > 0), do
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)
Comprobemos el resultado:
slope_rec
## class : SpatRaster
## dimensions : 958, 605, 1 (nrow, ncol, nlyr)
## resolution : 320, 320 (x, y)
## extent : 474221.7, 667821.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0
## max value : 6
9. Visualización interactiva
Ahora, queremos crear un mapa interactivo de pendientes.
Convirtamos el objeto sf en un objeto SpatVector:
(nmunic <- as(rep_munic, "SpatVector"))
## class : SpatVector
## geometry : polygons
## dimensions : 30, 9 (geometries, attributes)
## extent : 505849.9, 659441.8, 988021.7, 1254730 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
## names : DPTO_CCDGO MPIO_CCDGO MPIO_CNMBR MPIO_CRSLC MPIO_NAREA
## type : <chr> <chr> <chr> <chr> <num>
## values : 47 47541 PEDRAZA Decreto dptal ~ 325.5
## 47 47030 ALGARROBO Ordenanza 008 ~ 406.4
## 47 47058 ARIGUANÍ Ordenanza 14 B~ 1132
## MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
## <int> <chr> <num> <num>
## 2017 MAGDALENA 0.974 0.02685
## 2017 MAGDALENA 1.266 0.03354
## 2017 MAGDALENA 1.993 0.09328
Convirtamos el objeto ráster en un objeto terra:
(nslope <- as(slope, "SpatRaster"))
## class : SpatRaster
## dimensions : 3832, 2417, 1 (nrow, ncol, nlyr)
## resolution : 80, 80 (x, y)
## extent : 474221.7, 667581.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 76.59431
Ahora, calculemos las categorías de pendiente agregadas por
municipio:
## S4 method for signature 'SpatRaster,SpatVector'
munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)
Tenga en cuenta que el resultado es un objeto terra:
munic.slope
## class : SpatRaster
## dimensions : 3832, 2417, 1 (nrow, ncol, nlyr)
## resolution : 80, 80 (x, y)
## extent : 474221.7, 667581.7, 967699.7, 1274260 (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618)
## source(s) : memory
## name : slope
## min value : 0.2804244
## max value : 20.5705299
Ahora veamos el mapa
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
tm_shape(munic.slope) +
tm_raster(style="cont", n=7, alpha=0.6, title="Mean slope")+
tm_layout(legend.outside = TRUE)
## stars object downsampled to 794 by 1259 cells. See tm_shape manual (argument raster.downsample)
