Atributos Geomorfométricos del terreno, departamento de Boyacá
Yaneth Alexandra Arenas Garrido
2025-06-12
1. Introducción
En este cuaderno se muestra el resultado del cálculo de los atributos geomorfométricos del terreno, a partir de modelos digitales de elevación en cuadrícula y haciendo uso del paquete MultiscaleDTM en R.
El objetivo es evidenciar el conocimiento adquirido sobre modelos de elevación, generando un mapa que presente las características geomorfométricas del departamento de Boyacá, Colombia.
2. Configuración
Para dar inicio, realizamos la instalación de los paquetes desde la consola y limpiamos la memoria del sitio de trabajo:
rm(list=ls())Luego, cargamos las librerias necesarias para el desarrollo del código:
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)3. Leer salida de datos
En primer lugar, cargaremos un modelo de elevación (DEM) ráster para el departamento de Boyacá. Después leeremos el DEM con ayuda del paquete Terra:
(dem = terra::rast("boyaca/elev_boyaca_z10.tif"))## class : SpatRaster
## size : 4083, 4619, 1 (nrow, ncol, nlyr)
## resolution : 0.0006849771, 0.0006849771 (x, y)
## extent : -74.88281, -71.7189, 4.565705, 7.362467 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_boyaca_z10.tif
## name : elev_boyaca_z10
## min value : -272
## max value : 5335Ahora, vamos a disminuir la resolución del DEM para evitar ocupar mucho espacio en la memoria y facilitar el procesamiento de datos:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Seguimos con la lectura de los municipios del departamento usando la librería sf:
(munic <- sf::st_read("boyaca/Boyaca.gpkg"))## Reading layer `boyaca' from data source
## `C:\Users\arena\Desktop\cuaderno_geomorfometria\boyaca\Boyaca.gpkg'
## using driver `GPKG'
## Simple feature collection with 123 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: MAGNA-SIRGASHay que notar que el modelo de elevación digital (dem2) como los límites del departamento (munic) están en el mismo sistema de coordenadas geográficas. Vamos a recortar dem2 para obtener los valores solo dentro de los límites del departamento (munic):
dem3 = terra::crop(dem2,munic, mask=TRUE)## Warning: [crop] CRS do not match4. Transformar coordenadas
Para calcular los atributos del terreno, debemos transformar el DEM en coordenadas planas, ya que así lo indica el sistema de coordenadas nacional de Colombia. Para esto vamos a usar el código EPSG:9377, la función project () del paquete terra y guardamos la información en el objeto dem_plane:
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1753, 1988, 1 (nrow, ncol, nlyr)
## resolution : 151.5067, 151.5067 (x, y)
## extent : 4815384, 5116579, 2072369, 2337961 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_boyaca_z10
## min value : 96.33419
## max value : 5295.17871De igual manera, transformamos el objeto vectorial munic al mismo sistema de coordenadas planas usando la función st_transform() del paquete sf y guardamos el resultado en el objeto munic_plane:
(munic_plane = sf::st_transform(munic, "EPSG:9377"))5. Calcular atributos del terreno
Vamos a calcular la pendiente y la orientación del terreno a partir del modelo de elevación dem_plane. Para esto usamos la función SlpAsp() del paquete MultiscaleDTM. Comprendamos algunos parámetros del código:
w: es un vector de dos elementos que indica el tamaño de la ventana rectangular a utilizar, donde el primer valor representa la cantidad de filas y el segundo la cantidad de columnas. Ambos valores deben ser números impares.
method: el método indica qué celdas usar en los cálculos. En este caso, “queen” agrega cuatro celdas de esquina adicionales, es decir, considera las 8 celdas vecinas alrededor de cada celda central. Este método permite incorporar más información del entorno.
El resultado es un raster con dos capas: una para la pendiente y otra para la orientación de la pendiente:
(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 : 1753, 1988, 2 (nrow, ncol, nlyr)
## resolution : 151.5067, 151.5067 (x, y)
## extent : 4815384, 5116579, 2072369, 2337961 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 3.389264e-04
## max values : 67.50974, 3.599999e+02Ahora, vamos a dividir el objeto slp_asp en dos partes:
Primera parte, pendiente del terreno:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1753, 1988, 1 (nrow, ncol, nlyr)
## resolution : 151.5067, 151.5067 (x, y)
## extent : 4815384, 5116579, 2072369, 2337961 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 67.50974Luego, realizamos un histograma con la función hist del paquete terra. En este gráfico se observa la distribución del rango de pendientes en el departamento:
terra::hist(slope,
main = "Pendiente de Boyacá",
xlab = "Pendiente (en grados)",
ylab = "Frecuencia")## Warning: [hist] a sample of 29% of the cells was used (of which 71% was NA)
En el histograma podemos observar que el terreno presenta pendientes bajas, siendo las pendientes de 5 a 20 grados las más frecuentes. Esto sugiere que Boyacá tiene un relieve ondulado a moderadamente inclinado en gran parte de su territorio.
Segunda parte, orientación del terreno:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1753, 1988, 1 (nrow, ncol, nlyr)
## resolution : 151.5067, 151.5067 (x, y)
## extent : 4815384, 5116579, 2072369, 2337961 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 3.389264e-04
## max value : 3.599999e+02Después, al igual que en la primera parte, creamos un histograma para observar la distribución de los valores de la oreientación en el deparatmento:
terra::hist(aspect,
main = "Orientación de Boyacá",
xlab = "Orientación (en grados)",
ylab = "Frecuencia")## Warning: [hist] a sample of 29% of the cells was used (of which 71% was NA)
El histograma muestra que hay un leve aumento en orientaciones entre 90° y 150°, y otra entre 270° y 330°. Esto podría indicar que el relieve de Boyacá tiene una cierta preferencia de orientación probablemente por las formaciones montañosas, aunque cabe resaltar que la orientación varía en todo el territorio.
Como paso sigueinte, convertimos los valores de la pendiente de grados a porcentaje:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1753, 1988, 1 (nrow, ncol, nlyr)
## resolution : 151.5067, 151.5067 (x, y)
## extent : 4815384, 5116579, 2072369, 2337961 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 241.5375Ahora, representamos los valores de la pendiente en porcentaje en un histograma. En este caso el objeto de entreda es slope_perc:
terra::hist(slope_perc,
main = "Pendiente de Boyacá",
xlab = "Pendiente (en porcentaje)",
ylab = "Frecuencia")## Warning: [hist] a sample of 29% of the cells was used (of which 71% was NA)
Podemos observar que la mayoría del terreno tiene pendientes entre 0% y 40%. A medida que aumenta la pendiente, la frecuencia disminuye rápidamente, es decir, las Zonas con alta pendiente existen, pero son pocas.
6. Cálculo de estadísticas zonales
En primer lugar vamos a reclasificar los valores de la pendiente de acuerdo a su rango en porcentajes:
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
75, 250, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)Luego, con el paquete exactextractr, Para cada polígono en el shapefile munic, se calcula la pendiente media usando el raster slope_perc y el resultado se guarda como nueva columna mean_slope en el objeto munic:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 12.518128 21.144888 23.988548 22.717548 26.608459 22.994267 26.524553
## [8] 33.178539 22.019363 24.242296 23.947521 20.284788 13.821063 33.908302
## [15] 25.449083 27.301411 13.390441 35.120361 29.212452 21.158005 12.456381
## [22] 16.181677 16.712172 29.751104 22.443031 39.557938 32.825733 15.269210
## [29] 15.684027 25.852493 38.720631 22.886679 24.836580 30.142498 12.261699
## [36] 22.412165 16.538471 25.376986 23.486912 32.943840 24.244211 26.724607
## [43] 41.656013 16.976667 18.397699 29.534027 34.195229 26.950207 32.220211
## [50] 32.070927 17.376850 31.695009 28.869654 23.904516 29.410658 19.923866
## [57] 14.821651 12.593950 34.014755 20.795874 20.460703 7.433327 29.419920
## [64] 27.215298 33.505062 21.688156 32.729694 37.937840 30.436802 31.650448
## [71] 29.407154 21.642443 34.988293 6.843872 30.364162 19.803070 16.886784
## [78] 28.387653 12.643967 17.540489 15.273659 21.158850 14.760439 21.637327
## [85] 40.014954 7.708375 33.039330 14.597834 36.630981 22.220503 18.159672
## [92] 30.547890 28.861370 11.565438 30.643137 28.562035 27.823042 17.066616
## [99] 26.394609 17.215567 23.752096 11.932631 30.993214 15.616819 25.768589
## [106] 25.804178 20.701824 21.864691 13.512199 15.909994 33.270813 11.754683
## [113] 18.250748 22.767176 17.036268 22.818758 22.437641 10.370167 26.109570
## [120] 23.181223 15.892552 17.797779 28.281349Como se observa, el CRS del objeto munic fue transformado automáticamente al CRS del ráster, es decir, a EPSG:9377.
Como paso siguiente, realizamos un histograma para mostrar la distribución de los valores de la pendiente media:
hist(munic$mean_slope,
main = "Pendiente media de Boyacá",
xlab = "Pendiente (en porcentaje)",
ylab = "Frecuencia")
En el histograma vemos que la mayoría de los valores se concentran entre 15% y 30% de pendiente media. La distribución relativamente es simétrica, con una leve inclinación por valores menores de pendiente media pero en menor proporción. De igual manera, presenta pendientes superiores al 35% aunque en menor cantidad.
Después, con el paquete exactextractr creamos una nueva columna llamada class en el objeto munic, que contiene la moda (el valor más frecuente) de los valores del raster rc dentro de cada polígono del shapefile munic:
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |========= | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 4 4 4 4 5 4 5 5 4 4 4 4 4 5 5 5 4 5 5 4 4 4 4 5 4 5 5 4 4 4 5 5 4 5 4 5 4
## [38] 4 4 5 4 5 5 4 4 5 5 4 5 5 3 5 5 4 5 4 4 4 5 5 4 2 5 4 5 4 5 5 5 5 5 4 5 1
## [75] 5 4 4 5 4 4 4 4 4 4 5 1 5 4 5 4 4 5 5 4 5 5 5 4 5 4 4 4 5 4 4 5 4 4 4 4 5
## [112] 4 4 4 4 4 4 3 5 4 4 4 5Representamos el resultado en un histograma:
hist(munic$class,
main = "Pendiente reclasificada de Boyacá",
xlab = "Pendiente (como categoría)",
ylab = "Frecuencia")
Se observa que la pendiente más frecuente es la que corresponde a la categoría de 4.
Ahora transformamos el sistema de coordenadas de del ráster a coordenadas geográficas, es decir, EPSG:4326:
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1756, 1991, 1 (nrow, ncol, nlyr)
## resolution : 0.00136995, 0.00136995 (x, y)
## extent : -74.67211, -71.94454, 4.652932, 7.058565 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7Realizamos nuevamente la transformación del sistema de coordenadas, ahora para slope_perc:
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1756, 1991, 1 (nrow, ncol, nlyr)
## resolution : 0.00136995, 0.00136995 (x, y)
## extent : -74.67211, -71.94454, 4.652932, 7.058565 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 228.81437. Gráficar la pendiente del terreno
El objetivo es visualizar la pendiente (en porcentaje) del terreno sobre el deparatamento de Boyacá. Para esto, es necesario crear la paleta de colores y un mapa inetractivo:
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ahora creamos un mapa interactivo haciendo uso del paquete leaflet, que muestre la pendiente del terreno (en porcentaje) en Boyacá, superpuesta sobre los límites municipales, con una leyenda y ventanas emergentes informativas de cada lugar:
leaflet(munic) %>% addTiles() %>% setView(-72.9, 6.0, 8) %>%
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 Boyacá (%)")## 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'