Atributos geomorfométricos del terreno en R
Andres Felipe Ramirez Ramirez
23 - 07 - 2026
1. Introducción
Este cuaderno ilustra cómo calcular atributos geomorfométricos del terreno a partir de modelos digitales de elevación en cuadrícula. Su objetivo es enseñar conceptos y herramientas básicas de geoinformática a estudiantes de agronomía de la Universidad Nacional de Colombia.
Para calcular los atributos del terreno, utilizaremos el paquete Multiscale DTM de R.
Al crear tu cuaderno, asegúrate de incluir todo el texto necesario para: (I) explicar cada fragmento de código (II) describir el resultado correspondiente. Estos dos criterios se tendrán en cuenta para evaluar la calidad de tu cuaderno.
2. Setup
Los paquetes necesarios deben instalarse previamente desde la consola.
## CONFIGURACIÓN
# INSTALE ESTE PAQUETE DESDE LA CONSOLA, NO DESDE ESTE FRAGMENTO
# paquetes = c("MultiscaleDTM", "exactextractr")
# install.packages(paquetes)Comencemos por limpiar la memoria:
rm(list=ls())Ahora, carguemos las bibliotecas necesarias:
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)Al ejecutar un fragmento de código, 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
nombres idénticos (por ejemplo, terra::extract y
tidyr::extract). Para evitar estos problemas, es necesario
llamar a las funciones de forma completa, es decir, utilizando el
formato nombre_paquete::función.
3. Introducción a Multiscale DTM
Este paquete de R calcula atributos geomorfométricos del terreno a múltiples escalas a partir de modelos digitales del terreno (MDT; es decir, rásteres de elevación o batimetría) con cuadrícula regular, mediante una ventana específica.
4. Lectura de datos de entrada
Primero, cargaremos un modelo digital de elevación (MDE) ráster para nuestro departamento. En este cuaderno, usaremos un MDE obtenido previamente con este mismo cuaderno.
Leamos el MDE con Terra:
# Cambie la ruta y el nombre del archivo DEM según sea necesario.
(dem = terra::rast("elev_zona_cafetera_z10.tif"))## class : SpatRaster
## size : 3067, 2566, 1 (nrow, ncol, nlyr)
## resolution : 0.0006851336, 0.0006851336 (x, y)
## extent : -76.28906, -74.53101, 3.864449, 5.965754 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : elev_zona_cafetera_z10.tif
## name : elev_zona_cafetera_z10
## min value : -299
## max value : 5282Vamos a reducir la resolución del DEM para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Ahora, leeremos la información de nuestros municipios utilizando la biblioteca sf:
(munic <- sf::st_read("zonacafetera.gpkg") |>
dplyr::select(
dpto_ccdgo,
mpio_ccdgo,
mpio_cdpmp,
mpio_cnmbr,
geom
))## Reading layer `zonacafetera_munic' from data source
## `/Users/felipe/Documents/GB2/Proyecto4/Datos/zonacafetera.gpkg'
## using driver `GPKG'
## Simple feature collection with 53 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.21143 ymin: 4.075311 xmax: -74.62746 ymax: 5.779182
## Geodetic CRS: MAGNA-SIRGASTenga en cuenta que ambos conjuntos de datos están en coordenadas geográficas.
Ahora, recortemos el objeto dem2 a los límites de nuestro departamento:
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Transformación de coordenadas (10%)
Para calcular los atributos del terreno, el modelo digital de elevación (MDE) debe estar en coordenadas planas. Como recordará, el IGAC cambió en 2020 a un sistema de coordenadas planas con origen único nacional. ¿Cuál es su código EPSG?
Modifique los siguientes dos fragmentos de código para transformar los datos de entrada al sistema de coordenadas nacional de Colombia.
# Consulta los parámetros de esta función de Terra aquí
# https://rdrr.io/github/rspatial/terra/man/project.html
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1247, 1162, 1 (nrow, ncol, nlyr)
## resolution : 151.752, 151.752 (x, y)
## extent : 4643555, 4819891, 2008382, 2197617 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_zona_cafetera_z10
## min value : 143.713638
## max value : 5272.071289También necesitamos realizar una transformación similar para el objeto vectorial:
# Compruebe aquí los parámetros de esta función sf
# https://github.com/rstudio/cheatsheets/blob/main/sf.pdf
(munic_plane = sf::st_transform(munic, "EPSG:9377"))6. Cálculo de atributos del terreno (20%)
Veamos cómo calcular la pendiente y la orientación.
Ve a la consola y escribe ?SlpAsp. A continuación, en la ventana de ayuda de la derecha, se mostrará la información correspondiente:
Ahora, llamemos a la función SlpAsp:
# w se define como el tamaño de la ventana de análisis (o vecindario) que se utiliza para calcular los atributos del terreno.
# method es un parametro que indica cómo se consideran a los vecinos dentro de esta ventana, es llamado método que en porque se considera a la celda del centro como reina y al igual que la reina del ajedrez se mueve en las direcciones: arriba, abajo, izquíerda, derecha y las 4 diagonales.
# Change if needed
(slp_asp = MultiscaleDTM::SlpAsp(
dem_plane,
w = c(3, 3),
unit = "degrees",
method = "queen",
metrics = c("slope", "aspect"),
na.rm = TRUE,
include_scale = FALSE,
mask_aspect = TRUE
))## class : SpatRaster
## size : 1247, 1162, 2 (nrow, ncol, nlyr)
## resolution : 151.752, 151.752 (x, y)
## extent : 4643555, 4819891, 2008382, 2197617 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.003857, 0.003072
## max values : 81.275057, 359.999127Ahora, dividamos el objeto slp_asp en dos partes.
Capa 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1247, 1162, 1 (nrow, ncol, nlyr)
## resolution : 151.752, 151.752 (x, y)
## extent : 4643555, 4819891, 2008382, 2197617 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.003857
## max value : 81.275057¿Cuál es la distribución de los valores de pendiente?
terra::hist(slope,
main = "Pendiente Zona Cafetera",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 69% of the cells was used (of which 61% was NA)
Capa 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1247, 1162, 1 (nrow, ncol, nlyr)
## resolution : 151.752, 151.752 (x, y)
## extent : 4643555, 4819891, 2008382, 2197617 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.003072
## max value : 359.999127terra::hist(aspect,
main = "Aspecto Zona Cafetera",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 69% of the cells was used (of which 61% was NA)
Además, vamos a convertir los grados de pendiente a porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1247, 1162, 1 (nrow, ncol, nlyr)
## resolution : 151.752, 151.752 (x, y)
## extent : 4643555, 4819891, 2008382, 2197617 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.006732
## max value : 651.605618terra::hist(slope_perc,
main = "Pendiente Zona Cafetera",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 69% of the cells was used (of which 61% was NA)
#terra::hist(slope_perc)7. Cálculo de estadísticas zonales (20%)
Una operación de estadísticas zonales consiste en calcular estadísticas sobre los valores de las celdas de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En nuestro caso, nos interesa calcular el promedio de los valores de pendiente por municipio.
Primero, veamos la clasificación de pendientes del IGAC para fines agronómicos:
Ahora, ¡manos a la obra!
# Reclasificar el ráster de pendiente
#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
75, 160, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)Ahora, calculemos las estadísticas zonales usando exactextractr:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG: 9377) .
## To avoid on-the-fly coordinate transformation, ensure that st_crs() returns the
## same result for both raster and polygon inputs.## | | | 0% | |= | 2% | |=== | 4% | |==== | 6% | |===== | 8% | |======= | 9% | |======== | 11% | |========= | 13% | |=========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================= | 25% | |================== | 26% | |==================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |========================== | 38% | |============================ | 40% | |============================= | 42% | |============================== | 43% | |================================ | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |====================================== | 55% | |======================================== | 57% | |========================================= | 58% | |========================================== | 60% | |============================================ | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |================================================== | 72% | |==================================================== | 74% | |===================================================== | 75% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |=========================================================== | 85% | |============================================================= | 87% | |============================================================== | 89% | |=============================================================== | 91% | |================================================================= | 92% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%## [1] 24.700743 30.159378 23.075106 31.184830 23.271345 16.382748 26.344486
## [8] 4.198874 34.086853 36.136024 38.867756 34.436546 38.440758 27.183477
## [15] 15.886221 35.445564 17.174583 37.743206 24.032919 22.917028 30.877117
## [22] 30.857990 20.704926 28.687483 11.412414 28.958290 9.316138 4.428427
## [29] 23.261639 22.019220 7.772769 32.301636 10.626133 39.925201 6.030480
## [36] 5.754775 36.760571 6.032310 31.171371 13.832696 31.135393 20.318336
## [43] 24.945431 17.921057 22.805515 34.493271 14.171336 21.902081 41.084381
## [50] 37.790211 28.543644 25.663641 33.206814hist(munic$mean_slope,
main = "Pendiente media Zona Cafetera",
xlab = "Slope (in percentage)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG: 9377) .
## To avoid on-the-fly coordinate transformation, ensure that st_crs() returns the
## same result for both raster and polygon inputs.## | | | 0% | |= | 2% | |=== | 4% | |==== | 6% | |===== | 8% | |======= | 9% | |======== | 11% | |========= | 13% | |=========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================= | 25% | |================== | 26% | |==================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |========================== | 38% | |============================ | 40% | |============================= | 42% | |============================== | 43% | |================================ | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |====================================== | 55% | |======================================== | 57% | |========================================= | 58% | |========================================== | 60% | |============================================ | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |================================================== | 72% | |==================================================== | 74% | |===================================================== | 75% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |=========================================================== | 85% | |============================================================= | 87% | |============================================================== | 89% | |=============================================================== | 91% | |================================================================= | 92% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%## [1] 5 5 5 5 5 4 5 1 5 5 5 5 5 5 2 5 4 5 4 5 5 5 5 5 2 5 1 1 5 5 2 5 3 5 1 1 5 1
## [39] 5 2 5 5 5 4 4 5 4 4 5 5 5 4 5hist(munic$class,
main = "Pendiente reclasificada de Zona Cafetera",
xlab = "Slope (as a category)")
#terra::hist(slope_perc)Vamos a transformar la pendiente de nuevo en coordenadas geográficas:
# Esta es la pendiente reclasificada
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1253, 1165, 1 (nrow, ncol, nlyr)
## resolution : 0.00137025, 0.00137025 (x, y)
## extent : -76.21916, -74.62282, 4.069227, 5.786149 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 542.183655# Esta es la pendiente porcentual
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1253, 1165, 1 (nrow, ncol, nlyr)
## resolution : 0.00137025, 0.00137025 (x, y)
## extent : -76.21916, -74.62282, 4.069227, 5.786149 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.05912
## max value : 542.1836558. Trazado de la pendiente del terreno (20%)
Primero, necesitaremos una paleta de colores para la pendiente:
# Amplía el número de colores para mejorar tu visualización
# https://r-graph-gallery.com/42-colors-names.html
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ahora es el momento de graficar.
leaflet(munic) %>% addTiles() %>% setView(-75.5, 5.0, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Pendiente del terreno en Zona Cafetera (%)")## 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'Haz clic en diferentes sitios para obtener los nombres de los municipios y los valores de pendiente porcentual más frecuentes en cada uno.
9. Visualización alternativa de la pendiente del terreno (15%)
Visualice el terreno inclinado en grados y añada etiquetas que muestren los nombres de los municipios junto con sus correspondientes valores de pendiente media. Mantenga los valores emergentes tal como se muestran en el gráfico anterior.
# Calcular el centro de cada municipio
centros <- sf::st_centroid(munic)## Warning: st_centroid assumes attributes are constant over geometries# Paleta de colores para la pendiente en grados
palSlope <- colorNumeric(
palette = c("darkseagreen3","yellow","orange","brown","darkred"),
domain = values(slope),
na.color = "transparent"
)
# Mapa
leaflet(munic) %>%
addTiles() %>%
setView(-75.5, 5.0, 9) %>%
addRasterImage(slope,
colors = palSlope,
opacity = 0.8) %>%
addPolygons(
color = "gray",
weight = 1,
fillOpacity = 0.1,
popup = paste(
"Municipio: ", munic$mpio_cnmbr,
"<br>",
"Clase de pendiente: ", munic$class
)
) %>%
addLabelOnlyMarkers(
data = centros,
label = paste0(
munic$mpio_cnmbr,
"\n",
round(munic$mean_slope,1),
"%"
),
labelOptions = labelOptions(
noHide = TRUE,
direction = "center",
textOnly = TRUE,
style = list(
"font-size" = "10px",
"font-weight" = "bold"
)
)
) %>%
addLegend(
pal = palSlope,
values = values(slope),
title = "Pendiente (°)"
)## 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'## 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'