#1. Introducción

Este cuaderno ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación (MDE) en R. Su objetivo es ayudar a los estudiantes de Geomática Básica de la Universidad Nacional de Colombia a familiarizarse con los MDE y las variables geomorfométricas.

Para la gestión de datos geoespaciales, existen dos paquetes principales en R:

la biblioteca de características simples (sf), que proporciona funcionalidades para datos vectoriales

la biblioteca terra, que proporciona funcionalidades para datos ráster.

Además, para la cartografía interactiva, una buena opción es la biblioteca leaflet.

Antes de continuar, a continuación se ofrecen algunos consejos para escribir tu propio cuaderno:

Escribe y ejecuta cada fragmento paso a paso

En caso de errores, lee el mensaje, interprétalo e intenta solucionarlo

En caso de que el error persista, pregunta al profesor antes de la fecha límite

Una vez producido el archivo html, es hora de publicarlo en RPubs

Después de publicar el cuaderno, revísalo, corrige la redacción deficiente y vuelve a publicarlo utilizando la opción de actualización

En caso de problemas, revisa la documentación de R

Al crear su cuaderno, asegúrese de agregar tanto texto como sea necesario para: (i) explicar cada fragmento de código; y (ii) describir el resultado correspondiente. Estos dos criterios se tendrán en cuenta para evaluar la calidad de su cuaderno

#2. Configuración Los paquetes necesarios deben instalarse previamente desde la consola.

## SETUP
# INSTALL THESE PACKAGES FROM THE CONSOLE, NOT FROM THIS CHUNK
# paquetes = c("terra","elevatr","sf", "leaflet")
# install.packages(paquetes)

Es recomendable comenzar a limpiar la memoria:

rm(list=ls())

Ahora, carguemos las bibliotecas:

library(elevatr)
## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
## deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
## will be dropped in future versions, so please plan accordingly.
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE
library(leaflet)
library(terra)
## terra 1.8.54

Al ejecutar un fragmento, preste atención a cualquier mensaje que indique un error o advierta sobre posibles conflictos.

Los conflictos pueden deberse a funciones de diferentes paquetes con el mismo nombre (p. ej., terra::extract y tidyr::extract). Para evitar estos problemas, es necesario llamar a varias funciones utilizando la ruta larga, es decir, package_name::function.

#3. 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 modelado ecológico. Obtener acceso a estos datos en R no ha tenido una única 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

Existen varias fuentes de modelos digitales de elevación, como la Shuttle Radar Topography Mission (SRTM), el USGS National Elevation Dataset (NED), Global DEM (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 síntesis de elevación que utiliza los mejores datos de elevación disponibles para una región determinada con un nivel de zoom determinado. Aunque está cerrado, estos datos compilados por Mapzen son actualmente accesibles a través de Terrain Tiles en Amazon Web Services (AWS).

La función principal es get_elev_raster(), que necesita un objeto grpspatial, como un data.frame con las ubicaciones x (longitud) e y (latitud) como las dos primeras columnas, cualquier objeto de entidad simple (sf) o cualquier objeto ráster, y devuelve una RasterLayer de las teselas que se superponen al cuadro delimitador de la entrada. Si se recuperan varias teselas, la salida resultante es una Raster Layer fusionada

Uso de get_elev_raster() para acceder a los mosaicos de terreno en AWS.

Tenga en cuenta que el objeto geoespacial requiere un CRS (es decir, un sistema de referencia de coordenadas). El nivel de zoom (z) predeterminado es 9 (un equilibrio entre la resolución y el tiempo de descarga). Podría empezar probando con un nivel de zoom más alto (por ejemplo, 10).

#4. Obtener datos de elevación para su departamento Primero, cargaremos un shapefile o un geopaquete que represente los municipios de nuestro departamento. En este cuaderno, usaré un geopaquete con los límites oficiales de Córdoba.

Leamos el archivo geoespacial usando una función del paquete sf:

#4.1 Importar límites departamentales Como este cuaderno se encuentra en una carpeta llamada Córdoba, es sencillo obtener una lista de los datos disponibles.

list.files("C:/Users/menju/OneDrive/Documentos/GB2/P1/DATOS")
##  [1] "AGUA"                                 
##  [2] "COL_MUNICIPIOS"                       
##  [3] "ELEVACION"                            
##  [4] "license.txt"                          
##  [5] "MAPAS DPTO"                           
##  [6] "MAPAS QGIS.qgz"                       
##  [7] "MUNI_TOL.gpkg"                        
##  [8] "POBLACION_HUILA.csv"                  
##  [9] "POBLACION_TOL_FINAL.gpkg"             
## [10] "PROYECTO1.qgz"                        
## [11] "PROYECTO2.qgz"                        
## [12] "PT.csv"                               
## [13] "PT.xlsx"                              
## [14] "PTT.csv"                              
## [15] "simplemaps_worldcities_basicv1.77.zip"
## [16] "VÍAS"                                 
## [17] "worldcities.csv"                      
## [18] "worldcities.xlsx"

Ahora, leeremos el geopaquete necesario usando la biblioteca sf:

(munic <-  sf::st_read("C:/Users/menju/OneDrive/Documentos/GB2/P1/DATOS/MUNI_TOL.gpkg"))
## Reading layer `MUNICIPIOS_TOL' from data source 
##   `C:\Users\menju\OneDrive\Documentos\GB2\P1\DATOS\MUNI_TOL.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS:  MAGNA-SIRGAS

Tenga en cuenta el CRS de este conjunto de datos: MAGNA-SIRGAS, sin detalles. Como necesitamos más información sobre el CRS, preguntaremos:

(MUNI_TOL_crs = st_crs(st_geometry(munic)))
## Coordinate Reference System:
##   User input: MAGNA-SIRGAS 
##   wkt:
## GEOGCRS["MAGNA-SIRGAS",
##     DATUM["Marco Geocentrico Nacional de Referencia",
##         ELLIPSOID["GRS 1980",6378137,298.257222101,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
##         BBOX[-4.23,-84.77,15.51,-66.87]],
##     ID["EPSG",4686]]

Tenga en cuenta el atributo que almacena los nombres de los municipios. En mi caso, es mpio_cnmbr. Lo necesitaremos más adelante para el gráfico.

Tenga en cuenta también que, en mi caso, los municipios no tienen ningún atributo con el área correspondiente. Por lo tanto, conviene calcularlo:

# This code calculates area in square meters
# Can you tell us why?
munic$area_m2 = sf::st_area(munic)

Ahora, el área expresada en kilómetros cuadrados:

# This code calculates area in square kilometers
# Can you tell us why?
munic$area = munic$area_m2/1000000
# Now, let's round the area to two digits
munic$area = signif(munic$area, digits = 6)

Volvamos a llamar a nuestros municipios para que se opongan:

munic

#4.2 Obtener centroides de municipios A continuación, crearemos un nuevo objeto con centroides de municipios:

(centers <- st_centroid(munic))
## Warning: st_centroid assumes attributes are constant over geometries

Tenga en cuenta que el atributo geom codifica la geometría del punto. Por lo tanto, debemos dividir cada punto en sus coordenadas x e y:

centers$x = st_coordinates(centers)[,1]
centers$y = st_coordinates(centers)[,2]

Comprobemos el resultado:

centers

#4.3 Obtener los datos de elevación Ahora es el momento de obtener el DEM:

elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
## Mosaicing & Projecting
## Note: Elevation units are in meters.

¿Qué es el objeto de elevación?

elevation
## class      : RasterLayer 
## dimensions : 4092, 3078, 12595176  (nrow, ncol, ncell)
## resolution : 0.000685414, 0.000685414  (x, y)
## extent     : -76.28906, -74.17936, 2.811272, 5.615986  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source     : file19c073d0f52.tif 
## names      : file19c073d0f52

Toma nota de las características del objeto de elevación y comprende cada una. ¿Cuál es su resolución espacial (en metros)?

Observa algunos aspectos de 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 ráster.

resolution: 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.

extension: la extensión espacial del ráster. Este valor estará en las mismas unidades de coordenadas que el sistema de referencia de coordenadas del ráster.

crs: la cadena del sistema de referencia de coordenadas del ráster. Este ráster está en coordenadas geográficas con un datum WGS 84.

En tu cuaderno, anota el tamaño de celda de tu ráster (es decir, la resolución espacial), en metros. Tenga en cuenta que la elevación es una capa ráster temporal (es decir, solo existe en memoria). Sin embargo, es posible escribir el DEM en su directorio local mediante la función writeRaster de la biblioteca rgdal:

writeRaster(elevation, "C:/Users/menju/OneDrive/Documentos/GB2/Rstudio2/elev_ELEVACION_z10.tif", overwrite=TRUE)
(elevation2 <- terra::rast(elevation))
## class       : SpatRaster 
## size        : 4092, 3078, 1  (nrow, ncol, nlyr)
## resolution  : 0.000685414, 0.000685414  (x, y)
## extent      : -76.28906, -74.17936, 2.811272, 5.615986  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source      : file19c073d0f52.tif 
## name        : file19c073d0f52

5

pal <- colorNumeric(c("cyan", "forestgreen","yellow","tan","orange", "brown"), values(elevation),
  na.color = "transparent")
(elevation3 <- terra::aggregate(elevation2, fact = 2))
## |---------|---------|---------|---------|=========================================                                          
## class       : SpatRaster 
## size        : 2046, 1539, 1  (nrow, ncol, nlyr)
## resolution  : 0.001370828, 0.001370828  (x, y)
## extent      : -76.28906, -74.17936, 2.811272, 5.615986  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source(s)   : memory
## name        : file19c073d0f52 
## min value   :          -38.75 
## max value   :         5362.50
elevation4 <- terra::crop(elevation3, munic, mask=TRUE)
leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% 
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
                          "Km2: ", munic$area, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~mpio_cnmbr,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE, textsize = "10px")) %>%
  addRasterImage(elevation4, colors = pal, opacity = 0.9)  %>%
  addLegend(pal = pal, values = values(elevation),
    title = "Elevation data for Tolima (mts)")
## Warning: sf layer has inconsistent datum (+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs).
## Need '+proj=longlat +datum=WGS84'