Atributos geomorfométricos en R
23.07.2026
1. Introducción
En este cuaderno veremos cómo calcular atributos geomorfológicos de un tereno a partir de un DEM (Modelo de elevación digital)
2. Configuración
Primero, hay que limpiar la memoria del servidor.
rm(list = ls())Luego, cargamos las librerias necesarias para este cuaderno, es importante descargarlas desde la consola.
library(elevatr)## 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)## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.9.27library(MultiscaleDTM)## Warning: package 'MultiscaleDTM' was built under R version 4.6.1library(exactextractr)## Warning: package 'exactextractr' was built under R version 4.6.13. Lectura de los datos de entrada
Necesitamos un documento DEM de nuestro departamento (Santander), en este caso, nuestro DEM fue obtenido de la ejecución del cuaderno “Elevación y análisis de datos altimétricos en R”
(dem= terra::rast("./elev_stder_z10.tif"))## class : SpatRaster
## size : 4090, 3607, 1 (nrow, ncol, nlyr)
## resolution : 75.52857, 75.52857 (x, y)
## extent : 4791535, 5063967, 2178575, 2487487 (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_stder_z10.tif
## name : elev_stder_z10
## min value : -1458
## max value : 4511Reducimos la resolución del DEM para evitar problemas en la memoria
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Leemos nuestros municipios usando la libreria sf (simple features)
(munic <- sf::st_read("./Municipios__9377.gpkg"))## Reading layer `muncipios9377' from data source
## `C:\Users\isabe\OneDrive\Documentos\GB2\proyecto 2\Datos\Municipios__9377.gpkg'
## using driver `GPKG'
## Simple feature collection with 87 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4828392 ymin: 2188513 xmax: 5053451 ymax: 2455407
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-NacionalCon la ayuda de dem2 establecemos los límites de nuestro departamento
dem3 = terra::crop(dem2,munic, mask=TRUE)## Warning: [crop] CRS do not match##4. Transformación de coordenadas
Para este cuaderno es necesario tener coordenadas planas, por eso lo convertimos al EPSG:9377 MAGNA SIRGAS, sistema de referecnia de Colombia
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1767, 1490, 1 (nrow, ncol, nlyr)
## resolution : 151.0571, 151.0571 (x, y)
## extent : 4828393, 5053469, 2188545, 2455463 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_stder_z10
## min value : -138.25
## max value : 4413.25(munic_plane = sf::st_transform(munic, "EPSG:9377"))##5. Cálculo de atributos del terreno Los atributos a calcular son pendiente y orientación, definiendo 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 : 1767, 1490, 2 (nrow, ncol, nlyr)
## resolution : 151.0571, 151.0571 (x, y)
## extent : 4828393, 5053469, 2188545, 2455463 (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 : 68.999358, 359.986137Vamos a dividir el slp_asp en dos partes: Capa 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1767, 1490, 1 (nrow, ncol, nlyr)
## resolution : 151.0571, 151.0571 (x, y)
## extent : 4828393, 5053469, 2188545, 2455463 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 68.999358terra::hist(slope,
main = "Santander's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 38% of the cells was used (of which 50% was NA)
Capa 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1767, 1490, 1 (nrow, ncol, nlyr)
## resolution : 151.0571, 151.0571 (x, y)
## extent : 4828393, 5053469, 2188545, 2455463 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 359.986137terra::hist(aspect,
main = "Santander's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 38% of the cells was used (of which 50% was NA)
Además vamos a convertir la pendiente de grados a porcentaje:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1767, 1490, 1 (nrow, ncol, nlyr)
## resolution : 151.0571, 151.0571 (x, y)
## extent : 4828393, 5053469, 2188545, 2455463 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 260.500187terra::hist(slope_perc,
main = "Santander's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 38% of the cells was used (of which 50% was NA)
6. Cálculo de estadísticas computacionales
Vamos a clasificar los datos de pendiente de nuestro departamento conociendo la clasificación de pendiente hecha por el IGAC
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)Vamos a calcular las estadísticas zonales con exactextractr
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## | | | 0% | |= | 1% | |== | 2% | |== | 3% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |=================== | 26% | |=================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |=================================== | 49% | |=================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 31.321646 21.353863 30.627333 13.014132 19.820490 2.225787 29.110426
## [8] 14.965528 21.825090 20.713629 38.334518 32.808846 34.670158 44.055855
## [15] 28.349388 22.879078 37.937775 36.332558 20.610285 4.358469 28.797840
## [22] 14.830569 35.288429 37.843994 19.448462 19.399921 27.172073 28.501127
## [29] 31.664316 36.835224 31.420605 33.241192 21.725061 21.003988 15.499665
## [36] 40.921448 19.696968 37.789570 16.215914 16.230015 15.926795 40.528732
## [43] 18.842348 30.603006 26.678057 21.462376 24.676245 16.312748 28.232683
## [50] 30.218517 37.722324 36.439442 26.936869 43.165722 22.694925 14.593230
## [57] 39.865593 17.375193 24.758087 15.237010 29.483728 20.544241 15.354102
## [64] 3.138125 1.254872 13.708508 2.226363 41.695950 21.377132 20.223143
## [71] 36.614407 45.185490 35.478828 19.305529 46.783909 23.287184 14.771005
## [78] 16.493189 17.321098 23.217216 41.619083 31.110353 21.058552 18.440647
## [85] 40.554928 14.688760 32.418762hist(munic$mean_slope,
main = "Santander's mean slope",
xlab = "Slope (in percentage)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |= | 1% | |== | 2% | |== | 3% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |=================== | 26% | |=================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |======================= | 32% | |======================= | 33% | |======================== | 34% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |=================================== | 49% | |=================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |=================================================== | 72% | |=================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 5 4 5 4 4 1 5 1 4 4 5 5 5 5 5 4 5 5 4 1 5 4 5 5 4 4 4 5 5 5 5 5 5 4 4 5 4 5
## [39] 4 4 4 5 4 5 5 4 5 4 5 5 5 5 5 5 4 4 5 4 4 4 5 4 4 1 1 1 1 5 4 4 5 5 5 4 5 4
## [77] 1 4 4 4 5 5 4 4 5 4 5hist(munic$class,
main= "Santander's reclassified slope",
xlab= "Slope (as a category)")
Ahora, transformamos las coordenadas geográficas
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1768, 1495, 1 (nrow, ncol, nlyr)
## resolution : 0.001367307, 0.001367307 (x, y)
## extent : -74.55813, -72.51401, 5.704487, 8.121886 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 244.552322(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1768, 1495, 1 (nrow, ncol, nlyr)
## resolution : 0.001367307, 0.001367307 (x, y)
## extent : -74.55813, -72.51401, 5.704487, 8.121886 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 244.5523227. Plotteando la pendiente del terreno
Primero debemos definir nuestra paleta de colores para la pendiente
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ejecutamos nuestro plotteo
leaflet(munic) %>% addTiles() %>% setView(-73, 7, 9) %>% addPolygons(color="gray", weight = 1.0, smoothFactor = 0.5,
opacity=0.4, fillOpacity = 0.10,
popup = paste("Municipio:", munic$NAME_2, "<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 Santander(%)")## Warning: sf layer is not long-lat data## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'##8. Visualización alternativa
Con este código tenemos una alternativa para observar la pendiente de nuestro territorio
slope.geo <- subset(slp_asp, "slope")
leaflet(munic) %>%
addTiles() %>%
setView(-73, 7, 9) %>%
addPolygons(
color = "gray",
weight = 1,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste(
"Municipio:", munic$NAME_2
)
) %>%
addRasterImage(
slope.geo,
colors = palredgreen,
opacity = 0.8
) %>%
addLegend(
pal = palredgreen,
values = terra::values(slope.geo),
title = "Terrain slope (degrees)"
)## Warning: sf layer is not long-lat data## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'Visualización alternativa de la orientación
Con este código creamos una visualización alternativa de la orientación.
Creamos una paleta
aspect.geo <- subset(slp_asp, "aspect")
pal_aspect <- colorNumeric(
palette = hcl.colors(12, "Spectral"),
domain = terra::values(aspect.geo),
na.color = "transparent"
)Y ejecutamos
leaflet(munic) %>%
addTiles() %>%
setView(-73, 7, 9) %>%
addPolygons(
color = "gray",
weight = 1,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste(
"Municipio:", munic$NAME_2
)
) %>%
addRasterImage(
aspect.geo,
colors = pal_aspect,
opacity = 0.8
) %>%
addLegend(
pal = pal_aspect,
values = terra::values(aspect.geo),
title = "Terrain aspect (degrees)")## Warning: sf layer is not long-lat data## Warning: sf layer has inconsistent datum (+proj=tmerc +lat_0=4 +lon_0=-73 +k=0.9992 +x_0=5000000 +y_0=2000000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs).
## Need '+proj=longlat +datum=WGS84'