Atributos Geomorfológicos del Chocó
Laura Valentina Torres Hernández
23.Julio.2025
1. Introducción
Este notebook busca ilustrar los atributos geomorfológicos del departamento del Chocó, Colombia, usando modelos de elevación digital DEM. Se utiliza el paquete Multiscale DTM de R para el procesamiento y análisis de los datos a fin de obtener variables que describan el terreno para comprender de mejor manera la variabilidad topográfica del departamento.
Modelos de Elevación Digital DEMs
Puede definirse como la representación digital del relieve de la tierra, donde a cada punto (Vectorial) o cuadrícula (Raster) le corresponde un valor de elevación. Estos modelos se utilizan para visualizar y analizar la topografía de un terreno, puede usarse en distintas disciplinas, entre ellas la agronomía, ya que permite tomar decisiones informadas según las necesidades del usuario.
2. Configuración del cuaderno
2.1 Limpieza de memoria
Se debe limpiar la memoria para mejorar el rendimiento del código, evitando la confusión de variables o archivos dentro del notebook
rm(list=ls())2.2 Instalación y configuración de paquetes
Inicialmente, se deben instalar los siguientes paquetes, que serán usados a lo largo del código para el análisis de datos geomorfológicos
#install.packages("MultiscaleDTM")
#install.packages("exactextractr")Adicionalmente se requiere del uso de las siguientes librerías:
# Librería que permite el procesamiento de datos de elevación
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.# Manejo de datos espaciales
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUE# Visualización interactiva de los datos
library(leaflet)
# Librería para la manipulación de datos ráster
library(terra)## terra 1.8.54# Cálculo de atributos geomorfométricos del terreno
library(MultiscaleDTM)## Warning: package 'MultiscaleDTM' was built under R version 4.5.1# Extracción de valores de ráster dentro de áreas específicas
library(exactextractr)3. Datos de elevación del Chocó
A continuación, se usará la base de datos trabajada a lo largo del curso de Geomática Básica.
Inicialmente, se carga el DEM raster del Chocó obtenido de un documento anteriormente trabajado. Se lee DEM usando la librería Terra, este archivo se llamará “dem”
(dem = terra::rast("~/UNAL Valentina/GB2/p1/datos/elev_choco_z10.tif"))## class : SpatRaster
## size : 7161, 3089, 1 (nrow, ncol, nlyr)
## resolution : 0.0006829542, 0.0006829542 (x, y)
## extent : -78.04687, -75.93723, 3.86416, 8.754795 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_choco_z10.tif
## name : elev_choco_z10
## min value : -4485
## max value : 4684Propiedades del documento “dem”: Tiene clase SpatRaster, ideal para trabajar con datos de tipo espacial raster. Posee una tabla con un total de 7161 filas, 3089 columnas y una sola capa. Cada celda tiene una resolución de aproximadamente 0,00068 grados en (x,y), es decir, se tiene un nivel de detalle bastante alto.
El sistema de referencia espacial es geográfico, con un elipsoide GRS80.
Reducir la resolución del DEM, para evitar problemáticas en la memoria, ya que se agrupan las celdas del archivo DEM original en bloques de 2*2, promediando valores, así, se reduce la resolución y por lo tanto, la carga de la memoria del sistema:
dem2= terra::aggregate(dem,2,"mean")## |---------|---------|---------|---------|========================================= Se incluye el geopackage del departamento con las limitaciones del departamento. Se usa la librería sf para esta acción
(munic <- sf::st_read("~/UNAL Valentina/GB2/p1/datos/MGN2023_MPIO_POLITICO/choco.gpkg"))## Reading layer `choco' from data source
## `C:\Users\RAUL TORRES MOYA\Documents\UNAL Valentina\GB2\p1\datos\MGN2023_MPIO_POLITICO\choco.gpkg'
## using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.88378 ymin: 3.964883 xmax: -76.00185 ymax: 8.67773
## Geodetic CRS: MAGNA-SIRGASPreguardado de los archivos en un nuevo atributo
dem3 = terra::crop(dem2, munic, mask=TRUE)## Warning: [crop] CRS do not match#3.1 Transformación de coordenadas para calcular correctamente la geomorgfología del departamento el DEM debe tener coordenadas planas. Para el año 2020, el Instituto Geográfico Agustin Codazzi (IGAC) adopta el sistema de coordenadas plana con origen único nacional en Colombia, con código EPSG 9377. Se usa el siguiente código para la transformación de coordenadas:
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 3462, 1404, 1 (nrow, ncol, nlyr)
## resolution : 151.2992, 151.2992 (x, y)
## extent : 4457349, 4669773, 1996619, 2520417 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_choco_z10
## min value : -1084.250
## max value : 4306.362Como se observa anteriormente, la transformación de coordenadas fue exitosa.
Se transforma el objeto a un vectorial, con el siguiente código:
(munic_plane= sf::st_transform(munic, "EPSG:9377"))4.Cálculo de los atributos geomorfométricos.
Para calcular el aspecto y la pendiente, se debe usar la función SlpAsp de la librería MultiscaleDTM. Inicialmente, se define el parámetro w como el tamaño de la ventana usada para calcular los atributos. Para este caso, w=c(3,3), es decir, una ventana con celdas de 3*3. Continuando con e código, el parámetro method define eñ tipo de vecindad utilizado para los cálculos, donde “queen” es un tipo de vecindad que incluye aquellos datos que pueden llamarse vecinos directos o intermedios.
(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 : 3462, 1404, 2 (nrow, ncol, nlyr)
## resolution : 151.2992, 151.2992 (x, y)
## extent : 4457349, 4669773, 1996619, 2520417 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 0
## max values : 84.12251, 360Se deben separar los resultados del código anterior en dos capas: Capa 1: Capa de la pendiente.
(slope = subset(slp_asp,1))## class : SpatRaster
## size : 3462, 1404, 1 (nrow, ncol, nlyr)
## resolution : 151.2992, 151.2992 (x, y)
## extent : 4457349, 4669773, 1996619, 2520417 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 84.12251La capa de pendiente, puede visualizarse a continuación:
terra::hist(slope, main= "Pendiente del Chocó", xlab = "Pendiente (en grados)")## Warning: [hist] a sample of 21% of the cells was used (of which 56% was NA)
Capa 2: Orientación
(aspect = subset(slp_asp,2))## class : SpatRaster
## size : 3462, 1404, 1 (nrow, ncol, nlyr)
## resolution : 151.2992, 151.2992 (x, y)
## extent : 4457349, 4669773, 1996619, 2520417 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 360La capa orientación, puede visualizarse a continuación:
terra::hist(aspect, main= "Aspecto del Chocó", xlab= "Aspecto(en grados)")## Warning: [hist] a sample of 21% of the cells was used (of which 56% was NA)
Se debe convertir la pendiente de grados a porcentaje:
(slope_perc = tan(slope* (pi/180))*100)## class : SpatRaster
## size : 3462, 1404, 1 (nrow, ncol, nlyr)
## resolution : 151.2992, 151.2992 (x, y)
## extent : 4457349, 4669773, 1996619, 2520417 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 971.4124#gráfica con la nueva conversión
terra:: hist(slope_perc, main= "Pendiente del Chocó", xlab= "% Pendiente")## Warning: [hist] a sample of 21% of the cells was used (of which 56% was NA)
5.Cálculo de estadísticas zonales en el Chocó
Inicialmente, para analizar las pendientes en el departamento, se utiliza la clasificación diseñada por el IGAC con fines agrológicos, que puede observarse a continuación:
Esta clasificación divide los valores de pendiente en diferentes rangos,
representando la dificultad del manejo del terreno Para realizar la
computación de zonas estáticas parala transformación de los porcentajes
continuos de la pendiente a intervalos categorizados se utiliza la
siguiente matriz, a fin de obtener un nuevo ráster reclasificado que
permite identificar las áreas con pendientes similares y categorizarlas
según su inclinación.
#Reclasificación de la pendiente raster
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)Cálculo de las estadísticas zonales usando exactextractr
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## 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] 15.020624 12.928195 2.764322 30.254103 11.365789 33.327873 15.111613
## [8] 8.296690 14.929053 5.924431 15.434294 10.987267 19.116861 8.273375
## [15] 4.182716 6.839398 4.598071 2.769994 22.241301 5.588526 13.625473
## [22] 6.729395 11.251538 4.060682 3.106334 16.478111 13.420687 12.449834
## [29] 38.972328 25.426815Visualización gráfica mediante histograma de la pendiente media en el Chocó
hist(munic$mean_slope, main= "Pendiente media del Chocó", xlab="% Pendiente")
Pendiente predominante de los municipios en el Chocó
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## 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] 4 4 1 5 4 5 4 1 4 1 4 1 4 1 1 2 2 1 4 1 4 1 1 1 1 4 4 1 5 5hist(munic$class, main = "Pendiente reclasificada del Chocó", xlab="Pendiente (Y su categoría)")
Transformación de la pendiente a coordenadas geográficas
(rc.geo= project(rc,"EPSG:4326"))## class : SpatRaster
## size : 3472, 1431, 1 (nrow, ncol, nlyr)
## resolution : 0.001365887, 0.001365887 (x, y)
## extent : -77.92888, -75.9743, 3.95545, 8.697809 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 918.6257Luego de obtener las coordenadas geográficas, estas se convierten a porcentaje de pendiente
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 3472, 1431, 1 (nrow, ncol, nlyr)
## resolution : 0.001365887, 0.001365887 (x, y)
## extent : -77.92888, -75.9743, 3.95545, 8.697809 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 918.62576.Representación gráfica de la pendiente del terreno
Al querer graficar los procedimientos hechos anteriormente, se genera un “mapa” interactivo, con la intención de visualizar la pendiente de los municipios en el Chocó. Es bastante útil ya que permite comprender las características topográficas de mejor manera y generar un material esencial para tomar decisiones relacionadas con el uso de los suelos y la planificación de actividades en el territorio.
Para la configuración del mapa interactivo, inicialmente se definen los colores para los diferentes niveles de la pendiente, usando el más bajo como aguamarina, y el más alto como rojo oscuro. Se identifican los límites municipales con un borde de coor gris, y la opacidad del relleno es de un 40%, además, al interactuar el mapa haciendo clic en algún municipio, se desplegará el cuadro emergente con el nombre del municipio y la clase de pendiente más frecuente.
palredgreen <- colorNumeric(c("aquamarine","chartreuse","darkseagreen3","yellow2", "orange", "brown2","chocolate", "darkred"), values(slope.geo),
na.color = "transparent")
#Mapa interactivo del Chocó
leaflet(munic) %>% addTiles() %>% setView( -75.9743,8.697809 , 9) %>%
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 el Chocó (%)")## 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'