Geomorphometric terrain attributes in R
Mateo Bustos
20/07/2026
1.Introduccion
Este cuaderno ilustra cómo calcular atributos de terreno geomorfométricos 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 los estudiantes de pregrado de agronomía de la Universidad Nacional de Colombia.
Para calcular los atributos del terreno, utilizaremos el paquete de R MultiscaleDTM.
Al crear tu cuaderno, asegúrate de agregar tanto texto como sea necesario para: (i) explicar cada bloque de código; y (ii) describir el resultado correspondiente. Estos dos criterios serán considerados para evaluar la calidad de tu cuaderno.
2.Setup
Limpiamos memoria
rm(list=ls())Cargamos las librerias que se habian instalado
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)3. Introduction to MultiscaleDTM
Este paquete de R calcula atributos geomorfométricos del terreno a múltiples escalas a partir de modelos digitales de terreno en cuadrícula regular (como rásteres de elevación o batimetría) mediante el uso de una ventana móvil especificada.
Funciones Principales
Pendiente, Orientación y Curvatura (Slope, Aspect and Curvature)
- La función
SlpAspcalcula la pendiente y la orientación a múltiples escalas utilizando un subconjunto de celdas dentro de una ventana focal (por ejemplo, ventanas de 3x3 o de 5x5) para analizar las características topográficas del terreno.
4.Leemos los datos que tenemos
Primero, cargaremos un MDE (DEM) ráster para nuestro departamento. En este cuaderno, utilizaré un MDE obtenido previamente utilizando este cuaderno.
Leamos el MDE usando terra:
(dem = terra::rast("C:/Users/usuagro/Desktop/GB2R/CuadernoDEM/datos/elev_narino_z10.tif"))## class : SpatRaster
## size : 3581, 3604, 1 (nrow, ncol, nlyr)
## resolution : 76.36952, 76.36952 (x, y)
## extent : 4320038, 4595274, 1596968, 1870447 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +units=m +no_defs
## source : elev_narino_z10.tif
## name : elev_narino_z10
## min value : -2298
## max value : 4860Reduzcamos la resolución del MDE para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Ahora, leeremos nuestros municipios utilizando la librería
sf:
(munic <- sf::st_read("datos/MUNICIPIOS_9377.gpkg"))## Reading layer `Municipios' from data source
## `C:\Users\usuagro\Desktop\GB2R\CuadernoDEM\datos\MUNICIPIOS_9377.gpkg'
## using driver `GPKG'
## Simple feature collection with 64 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4330516 ymin: 1598084 xmax: 4573425 ymax: 1855793
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-NacionalAhora, recortemos el objeto dem2 a los límites de
nuestro departamento:
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Trasformar coordenadas (10%)
Para calcular los atributos del terreno, el MDE debe estar en coordenadas planas. Como recordarán, el IGAC cambió en 2020 a un sistema de coordenadas planas con origen único nacional. El código EPSG correspondiente al sistema de coordenadas planas de Origen Nacional de Colombia es 9377.
Modifique los siguientes dos bloques de código para transformar los datos de entrada al sistema de coordenadas nacional de Colombia.
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1687, 1590, 1 (nrow, ncol, nlyr)
## resolution : 152.739, 152.739 (x, y)
## extent : 4330577, 4573432, 1598114, 1855785 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_narino_z10
## min value : -220.5
## max value : 4724.5También necesitamos realizar una transformación similar para el objeto vectorial:
(munic_plane = sf::st_transform(munic, "EPSG:9377"))6.Computing terrain attributes (20%)
Busquemos cómo calcular la pendiente (slope) y la orientación (aspect).
Al escribir “?SlpAsp” nos dice efectivamente “Multiscale Slope and Aspect Description Calculates multiscale slope and aspect based on a modified version of the algorithm from Misiuk et al (2021) which extends classical formulations of slope restricted to a 3x3 window to multiple scales by using only cells on the edges of the focal window (see details for more information).”.
Ahora, llamemos a la función SlpAsp:
(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 : 1687, 1590, 2 (nrow, ncol, nlyr)
## resolution : 152.739, 152.739 (x, y)
## extent : 4330577, 4573432, 1598114, 1855785 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0, 0
## max values : 63.7891, 360Partimos el objeto slp_asp en 2 partes
Layer 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1687, 1590, 1 (nrow, ncol, nlyr)
## resolution : 152.739, 152.739 (x, y)
## extent : 4330577, 4573432, 1598114, 1855785 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 63.7891Descripción del resultado y su significado: El histograma muestra que la mayor parte del territorio es de pendientes bajas, entre 0° y 10°, con una frecuencia que disminuye a medida que aumentan los grados de inclinación. Esto revela un predominio de zonas planas a onduladas, con presencia menor de zonas montañosas con pendientes pronunciadas de maximo 60°.
terra::hist(slope,
main = "Narino's slope",
xlab = "Slope (in degrees)")
Layer 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1687, 1590, 1 (nrow, ncol, nlyr)
## resolution : 152.739, 152.739 (x, y)
## extent : 4330577, 4573432, 1598114, 1855785 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 360terra::hist(aspect,
main = "Narino's aspect",
xlab = "Aspect (in degrees)")
los grados de pendiente a porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1687, 1590, 1 (nrow, ncol, nlyr)
## resolution : 152.739, 152.739 (x, y)
## extent : 4330577, 4573432, 1598114, 1855785 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 203.12928terra::hist(slope_perc,
main = "Narino's slope",
xlab = "Slope (in percentage)")
terra::hist(slope_perc)
7. Computing zonal statistics (20%)
Una operación de estadística zonal es aquella que calcula estadísticas sobre los valores de celda de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En nuestro caso, estamos interesados en calcular el promedio de los valores de pendiente por municipio.
Primero, conozcamos la clasificación de pendientes del IGAC para fines agrológicos:
Ahora, pongamos manos a la obra:
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'))## | | | 0% | |= | 2% | |== | 3% | |=== | 5% | |==== | 6% | |===== | 8% | |======= | 9% | |======== | 11% | |========= | 12% | |========== | 14% | |=========== | 16% | |============ | 17% | |============= | 19% | |============== | 20% | |=============== | 22% | |================ | 23% | |================== | 25% | |=================== | 27% | |==================== | 28% | |===================== | 30% | |====================== | 31% | |======================= | 33% | |======================== | 34% | |========================= | 36% | |========================== | 38% | |=========================== | 39% | |============================ | 41% | |============================== | 42% | |=============================== | 44% | |================================ | 45% | |================================= | 47% | |================================== | 48% | |=================================== | 50% | |==================================== | 52% | |===================================== | 53% | |====================================== | 55% | |======================================= | 56% | |======================================== | 58% | |========================================== | 59% | |=========================================== | 61% | |============================================ | 62% | |============================================= | 64% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 69% | |================================================= | 70% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 75% | |====================================================== | 77% | |======================================================= | 78% | |======================================================== | 80% | |========================================================= | 81% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 86% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 91% | |================================================================= | 92% | |================================================================== | 94% | |=================================================================== | 95% | |==================================================================== | 97% | |===================================================================== | 98% | |======================================================================| 100%## [1] 21.080587 33.400925 9.423068 42.493084 33.332905 16.397694 28.878056
## [8] 28.337721 32.619019 39.599083 26.321308 33.884186 10.676829 30.390375
## [15] 38.427006 32.162880 19.320099 38.239620 37.415833 34.803310 36.027485
## [22] 37.178776 12.075477 29.880524 18.173370 24.162079 33.497238 31.407728
## [29] 30.082270 31.452751 32.958107 5.231098 23.030230 39.376263 46.487659
## [36] 30.090084 7.248667 32.195621 5.717163 24.430866 3.339289 23.524950
## [43] 4.156869 38.465538 29.589329 34.138847 37.771858 12.983334 28.781101
## [50] 3.360935 31.774059 33.629314 31.390396 32.647106 32.124718 26.845453
## [57] 9.088976 33.606876 11.302633 29.685736 27.069786 3.189805 22.982910
## [64] 28.912632hist(munic$mean_slope,
main = "Narino's mean slope",
xlab = "Slope (in percentage)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |= | 2% | |== | 3% | |=== | 5% | |==== | 6% | |===== | 8% | |======= | 9% | |======== | 11% | |========= | 12% | |========== | 14% | |=========== | 16% | |============ | 17% | |============= | 19% | |============== | 20% | |=============== | 22% | |================ | 23% | |================== | 25% | |=================== | 27% | |==================== | 28% | |===================== | 30% | |====================== | 31% | |======================= | 33% | |======================== | 34% | |========================= | 36% | |========================== | 38% | |=========================== | 39% | |============================ | 41% | |============================== | 42% | |=============================== | 44% | |================================ | 45% | |================================= | 47% | |================================== | 48% | |=================================== | 50% | |==================================== | 52% | |===================================== | 53% | |====================================== | 55% | |======================================= | 56% | |======================================== | 58% | |========================================== | 59% | |=========================================== | 61% | |============================================ | 62% | |============================================= | 64% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 69% | |================================================= | 70% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 75% | |====================================================== | 77% | |======================================================= | 78% | |======================================================== | 80% | |========================================================= | 81% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 86% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 91% | |================================================================= | 92% | |================================================================== | 94% | |=================================================================== | 95% | |==================================================================== | 97% | |===================================================================== | 98% | |======================================================================| 100%## [1] 4 5 3 5 5 4 5 5 5 5 4 5 4 5 5 5 1 5 5 5 5 5 2 5 4 4 4 5 5 5 5 1 4 5 6 5 1 5
## [39] 1 4 1 4 1 5 4 5 5 4 5 1 5 5 5 5 5 5 1 5 2 5 5 1 4 4hist(munic$class,
main = "Narino's reclassified slope",
xlab = "Slope (as a category)")
terra::hist(slope_perc)
Transdormamos las coordenadas:
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1695, 1589, 1 (nrow, ncol, nlyr)
## resolution : 0.00137313, 0.00137313 (x, y)
## extent : -79.01379, -76.83188, 0.3612346, 2.68869 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 174.919952(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1695, 1589, 1 (nrow, ncol, nlyr)
## resolution : 0.00137313, 0.00137313 (x, y)
## extent : -79.01379, -76.83188, 0.3612346, 2.68869 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 189.1986858. Plotting terrain slope (20%)
Primero, necesitaremos una paleta de colores para la pendiente (slope):
Ahora graficámos:
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Hacemos los cambios necesarios:
slope_leaflet <- terra::project(slope_perc, "EPSG:4326")munic_wgs84 <- sf::st_transform(munic, "EPSG:4326")slope_leaflet_opt <- terra::aggregate(slope_leaflet, fact = 2, fun = mean)
leaflet(munic_wgs84) %>%
addTiles() %>%
setView(-77.3, 1.2, 9) %>%
addRasterImage(slope_leaflet, colors = palredgreen, opacity = 0.8) %>%
addPolygons(
color = "black",
weight = 1.5,
smoothFactor = 0.5,
opacity = 1.0,
fillOpacity = 0.0,
popup = paste("Municipio: ", munic_wgs84$mpio_cnmbr, "<br>",
"slope class: ", munic_wgs84$class, "<br>")
) %>%
addLegend(pal = palredgreen, values = values(slope_leaflet),
title = "Terrain slope in Narino (%)")