Cálculo de variables geomorfométricas en el departamento de La Guajira
Julian Santiago Gonzalez Malpica
04.11.2025
1. Introducción
Este cuaderno muestra cómo calcular atributos geomorfométricos del terreno a partir de modelos digitales de elevación en formato de cuadrícula (DEM) para el departamento de La Guajira. Para el cálculo de los atributos del terreno, el paquete principal que utilizaremos es MultiscaleDTM en R.
2. Configuración
Limpiamos la memoria:
rm(list=ls())Cargamos las librerías necesarias:
library(sp)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)Librerías utilizadas:
- terra: lectura, procesamiento y análisis de datos espaciales ráster
- sf: trabajar con datos espaciales vectoriales
- leaflet: creación de mapas interactivos
- MultiscaleDTM: análisis de terreno a diferentes escalas espaciales
- exactextractr: extracción precisa de valores ráster dentro de polígonos
3. Introducción a MultiscaleDTM
Este paquete permite calcular atributos geomorfométricos del terreno
a múltiples escalas. La función SlpAsp() calcula la
pendiente (slope) y la exposición
(aspect) usando modificaciones de algoritmos tradicionales.
Métodos de análisis espacial:
- Método Rook: usa 4 celdas ortogonales (N, S, E, O)
- Método Queen: usa 8 celdas vecinas (incluye diagonales)
- Método Boundary: usa todas las celdas del borde exterior
4. Leer datos de entrada
Verificamos los archivos disponibles:
list.files()## [1] "COL_msk_alt.tif" "geomorfoGuajira.html" "geomorfoGuajira.Rmd"
## [4] "geomorfoGuajira_files" "MunicGuajira.gpkg" "rsconnect"Cargamos el DEM de Colombia:
dem <- terra::rast("COL_msk_alt.tif")
dem## class : SpatRaster
## size : 2136, 1800, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -81.8, -66.8, -4.3, 13.5 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source : COL_msk_alt.tif
## name : COL_msk_alt
## min value : -28
## max value : 5453Reducimos la resolución del DEM para optimizar memoria:
dem2 <- terra::aggregate(dem, 2, "mean")## |---------|---------|---------|---------|========================================= dem2## class : SpatRaster
## size : 1068, 900, 1 (nrow, ncol, nlyr)
## resolution : 0.01666667, 0.01666667 (x, y)
## extent : -81.8, -66.8, -4.3, 13.5 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : COL_msk_alt
## min value : -3.00
## max value : 5222.75Cargamos los municipios de La Guajira:
munic <- sf::st_read("MunicGuajira.gpkg")## Reading layer `municipios' from data source
## `C:\Users\jsmal\OneDrive\Documents\GEOMATICA\GB2\Proyecto 4\MunicGuajira.gpkg'
## using driver `GPKG'
## Simple feature collection with 10 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.6846 ymin: 10.3812 xmax: -70.99986 ymax: 12.45847
## Geodetic CRS: WGS 84head(munic)## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.6846 ymin: 10.53404 xmax: -71.99616 ymax: 11.78986
## Geodetic CRS: WGS 84
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2 ENGTYPE_2
## 1 53 COL Colombia 18 La Guajira 653 Barrancas Municipio Municipality
## 2 53 COL Colombia 18 La Guajira 654 El Molino Municipio Municipality
## 3 53 COL Colombia 18 La Guajira 655 Fonseca Municipio Municipality
## 4 53 COL Colombia 18 La Guajira 656 Maicao Municipio Municipality
## 5 53 COL Colombia 18 La Guajira 657 Manaure Municipio Municipality
## 6 53 COL Colombia 18 La Guajira 658 Riohacha Municipio Municipality
## NL_NAME_2 VARNAME_2 geom
## 1 <NA> <NA> MULTIPOLYGON (((-72.67934 1...
## 2 <NA> <NA> MULTIPOLYGON (((-72.83597 1...
## 3 <NA> <NA> MULTIPOLYGON (((-72.69444 1...
## 4 <NA> <NA> MULTIPOLYGON (((-71.99616 1...
## 5 <NA> <NA> MULTIPOLYGON (((-72.92757 1...
## 6 <NA> <NA> MULTIPOLYGON (((-72.6719 11...5. Transformación de coordenadas
Consultamos el CRS de los municipios:
(Guajira_crs <- crs(munic))## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Consultamos el CRS del DEM:
crs(dem2)## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"Asignamos el CRS de La Guajira al DEM:
crs(dem2) <- Guajira_crs
crs(dem2)## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Recortamos el DEM a los límites de La Guajira:
dem3 <- terra::crop(dem2, munic, mask=TRUE)6. Reproyección a coordenadas planas
Transformamos el DEM a coordenadas planas (EPSG:9377 - MAGNA-SIRGAS 2018 / Origen-Nacional):
dem_plane <- project(dem3, "EPSG:9377")
dem_plane## class : SpatRaster
## size : 126, 161, 1 (nrow, ncol, nlyr)
## resolution : 1826.923, 1826.923 (x, y)
## extent : 4925227, 5219362, 2706282, 2936475 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : COL_msk_alt
## min value : -2.065217
## max value : 4866.989746Transformamos también el objeto vectorial:
munic_plane <- sf::st_transform(munic, "EPSG:9377")
head(munic_plane)## Simple feature collection with 6 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4925236 ymin: 2722047 xmax: 5109398 ymax: 2860911
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2 ENGTYPE_2
## 1 53 COL Colombia 18 La Guajira 653 Barrancas Municipio Municipality
## 2 53 COL Colombia 18 La Guajira 654 El Molino Municipio Municipality
## 3 53 COL Colombia 18 La Guajira 655 Fonseca Municipio Municipality
## 4 53 COL Colombia 18 La Guajira 656 Maicao Municipio Municipality
## 5 53 COL Colombia 18 La Guajira 657 Manaure Municipio Municipality
## 6 53 COL Colombia 18 La Guajira 658 Riohacha Municipio Municipality
## NL_NAME_2 VARNAME_2 geom
## 1 <NA> <NA> MULTIPOLYGON (((5035043 274...
## 2 <NA> <NA> MULTIPOLYGON (((5017940 272...
## 3 <NA> <NA> MULTIPOLYGON (((5033394 274...
## 4 <NA> <NA> MULTIPOLYGON (((5109398 283...
## 5 <NA> <NA> MULTIPOLYGON (((5007894 283...
## 6 <NA> <NA> MULTIPOLYGON (((5035816 278...7. Cálculo de atributos del terreno
Calculamos pendiente y exposición usando
MultiscaleDTM::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
)
slp_asp## class : SpatRaster
## size : 126, 161, 2 (nrow, ncol, nlyr)
## resolution : 1826.923, 1826.923 (x, y)
## extent : 4925227, 5219362, 2706282, 2936475 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.004260348, 7.965948e-03
## max values : 20.650416231, 3.599769e+02Separamos las capas:
Capa 1 - Pendiente:
slope <- subset(slp_asp, 1)
slope## class : SpatRaster
## size : 126, 161, 1 (nrow, ncol, nlyr)
## resolution : 1826.923, 1826.923 (x, y)
## extent : 4925227, 5219362, 2706282, 2936475 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.004260348
## max value : 20.650416231Histograma de pendientes:
terra::hist(slope,
main = "Pendientes de La Guajira",
xlab = "Pendiente (en grados)")
Capa 2 - Exposición:
aspect <- subset(slp_asp, 2)
aspect## class : SpatRaster
## size : 126, 161, 1 (nrow, ncol, nlyr)
## resolution : 1826.923, 1826.923 (x, y)
## extent : 4925227, 5219362, 2706282, 2936475 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 7.965948e-03
## max value : 3.599769e+02Histograma de exposición:
terra::hist(aspect,
main = "Exposición de La Guajira",
xlab = "Exposición (en grados)")
Convertimos pendiente de grados a porcentaje:
slope_perc <- tan(slope * (pi/180)) * 100
slope_perc## class : SpatRaster
## size : 126, 161, 1 (nrow, ncol, nlyr)
## resolution : 1826.923, 1826.923 (x, y)
## extent : 4925227, 5219362, 2706282, 2936475 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00743571
## max value : 37.68798690Histograma de pendientes en porcentaje:
terra::hist(slope_perc,
main = "Pendientes de La Guajira",
xlab = "Pendiente (en porcentaje)")
8. Cálculo de estadísticas zonales
Clasificación de pendientes según IGAC:
| Código | Clase de pendiente | Pendiente (%) |
|---|---|---|
| a (1) | Ligeramente plano | 0 - 3 |
| b (2) | Ligeramente inclinado | 3 - 7 |
| c (3) | Moderadamente inclinado | 7 - 12 |
| d (4) | Fuertemente inclinado | 12 - 25 |
| e (5) | Moderadamente escarpado | 25 - 50 |
| f (6) | Escarpado | 50 - 75 |
| g (7) | Muy escarpado | 75 - 160 |
Reclasificamos el ráster de pendientes:
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 por municipio:
munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean')## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%munic$mean_slope## [1] 5.1060228 13.7118788 4.6819715 1.0663840 0.2069623 7.0002050
## [7] 4.4011941 1.0225923 9.3844337 5.8934627Histograma de pendientes promedio por municipio:
hist(munic$mean_slope,
main = "La Guajira - Pendiente promedio por municipio",
xlab = "Pendiente (en grados)")
Clase de pendiente dominante por municipio:
munic$class <- exactextractr::exact_extract(rc, munic, 'mode')## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |======= | 10% | |============== | 20% | |===================== | 30% | |============================ | 40% | |=================================== | 50% | |========================================== | 60% | |================================================= | 70% | |======================================================== | 80% | |=============================================================== | 90% | |======================================================================| 100%munic$class## [1] 1 4 1 1 1 1 1 1 4 1Histograma de clases de pendiente:
hist(munic$class,
main = "Pendiente reclasificada de La Guajira",
xlab = "Pendiente (como categoría)")
9. Reproyección a coordenadas geográficas
Transformamos los rásters de vuelta a WGS84:
rc.geo <- project(rc, "EPSG:4326")
rc.geo## class : SpatRaster
## size : 125, 162, 1 (nrow, ncol, nlyr)
## resolution : 0.01666694, 0.01666694 (x, y)
## extent : -73.68836, -70.98831, 10.39071, 12.47407 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 5slope.geo <- project(slope_perc, "EPSG:4326")
slope.geo## class : SpatRaster
## size : 125, 162, 1 (nrow, ncol, nlyr)
## resolution : 0.01666694, 0.01666694 (x, y)
## extent : -73.68836, -70.98831, 10.39071, 12.47407 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.009176374
## max value : 35.33549499510. Mapa interactivo de pendientes
Creamos paleta de colores:
palredgreen <- colorNumeric(
c("darkseagreen3", "yellow2", "orange", "brown2", "darkred"),
values(slope.geo),
na.color = "transparent"
)Mapa interactivo con Leaflet:
# Coordenadas calculadas del extent de La Guajira:
# Centro X: -72.34 (promedio de -73.69 y -70.99)
# Centro Y: 11.43 (promedio de 10.39 y 12.47)
leaflet(munic) %>%
addTiles() %>%
setView(-72.34, 11.43, zoom = 9) %>% # Coordenadas centradas en La Guajira
addPolygons(
color = "gray",
weight = 1.0,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste(
"<b>Municipio:</b> ", munic$NAME_2, "<br>",
"<b>Clase de pendiente:</b> ", munic$class, "<br>",
"<b>Pendiente promedio:</b> ", round(munic$mean_slope, 2), "%<br>"
)
) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(
pal = palredgreen,
values = values(slope.geo),
title = "Pendientes en La Guajira (%)",
position = "bottomright"
)11. Visualización estática de pendientes
Recorte del DEM original para La Guajira:
# Reproyectamos munic de vuelta a WGS84 para que coincida con dem
munic_wgs84 <- sf::st_transform(munic, "EPSG:4326")
# Ahora sí hacemos el recorte sin warnings
dem_guajira <- terra::mask(dem, munic_wgs84)## Warning: [mask] CRS do not matchdem_guajira <- terra::crop(dem_guajira, munic_wgs84)Cálculo de pendiente en grados:
slope_angle_guajira <- terra::terrain(
dem_guajira,
v = "slope",
unit = "degrees"
)Mapa estático de pendientes:
rgb.palette <- colorRampPalette(c("blue", "green", "yellow", "orange", "red"))
spplot(
slope_angle_guajira,
main = "Pendiente del terreno en grados - Departamento de La Guajira",
col.regions = rgb.palette(100),
sp.layout = list(
list("sp.polygons", as_Spatial(munic), col = "black", lwd = 1.2)
),
xlim = c(ext(dem_guajira)[1], ext(dem_guajira)[2]),
ylim = c(ext(dem_guajira)[3], ext(dem_guajira)[4])
)
12. Visualización de la orientación del terreno
Cálculo de aspect (exposición):
aspect_guajira <- terra::terrain(
dem_guajira,
v = "aspect",
unit = "degrees"
)Mapa estático de exposición:
rgb.palette <- colorRampPalette(
c("blue", "cyan", "green", "yellow", "orange", "red", "magenta", "blue")
)
spplot(
aspect_guajira,
main = "Exposición del terreno en grados - Departamento de La Guajira",
col.regions = rgb.palette(100),
sp.layout = list(
list("sp.polygons", as_Spatial(munic), col = "black", lwd = 1.2)
),
xlim = c(ext(munic)[1], ext(munic)[2]),
ylim = c(ext(munic)[3], ext(munic)[4])
)
13. Bibliografía
Cuaderno original:
Lizarazo, I., 2025. Geomorphometric terrain attributes in R. Disponible en: https://rpubs.com/ials2un/geomorphometric
Paquetes R utilizados:
- Baston, D. (2023). exactextractr: Fast extraction from raster datasets using polygons. R package version 0.10.1.
- Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied spatial data analysis with R (2nd ed.). Springer.
- Cheng, J., Karambelkar, B., & Xie, Y. (2023). leaflet: Create interactive web maps with the JavaScript Leaflet library.
- Hijmans, R. J. (2024). terra: Spatial data analysis. R package version 1.7-78.
- Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446.
- Waga, K., Drăguţ, L., Eisank, C., & Blaschke, T. (2022). MultiscaleDTM: Multiscale analysis of digital terrain models.
Información adicional:
Zia, M. (2023). Digital Terrain Analysis. En Digital Soil Mapping with R. Recuperado de: https://zia207.quarto.pub/digital-terrain-analysis.html
sessionInfo()## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_United States.utf8
## [2] LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] exactextractr_0.10.0 MultiscaleDTM_1.0 terra_1.8-70
## [4] leaflet_2.2.3 sf_1.0-21 sp_2.2-0
##
## loaded via a namespace (and not attached):
## [1] sass_0.4.10 generics_0.1.4 class_7.3-23 KernSmooth_2.23-26
## [5] lattice_0.22-7 digest_0.6.37 magrittr_2.0.4 rgl_1.3.24
## [9] RColorBrewer_1.1-3 evaluate_1.0.5 grid_4.5.1 fastmap_1.2.0
## [13] jsonlite_2.0.0 e1071_1.7-16 DBI_1.2.3 promises_1.5.0
## [17] scales_1.4.0 crosstalk_1.2.2 codetools_0.2-20 jquerylib_0.1.4
## [21] cli_3.6.5 shiny_1.11.1 rlang_1.1.6 units_1.0-0
## [25] base64enc_0.1-3 cachem_1.1.0 yaml_2.3.10 otel_0.2.0
## [29] raster_3.6-32 tools_4.5.1 dplyr_1.1.4 httpuv_1.6.16
## [33] png_0.1-8 vctrs_0.6.5 R6_2.6.1 mime_0.13
## [37] proxy_0.4-27 lifecycle_1.0.4 classInt_0.4-11 htmlwidgets_1.6.4
## [41] pkgconfig_2.0.3 bslib_0.9.0 pillar_1.11.1 later_1.4.4
## [45] glue_1.8.0 Rcpp_1.1.0 xfun_0.53 tibble_3.3.0
## [49] tidyselect_1.2.1 rstudioapi_0.17.1 knitr_1.50 farver_2.1.2
## [53] xtable_1.8-4 htmltools_0.5.8.1 rmarkdown_2.30 compiler_4.5.1Adaptado para La Guajira por: [Julian Santiago
Gonzalez Malpica]
Fecha: 2025-11-04