La geomorfometría, según Olaya, V. (2025) en su libro Sistemas de Información Geográfica, se define como “la ciencia del análisis cuantitativo del relieve”. Es la rama de la geomorfología que se encarga del análisis cuantitativo del relieve terrestre, utilizando principalmente modelos digitales de elevación (DEM). Permite extraer variables como la pendiente, orientación, curvatura, entre otras, que describen la forma y estructura del terreno. Estas variables se emplean en múltiples campos como la hidrología, ecología, geología, planificación territorial y análisis de riesgos.
Este cuaderno busca familiarizar al lector en el análisis geomorfométrico mediante funciones del lenguaje de programación R, destacando la facilidad que ofrecen sus herramientas para extraer, organizar, calcular y visualizar datos. A partir de la replicación y adaptación del cuaderno Geomorphometric Terrain Attributes in R de Lizarazo, Iván (2025), la metodología se aplicó al departamento de Casanare, Colombia, permitiendo así observar características representadas de su morfología y topografía.
. Con este primer código se limpia el entorno de trabajo, asegurando así que el procesamiento de datos no tenga interferencia de algunos trabajos anteriores.
rm(list=ls())
.Para el desarrollo del análisis se utilizó diferentes librerías del lenguaje R, cada una cumple una función específica dentro del entorno de trabajo espacial. El paquete elevatr permite la descarga directa del modelo digital de elevación (DEM) del área de estudio. sf gestiona datos espaciales vectoriales, como los límites administrativos del departamento. terra cimpler como el núcleo del procesamiento raster, el cúal calcula variables geomorfométricas como pendiente, orientación, curvatura y rugosidad. MultiscaleDTM funciona para análisis con funciones específicas para derivar atributos del terreno a distintas escalas. Para la visualización interactiva de los resultados, se utila leaflet, mientras que exactextractr realiza extracciones precisas de datos raster en función de polígonos, mejorando la precisión espacial en el análisis por zonas.
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)
library(dplyr)
list.files("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE")
## [1] "CASANARE.gpkg" "elev_CASANARE_z10.tif" "soc_csnr.gpkg"
Aquí se carga el modelo digital de elevación en formato raster desde un archivo tif. Donde están los dato altitud de la zona de estudio.
(dem = terra::rast("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/elev_CASANARE_z10.tif"))
## class : SpatRaster
## size : 3573, 5128, 1 (nrow, ncol, nlyr)
## resolution : 0.0006856127, 0.0006856127 (x, y)
## extent : -73.125, -69.60918, 4.214913, 6.664608 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_CASANARE_z10.tif
## name : elev_CASANARE_z10
## min value : -789
## max value : 5332
Aquí se reduce la resolución espacial del DEM y se agrupa en bloques de 2*2 pixeles y calcula su valor promedio
dem2 = terra::aggregate(dem,2,"mean")
## |---------|---------|---------|---------|=========================================
|———|———|———|———|=========================================
(munic <- sf::st_read("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/CASANARE.gpkg"))
## Reading layer `municipios' from data source
## `C:\Users\Janus\Desktop\GB2\proyecto3\datos\departamentos\CASANARE\CASANARE.gpkg'
## using driver `GPKG'
## Simple feature collection with 19 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
## Geodetic CRS: MAGNA-SIRGAS
En este bloque se recortan y transforman los datos para ajustarlos al área de estudio (Casanare). Primero, se utiliza la función crop() de terra que recorta el modelo digital de elevación (dem2) al contorno del municipio, eliminando así el área sobrante. Luego, se reproyecta el ráster a un sistema de coordenadas proyectadas más adecuado (EPSG:9377) mediante project(). De forma complementaria, también se transforma el objeto vectorial munic a la misma proyección con st_transform() para asegurar consistencia espacial entre capas ráster y vectoriales. (Ciencia de datos espaciales. 2025).
(dem3 = terra::crop(dem2,munic, mask=TRUE))
## Warning: [crop] CRS do not match
## class : SpatRaster
## size : 1502, 2365, 1 (nrow, ncol, nlyr)
## resolution : 0.001371225, 0.001371225 (x, y)
## extent : -73.07838, -69.83543, 4.286903, 6.346483 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source(s) : memory
## name : elev_CASANARE_z10
## min value : -89.5
## max value : 4423.0
(dem_plane = project(dem3, "EPSG:9377"))
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_CASANARE_z10
## min value : -8.611589
## max value : 4416.087891
(munic_plane = sf::st_transform(munic, "EPSG:9377"))
Aqui se va a calcular la pendiente del terreno (slope). Con multiscale se extraen dos medidas a través del modelo digital de elevación (dem_plane). -1 Slope que como ya se dijo extrae la pendiente en grados, y -2 aspect, que nos muestra la orientación cardinal 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 : 1506, 2370, 2 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 2.159213e-05, 4.141929e-04
## max values : 5.750944e+01, 3.599996e+02
Subset selecciona únicamente la capa del objeto generado que corresponde -1 a la pendiente y -2 a la orientación. mientras que terra::hist genera un histograma donde se ve la distribución de los valores de pendiente en el departamento de Casanare
(slope = subset(slp_asp, 1))
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 2.159213e-05
## max value : 5.750944e+01
terra::hist(slope,
main = "Casanare's slope",
xlab = "Slope (in degrees")
## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
(aspect =subset(slp_asp, 2))
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 4.141929e-04
## max value : 3.599996e+02
terra::hist(aspect,
main = "Casanare's aspect",
xlab = " Aspect (in degrees)")
## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
A partir del raster (slope), se calcula la pendiente expresada como porcentaje mediante la función trigonométrica tan(θ) * 100, donde θ se convierte a radianes. Esta nueva capa (slope_perc), permite una interpretación más fácil del relieve. Luego, se genera un histograma con terra::hist().
(slope_perc = tan(slope*(pi/180))*100)
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 3.768537e-05
## max value : 1.570256e+02
terra::hist(slope_perc,
main = "Casanare's slope",
xlab = "Slope (in percentage)")
## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
En esta estructura se clasifican estas pendientes en siete clases mediante la función classify(), lo que permite categorizar el relieve de acuerdo a intervalos establecidos. Con exact_extract(), se calcula la pendiente media dentro de cada municipio.
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)
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%
## [1] 5.9245358 9.2382898 30.9212017 1.4186283 40.3001289 0.5959594
## [7] 9.1764860 4.7089849 0.5160555 0.7966918 2.8408487 28.1054153
## [13] 15.5886774 28.3318501 0.5730320 25.7137718 4.7331200 0.6306346
## [19] 1.4323320
Warning in .local(x, y, …): Polygons transformed to raster CRS (EPSG:9377)
| | | 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] 2.5160389 1.2292731 1.4036171 5.2911134 0.9474660 0.5868814 [7] 1.5853550 3.0762432 0.3953781 1.7273880 2.0587800 4.6855459 [13] 1.4440339 10.1471510 3.8976312 3.0785139 2.0715215 3.6382058 [19] 16.4352589 2.8243833 2.1665022 2.3823540 3.1995339 1.3972237 [25] 3.4256206 11.3830862 1.5984520 14.8620186 2.7408032 5.6827564
hist(munic$mean_slope,
main = "Casanare´s slope",
xlab = "Slope (in percentage")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%
## [1] 1 1 5 1 5 1 1 1 1 1 1 5 1 5 1 5 1 1 1
Warning in .local(x, y, …): Polygons transformed to raster CRS (EPSG:9377) | | | 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] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1
hist(munic$class,
main = "Casanare´s reclassified slope",
xlab = "Slope (as a category)")
Este codigo proyecta el raster clasificado al sistema de coordenadas geográficas.Esta transoformación es necesaria para visualizar mapas interactivos (leaflet).
(rc.geo = project(rc, "EPSG:4326"))
## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 7
Para finalizar se proyecta el raster de pendientes al sistema de coordenadas EPSG:4326 y se genera una paleta de colores personalizada (Colornumeric) que facilite la visualización de la pendiente.
(slope.geo = project(slope_perc, "EPSG:4326"))
## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 2.620884e-03
## max value : 1.570256e+02
palredgreen <- colorNumeric(c("darkolivegreen3", "goldenrod1", "darkgoldenrod", "peru", "firebrick4"), values(slope.geo),
na.color = "transparent")
leaflet(munic) %>% addTiles() %>% setView(-71.8, 5.2, 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 departamento de Casanar (%)")
(aspect =subset(slp_asp, 2))
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 4.141929e-04
## max value : 3.599996e+02
terra::hist(aspect,
main = "Casanare´s slope",
xlab = "slope (in degrees)")
## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
(slope_perc = tan(slope*(pi/180))*100)
## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 3.768537e-05
## max value : 1.570256e+02
(munic$mean_slope_deg <- exactextractr::exact_extract(slope, munic, 'mean'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%
## [1] 3.2860961 5.1417561 16.8985863 0.8072063 21.3573971 0.3413945
## [7] 5.1038074 2.6620653 0.2956532 0.4544332 1.6183466 15.4698286
## [13] 8.5808420 15.3351622 0.3282120 14.1004839 2.6153166 0.3613016
## [19] 0.8187429
hist(munic$mean_slope_deg,
main = "Pendiente media del Casanare",
xlab = "Pendiente (en grados)")
(munic$class_deg <- exactextractr::exact_extract(rc, munic, 'mode'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%
## [1] 1 1 5 1 5 1 1 1 1 1 1 5 1 5 1 5 1 1 1
hist(munic$class_deg,
main = "Pendiente reclasificada del Casanare",
xlab = "Pendiente (as a category)")
(rc.geo_deg = project(rc, "EPSG:4326"))
## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 7
(slope.geo_dg = project(slope, "EPSG:4326"))
## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.001501656
## max value : 57.509437561
palredgreen_dg <- colorNumeric(c("darkolivegreen3", "goldenrod1", "darkgoldenrod", "peru", "firebrick4"), values(slope.geo),
na.color = "transparent")
leaflet(munic) %>% addTiles() %>% setView(-71.8, 5.2, 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_dg, "<br>")) %>%
addRasterImage(slope.geo_dg, colors = palredgreen_dg, opacity = 0.8) %>%
addLegend(pal = palredgreen_dg, values = values(slope.geo_dg),
title = "Pendiente del departamento de Casanare (in degrees)")
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## 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'
(munic$mean_aspect_dg <- exactextractr::exact_extract(aspect, munic, 'mean'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%
## [1] 154.8160 156.3739 170.9795 160.0436 140.0225 167.1255 165.6081 157.5163
## [9] 167.3810 167.9451 154.0796 158.3428 209.9635 153.1616 165.0070 167.1664
## [17] 157.4100 168.6098 164.1301
(slope.geo_aspdg = project(aspect, "EPSG:4326"))
## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : aspect
## min value : 0.1603621
## max value : 359.9807739
anotherpal <- colorNumeric(c("darkblue", "blue", "cyan", "green", "yellow", "orange", "red"), values(slope.geo),
na.color = "transparent")
leaflet(munic) %>% addTiles() %>% setView(-71.8, 5.2, 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_dg, "<br>")) %>%
addRasterImage(slope.geo_aspdg, colors = anotherpal, opacity = 0.8) %>%
addLegend(pal = anotherpal, values = values(slope.geo_dg),
title = "Aspecto del Terreno para el Departamento de Casanare")
## Warning in pal(c(r[1], cuts, r[2])): Some values were outside the color scale
## and will be treated as NA
## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA
## 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'
```
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.