Geomorphometry of Magdalena in R
Santiago Suarez
2026-07-22
1. Introducción.
En esta práctica se realizó el análisis de algunos atributos del terreno, como la pendiente y la orientación, utilizando un Modelo Digital de Elevación (DEM). Para ello se trabajó con información del departamento del Magdalena y se empleó el software R, junto con el paquete MultiscaleDTM, con el fin de obtener una mejor descripción de las características topográficas de la zona.
2. Configuración inicial.
Primero limpiamos la memoria de R para asegurar que no tengamos variables antiguas rondando en la sesión.
rm(list=ls())Luego cargamos las librerías necesarias para la lectura de datos espaciales, cálculos de pendientes, estadísticas por municipio y mapas interactivos.
library(elevatr)## Warning: package 'elevatr' was built under R version 4.5.3## 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)## Warning: package 'sf' was built under R version 4.5.3## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUElibrary(leaflet)## Warning: package 'leaflet' was built under R version 4.5.3library(terra)## Warning: package 'terra' was built under R version 4.5.3## terra 1.9.27library(MultiscaleDTM)## Warning: package 'MultiscaleDTM' was built under R version 4.5.3library(exactextractr)## Warning: package 'exactextractr' was built under R version 4.5.33. Sobre la librería MultiscaleDTM.
Para procesar el relieve usamos MultiscaleDTM, un paquete diseñado para extraer variables topográficas mediante el análisis de celdas vecinas (ventanas focales). En particular, la función SlpAsp evalúa las diferencias de altura en matrices de píxeles (por ejemplo de \(3 \times 3\)) para determinar tanto qué tan inclinado está el terreno como hacia qué dirección mira la pendiente.
4. Carga y preparación de datos.
Comenzamos leyendo el Modelo Digital de Elevación (DEM) en formato ráster que preparamos previamente para el Magdalena.
(dem = terra::rast("C://Users//suare//OneDrive//Documentos//GB2//Proyecto1//Datos//elev_Magdalena_z10.tif"))## class : SpatRaster
## size : 4078, 2589, 1 (nrow, ncol, nlyr)
## resolution : 0.0006789017, 0.0006789017 (x, y)
## extent : -75.23438, -73.4767, 8.754526, 11.52309 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev_Magdalena_z10.tif
## name : elev_Magdalena_z10
## min value : -2439
## max value : 5676Para acelerar el procesamiento y optimizar el rendimiento del computador, agrupamos los píxeles de 2 en 2 promediando sus alturas.
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Leemos los municipios del departamento con sf.
munic <- sf::st_read("C://Users//suare//OneDrive//Documentos//GB2//Proyecto1//Datos//Municipios_Magdalena.gpkg",layer = "munic")## Reading layer `munic' from data source
## `C:\Users\suare\OneDrive\Documentos\GB2\Proyecto1\Datos\Municipios_Magdalena.gpkg'
## using driver `GPKG'
## Simple feature collection with 30 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.9466 ymin: 8.936489 xmax: -73.54184 ymax: 11.34891
## Geodetic CRS: WGS 84Recortamos para que coincida exactamente con los limites del Magdalena.
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Transformación de coordenadas a sistema plano.
Para poder calcular pendientes e inclinaciones reales, el DEM debe estar en metros y no en grados geográficos. En Colombia utilizamos el sistema oficial de Origen Nacional (EPSG:9377).
Convertimos tanto el ráster recortado como la capa de municipios a este sistema plano:
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_Magdalena_z10
## min value : -158.965286
## max value : 5653.640625(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 30 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4786799 ymin: 2545597 xmax: 4940782 ymax: 2812311
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 47 001 47001 MAGDALENA SANTA MARTA
## 2 47 030 47030 MAGDALENA ALGARROBO
## 3 47 053 47053 MAGDALENA ARACATACA
## 4 47 058 47058 MAGDALENA ARIGUANÍ
## 5 47 161 47161 MAGDALENA CERRO DE SAN ANTONIO
## 6 47 170 47170 MAGDALENA CHIVOLO
## 7 47 189 47189 MAGDALENA CIÉNAGA
## 8 47 205 47205 MAGDALENA CONCORDIA
## 9 47 245 47245 MAGDALENA EL BANCO
## 10 47 258 47258 MAGDALENA EL PIÑÓN
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1525 MUNICIPIO 2343.9033 2025
## 2 Ordenanza 008 de Junio 24 de 1999 MUNICIPIO 405.8850 2024
## 3 Ordenanza 47 del 28 de Abril de 1915 MUNICIPIO 1738.4243 2025
## 4 Ordenanza 14 Bis de Noviembre 30 de 1966 MUNICIPIO 1130.3299 2024
## 5 Ordenanza 3038 de 1912 MUNICIPIO 177.0656 2024
## 6 Decreto 107 de Marzo 8 de 1974 MUNICIPIO 536.0309 2024
## 7 Ley 0339 de 1876 MUNICIPIO 1329.3560 2025
## 8 Ordenanza 007 del 24 de Junio de 1999 MUNICIPIO 109.3953 2024
## 9 Ley 0182 de 1871 MUNICIPIO 813.8018 2025
## 10 Ordenanza 032 del 20 de Abril de 1815 MUNICIPIO 556.6568 2024
## shape_Leng shape_Area AREA geom
## 1 3.1304030 0.194233447 2347.0954 MULTIPOLYGON (((4867382 278...
## 2 1.2656799 0.033536941 406.3880 MULTIPOLYGON (((4881075 267...
## 3 2.6435461 0.143834703 1740.7962 MULTIPOLYGON (((4934596 275...
## 4 1.9928810 0.093278409 1131.7403 MULTIPOLYGON (((4886737 266...
## 5 0.7004114 0.014623054 177.1748 MULTIPOLYGON (((4797435 270...
## 6 1.3253970 0.044254975 536.5174 MULTIPOLYGON (((4834685 266...
## 7 2.6311586 0.110066127 1331.0537 MULTIPOLYGON (((4861849 277...
## 8 0.5443263 0.009033178 109.4677 MULTIPOLYGON (((4801643 269...
## 9 1.6767358 0.067024887 814.8694 MULTIPOLYGON (((4912230 255...
## 10 1.1771096 0.045985558 557.0816 MULTIPOLYGON (((4816767 271...6. Cálculo de atributos del terreno: Pendiente y Orientación.
Ejecutamos la función SlpAsp. Configuramos el parámetro w = c(3, 3) para indicarle a R que analice ventanas focales de 3x3 celdas, y usamos method = “queen” para considerar los 8 píxeles vecinos (lados y diagonales) en el cálculo.
(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 : 1785, 1034, 2 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (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 : 61.324694, 359.999796Separamos el resultado en dos capas individuales: # Capa 1: Pendiente (grados).
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 61.324694Visualizamos la distribución de las pendientes en grados:
terra::hist(slope,
main = "Magdalena's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
Análisis del histograma:El gráfico muestra que la gran mayoría de las
celdas tienen pendientes muy bajas, casi todas concentradas entre 0° y
5° (con más de 400.000 celdas). Esto tiene todo el sentido con la
geografía del Magdalena, ya que refleja las zonas planas y bajas como el
valle del río Magdalena y la Ciénaga Grande.Hacia la derecha se ve una
cola más larga pero con frecuencias mucho más bajas, que va disminuyendo
hasta llegar cerca de los 40° o 50 °. Esta parte representa las zonas
con pendientes inclinadas y terrenos escarpados de la Sierra Nevada de
Santa Marta. # Capa 2: Orientación / Aspecto (grados azimutales).
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 359.999796Visualizamos la distribución de la orientación de las laderas:
terra::hist(aspect,
main = "Magdalena's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 54% of the cells was used (of which 45% was NA)
El histograma confirma que las montañas del Magdalena miran hacia todos
los puntos cardinales. Sin embargo, no es una distribución perfectamente
plana. Las barras muestran un “valle” en el sector del Noreste (50°-
75°), indicando que es la orientación menos común. A partir de ahí, la
frecuencia aumenta gradualmente hasta alcanzar su punto máximo en el
Noroeste y Oeste (275°- 300°), donde se concentra la mayor cantidad de
laderas. # Conversión de grados a porcentaje. Transformamos la pendiente
a porcentaje usando la relación trigonométrica tan(grados)x100:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 182.840775Revisamos su histograma en porcentaje:
terra::hist(slope_perc,
main = "Magdalena's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
terra::hist(slope_perc)## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
7. Estadísticas zonales por municipio.
Queremos resumir qué tipo de pendiente predomina en cada municipio del departamento. Para esto, reclasificamos la pendiente en porcentaje siguiendo los rangos agrológicos del IGAC.
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)Calculamos la pendiente promedio en porcentaje para cada municipio 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% | |== | 3% | |===== | 7% | |======= | 10% | |========= | 13% | |============ | 17% | |============== | 20% | |================ | 23% | |=================== | 27% | |===================== | 30% | |======================= | 33% | |========================== | 37% | |============================ | 40% | |============================== | 43% | |================================= | 47% | |=================================== | 50% | |===================================== | 53% | |======================================== | 57% | |========================================== | 60% | |============================================ | 63% | |=============================================== | 67% | |================================================= | 70% | |=================================================== | 73% | |====================================================== | 77% | |======================================================== | 80% | |========================================================== | 83% | |============================================================= | 87% | |=============================================================== | 90% | |================================================================= | 93% | |==================================================================== | 97% | |======================================================================| 100%## [1] 34.3125420 0.9711066 32.4357719 1.3878907 0.9490672 1.5636197
## [7] 34.2732773 1.8798674 0.9540940 0.7222221 0.5690165 21.9244385
## [13] 0.9261801 1.3292797 2.5145578 1.1275929 1.0748733 1.7687278
## [19] 0.4968248 0.5313346 1.6415391 0.4778264 1.2094046 0.9828508
## [25] 1.0238775 1.5891575 0.5086805 2.3828094 1.4427539 1.9263012Graficamos la pendiente media:
hist(munic$mean_slope,
main = "Magdalena's mean slope",
xlab = "Slope (in percentage)")
Análisis del histograma: La gran mayoría de las unidades analizadas (26
de ellas) presentan una pendiente media muy baja, entre el 0% y el 5%,
lo que reafirma el carácter predominantemente plano del departamento.
Solo una minoría registra promedios más altos, con un caso aislado entre
20% y 25%, y 3 unidades entre 30% y 35%, correspondientes a las zonas de
mayor relieve montañoso.
Y obtenemos la clase de pendiente más frecuente (moda) dentro de cada territorio:
(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% | |== | 3% | |===== | 7% | |======= | 10% | |========= | 13% | |============ | 17% | |============== | 20% | |================ | 23% | |=================== | 27% | |===================== | 30% | |======================= | 33% | |========================== | 37% | |============================ | 40% | |============================== | 43% | |================================= | 47% | |=================================== | 50% | |===================================== | 53% | |======================================== | 57% | |========================================== | 60% | |============================================ | 63% | |=============================================== | 67% | |================================================= | 70% | |=================================================== | 73% | |====================================================== | 77% | |======================================================== | 80% | |========================================================== | 83% | |============================================================= | 87% | |=============================================================== | 90% | |================================================================= | 93% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 1 5 1 1 1 5 1 1 1 1 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1Graficamos la distribución:
hist(munic$class,
main = "Magdalena's reclassified slope",
xlab = "Slope (as a category)")
Finalmente, reproyectamos las capas de nuevo a coordenadas geográficas
(WGS84) para poder cargarlas en los mapas web interactivos:
terra::hist(slope_perc)## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1785, 1048, 1 (nrow, ncol, nlyr)
## resolution : 0.001357862, 0.001357862 (x, y)
## extent : -74.96115, -73.53811, 8.931431, 11.35521 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 167.753052(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1785, 1048, 1 (nrow, ncol, nlyr)
## resolution : 0.001357862, 0.001357862 (x, y)
## extent : -74.96115, -73.53811, 8.931431, 11.35521 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 167.7530528. Visualización de la pendiente en porcentaje.
Definimos una rampa de color continua para resaltar la inclinación (de verde claro para zonas planas a rojo oscuro para zonas de alta pendiente):
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Creamos el primer mapa interactivo en Leaflet:
leaflet(munic) %>% addTiles() %>% setView(-74.1, 10.4, 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 = "Terrain slope in Magdalena (%)")