Atributos geomorfométricos del terreno en el Putumayo

1 1 Introducción

Este cuaderno tiene como objetivo calcular atributos geomorfométricos del terreno a partir de modelos digitales de elevación (MDE) para el departamento del Putumayo, Colombia. Estos atributos permiten describir las características físicas del relieve, como la pendiente, el aspecto y otras formas del terreno.

El análisis se realiza utilizando herramientas del lenguaje R, especialmente el paquete MultiscaleDTM, que permite derivar atributos a diferentes escalas espaciales. Este trabajo está orientado a fortalecer el aprendizaje práctico de conceptos en geoinformática, como parte de la formación en ciencias agronómicas de la Universidad Nacional de Colombia.

2 2 Configuración del entorno de trabajo

En esta sección se prepara el entorno de trabajo necesario para calcular atributos geomorfométricos a partir de modelos digitales de elevación (DEM). El proceso inicia con la limpieza de la memoria y continúa con la carga de las bibliotecas que permiten la manipulación, análisis y visualización de datos espaciales.

2.0.0.1 2.1 Limpieza de memoria

Esta línea elimina todos los objetos almacenados en la memoria de trabajo de R. Esto es útil para comenzar con un entorno limpio.

rm(list = ls())

2.0.0.2 2.2 Carga de bibliotecas necesarias

A continuación, se cargan los paquetes necesarios para ejecutar los procedimientos del cuaderno. Cada uno de estos paquetes cumple una función específica dentro del análisis: elevatr: permite descargar modelos digitales de elevación (DEM) directamente desde fuentes en línea como AWS o USGS, facilitando el acceso a datos topográficos.

¨sf: proporciona herramientas para manipular objetos espaciales vectoriales con geometría simple (shapefiles, puntos, líneas, polígonos), utilizando una sintaxis moderna y eficiente.

leaflet: se utiliza para la creación de mapas interactivos directamente en R, integrando fácilmente capas ráster y vectoriales.

terra: sirve para la lectura, manipulación y análisis de datos ráster, siendo esencial para trabajar con modelos digitales de elevación y derivar atributos del terreno como pendiente y aspecto.

MultiscaleDTM: es el paquete clave en este cuaderno, ya que permite calcular atributos geomorfométricos del terreno a múltiples escalas, como rugosidad, pendiente o curvatura.

exactextractr: permite extraer estadísticas de valores ráster sobre polígonos vectoriales con alta precisión, algo fundamental para el análisis zonal de atributos del terreno.

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.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE

library(leaflet)
library(terra)

## terra 1.8.29

library(MultiscaleDTM)
library(exactextractr)

3 3 Introduction to MultiscaleDTM

El paquete MultiscaleDTM en R permite calcular atributos geomorfométricos del terreno (como pendiente, orientación y curvatura)parámetros importantes para describir la forma de una superficie, a múltiples escalas, a partir de Modelos Digitales de Terreno DTM . Estos atributos son calculados utilizando ventanas móviles centradas en cada celda del raster de elevación.

El paquete MultiscaleDTM permite aplicar algoritmos sobre una celda focal (Z₀₀) y un conjunto de celdas vecinas, cuya selección depende del método usado:

Método rook: utiliza únicamente las celdas cardinales.

Método queen: incluye además las celdas diagonales, es decir, todas las 8 celdas vecinas.

Método boundary: considera todas las celdas ubicadas en el borde externo de la ventana.

El cálculo multiescala implica repetir este análisis usando diferentes tamaños de ventana (por ejemplo, 3x3, 5x5, 7x7, etc.), lo que permite capturar patrones topográficos a distintas resoluciones. Este enfoque es útil en aplicaciones como clasificación de uso del suelo, delimitación de formas del terreno, modelado hidrológico y análisis de conectividad ecológica.

En este cuaderno, usaremos MultiscaleDTM para derivar atributos como la pendiente, el aspecto (dirección de exposición) y la curvatura, observando cómo cambian a diferentes escalas.

knitr::include_graphics("Z(0,0)_.jpg")

4 4. Lectura de los datos de entrada

En esta sección se cargan los datos necesarios para calcular atributos del terreno. Se utiliza un Modelo Digital de Elevación (DEM) previamente obtenido y una capa vectorial de municipios del departamento de Córdoba. A continuación, se explican cada uno de los pasos realizados en el código:

4.0.0.1 4.1 Cargar el DEM usando el paquete terra

La función list.files() permite visualizar la lista de archivos disponibles en una carpeta específica, útil para verificar qué datos están disponibles antes de cargarlos.

list.files("./putumayo")

## [1] "elev_putumayo_z10.tif" "putumayo.gpkg"

Con este comando se carga un archivo raster (.tif) que representa el Modelo Digital de Elevación (DEM) del departamento de Córdoba. El paquete terra permite manipular este tipo de datos espaciales raster. La función rast() carga el archivo y lo almacena como un objeto tipo SpatRaster. para que se carga que funcion cumple este archivo

(dem = terra::rast("./putumayo/elev_putumayo_z10.tif"))

## class       : SpatRaster 
## dimensions  : 3584, 5120, 1  (nrow, ncol, nlyr)
## resolution  : 0.0006866186, 0.0006866186  (x, y)
## extent      : -77.34375, -73.82826, -0.7033042, 1.757537  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs 
## source      : elev_putumayo_z10.tif 
## name        : elev_putumayo_z10 
## min value   :              -540 
## max value   :              4186

4.0.0.2 4.2. Reducir la resolución del DEM

Con este código se reduce la resolución del raster.La función aggregate() toma bloques de celdas (en este caso, bloques de 2x2 píxeles) y calcula el valor promedio (“mean”), generando así un nuevo raster (dem2) con menor resolución espacial.

dem2 = terra::aggregate(dem,2, "mean")

## 
|---------|---------|---------|---------|
=========================================
                                          

4.0.0.3 4.3 Cargar el shapefile de municipios usando sf

Este código carga un archivo vectorial en formato .gpkg (Geopackage) que contiene la información geográfica de los municipios del departamento de Córdoba. Para que se carga este arcgivo?

La función st_read() del paquete sf lee los datos espaciales y los convierte en un objeto sf (Simple Features), permitiendo su manipulación en R.

En la tabla resultante se observan los códigos y nombres de cada municipio, junto con su geometría (polígonos que representan sus límites).

(munic <-  sf::st_read("putumayo/putumayo.gpkg"))

## Reading layer `putumayo' from data source 
##   `D:\GB2\P5_semana10_geomorfómetria\putumayo\putumayo.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 13 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.18681 ymin: -0.5622776 xmax: -73.84132 ymax: 1.467315
## Geodetic CRS:  MAGNA-SIRGAS

## Simple feature collection with 13 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.18681 ymin: -0.5622776 xmax: -73.84132 ymax: 1.467315
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr       mpio_cnmbr
## 1          86        001      86001   PUTUMAYO            MOCOA
## 2          86        219      86219   PUTUMAYO            COLÓN
## 3          86        320      86320   PUTUMAYO            ORITO
## 4          86        568      86568   PUTUMAYO      PUERTO ASÍS
## 5          86        569      86569   PUTUMAYO   PUERTO CAICEDO
## 6          86        571      86571   PUTUMAYO    PUERTO GUZMÁN
## 7          86        573      86573   PUTUMAYO PUERTO LEGUÍZAMO
## 8          86        749      86749   PUTUMAYO         SIBUNDOY
## 9          86        755      86755   PUTUMAYO    SAN FRANCISCO
## 10         86        757      86757   PUTUMAYO       SAN MIGUEL
##                              mpio_crslc mpio_tipo  mpio_narea mpio_nano
## 1                                  1958 MUNICIPIO  1304.63788      2023
## 2   Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO    64.27462      2023
## 3  Decreto 2891 de Diciembre 28 de 1978 MUNICIPIO  1939.35399      2023
## 4    Decreto 1951 de Octubre 24 de 1967 MUNICIPIO  2819.15478      2023
## 5  Ordenanza 12 de Noviembre 24 de 1992 MUNICIPIO   926.50246      2023
## 6  Ordenanza 13 de Noviembre 24 de 1992 MUNICIPIO  4576.59156      2023
## 7   Decreto 698 de Noviembre 13 de 1953 MUNICIPIO 10906.87681      2023
## 8       Decreto 1871 de Julio 1 de 1982 MUNICIPIO    97.73462      2023
## 9   Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO   407.35674      2023
## 10     Ordenanza 45 de Abril 29 de 1994 MUNICIPIO   379.74213      2023
##    shape_Leng  shape_Area                           geom
## 1   2.4752281 0.105792947 MULTIPOLYGON (((-76.6705 1....
## 2   0.4548864 0.005209533 MULTIPOLYGON (((-76.96835 1...
## 3   2.1523001 0.157168461 MULTIPOLYGON (((-76.76607 0...
## 4   3.6444021 0.228679089 MULTIPOLYGON (((-76.2263 0....
## 5   1.7199295 0.075142214 MULTIPOLYGON (((-76.46575 0...
## 6   4.7346323 0.371458749 MULTIPOLYGON (((-75.94095 1...
## 7   7.5377279 0.885757929 MULTIPOLYGON (((-75.20032 0...
## 8   0.5113819 0.007922269 MULTIPOLYGON (((-76.9043 1....
## 9   1.1754950 0.033022563 MULTIPOLYGON (((-76.87345 1...
## 10  1.3275843 0.030777834 MULTIPOLYGON (((-76.99677 0...

4.0.0.4 4.4 Recortar el DEM según los límites del departamento

Este paso permite recortar el raster (dem2) para que solo contenga información dentro de los límites de los municipios de Córdoba.

La función crop() utiliza el shapefile (munic) como máscara.

El argumento mask=TRUE indica que se excluya todo lo que esté por fuera del área delimitada por los municipios.

dem3 = terra::crop(dem2,munic, mask=TRUE)

Tanto el raster (dem2) como el shapefile (munic) se encuentran en coordenadas geográficas, lo que permite su compatibilidad para ser cruzados y utilizados en conjunto.

5 5 Transforming coordinates (10%)

Para continuar con el análisis geomorfométrico de atributos del terreno, es fundamental que el sistema de coordenadas de los datos utilizados esté en coordenadas planas, ya que muchas funciones de cálculo (como la pendiente y la curvatura) requieren trabajar en unidades métricas (metros y no grados).

Desde el año 2020, el IGAC adoptó oficialmente el sistema de coordenadas planas con origen único nacional, que corresponde al EPSG:9377, también conocido como MAGNA-SIRGAS / Origen Nacional. Por lo tanto, es necesario transformar tanto el raster (DEM) como el vector de municipios a este sistema de referencia.

5.0.1 5.1 Proyección del modelo digital de elevación (DEM)

Aquí se utiliza la función project() del paquete terra para transformar el sistema de coordenadas del objeto dem3 (DEM recortado para el departamento) al EPSG:9377. Este paso garantiza que las mediciones y análisis posteriores se realicen en metros esencial para obtener resultados precisos al calcular atributos del terreno. Eue se obtiene de este codigo?, describir resultados.

(dem_plane = project(dem3, "EPSG:9377"))

## class       : SpatRaster 
## dimensions  : 1473, 2441, 1  (nrow, ncol, nlyr)
## resolution  : 152.6319, 152.6319  (x, y)
## extent      : 4533838, 4906412, 1495829, 1720656  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        : elev_putumayo_z10 
## min value   :         -61.38813 
## max value   :        3892.48047

5.0.2 5.2 Transformación del shapefile de municipios

En este segundo paso se transforma el archivo vectorial que contiene los municipios (munic) utilizando la función st_transform() del paquete sf. Esto asegura que tanto los datos raster como vectoriales estén en el mismo sistema de coordenadas, un requisito indispensable para realizar operaciones como recortes y visualizaciones. Describir resultados.

(munic_plane = sf::st_transform(munic, "EPSG:9377"))

## Simple feature collection with 13 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4533910 ymin: 1495907 xmax: 4906419 ymax: 1720503
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr       mpio_cnmbr
## 1          86        001      86001   PUTUMAYO            MOCOA
## 2          86        219      86219   PUTUMAYO            COLÓN
## 3          86        320      86320   PUTUMAYO            ORITO
## 4          86        568      86568   PUTUMAYO      PUERTO ASÍS
## 5          86        569      86569   PUTUMAYO   PUERTO CAICEDO
## 6          86        571      86571   PUTUMAYO    PUERTO GUZMÁN
## 7          86        573      86573   PUTUMAYO PUERTO LEGUÍZAMO
## 8          86        749      86749   PUTUMAYO         SIBUNDOY
## 9          86        755      86755   PUTUMAYO    SAN FRANCISCO
## 10         86        757      86757   PUTUMAYO       SAN MIGUEL
##                              mpio_crslc mpio_tipo  mpio_narea mpio_nano
## 1                                  1958 MUNICIPIO  1304.63788      2023
## 2   Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO    64.27462      2023
## 3  Decreto 2891 de Diciembre 28 de 1978 MUNICIPIO  1939.35399      2023
## 4    Decreto 1951 de Octubre 24 de 1967 MUNICIPIO  2819.15478      2023
## 5  Ordenanza 12 de Noviembre 24 de 1992 MUNICIPIO   926.50246      2023
## 6  Ordenanza 13 de Noviembre 24 de 1992 MUNICIPIO  4576.59156      2023
## 7   Decreto 698 de Noviembre 13 de 1953 MUNICIPIO 10906.87681      2023
## 8       Decreto 1871 de Julio 1 de 1982 MUNICIPIO    97.73462      2023
## 9   Decreto 2830 de Diciembre 2 de 1989 MUNICIPIO   407.35674      2023
## 10     Ordenanza 45 de Abril 29 de 1994 MUNICIPIO   379.74213      2023
##    shape_Leng  shape_Area                           geom
## 1   2.4752281 0.105792947 MULTIPOLYGON (((4591580 172...
## 2   0.4548864 0.005209533 MULTIPOLYGON (((4558354 170...
## 3   2.1523001 0.157168461 MULTIPOLYGON (((4580834 164...
## 4   3.6444021 0.228679089 MULTIPOLYGON (((4640969 162...
## 5   1.7199295 0.075142214 MULTIPOLYGON (((4614310 165...
## 6   4.7346323 0.371458749 MULTIPOLYGON (((4672785 167...
## 7   7.5377279 0.885757929 MULTIPOLYGON (((4755206 161...
## 8   0.5113819 0.007922269 MULTIPOLYGON (((4565496 170...
## 9   1.1754950 0.033022563 MULTIPOLYGON (((4568933 170...
## 10  1.3275843 0.030777834 MULTIPOLYGON (((4555084 159...

6 6. Computing terrain attributes (20%)

En esta sección del cuaderno, se calcula la pendiente (slope) y la exposición (aspect) del terreno a partir del Modelo de Elevación Digital (DEM) previamente reproyectado a coordenadas planas. Estos atributos son fundamentales en análisis geomorfológicos y estudios ambientales, ya que permiten caracterizar la forma del relieve y su influencia sobre procesos naturales y antrópicos (como escorrentía, erosión, orientación solar, entre otros).

6.0.0.1 6.1. Cálculo multiescala de pendiente y exposición

Se utiliza la función SlpAsp() del paquete MultiscaleDTM para calcular ambos atributos de forma multiescala, es decir, considerando diferentes tamaños de ventana en este caso, una matriz 3x3.

El resultado es un objeto tipo SpatRaster con dos capas: una de pendiente y otra de exposición. dDscribir resultados que se generan

(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 
## dimensions  : 1473, 2441, 2  (nrow, ncol, nlyr)
## resolution  : 152.6319, 152.6319  (x, y)
## extent      : 4533838, 4906412, 1495829, 1720656  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## names       :    slope,   aspect 
## min values  :  0.00000,   0.0000 
## max values  : 54.96119, 359.9999

6.0.0.2 6.2. División del objeto en dos capas

A continuación, se separan las capas del objeto `slp_asp en slope y aspect para analizarlas por separado.

Esto permite realizar análisis específicos sobre cada variable y visualizar su distribución.

Capa 1

Este código nos ayuda a aislar la información de la pendiente del objeto multicapas slp_asp, para poder estudiarla y procesarla de forma independiente.La función subset() se utiliza aquí para seleccionar únicamente la primera capa, que en este caso corresponde a la pendiente del terreno en grados, slp_asp es el objeto multicapas generado por la función SlpAsp() en el punto anterior, 1 indica que se selecciona la primera capa del objeto.

(slope = subset(slp_asp, 1))

## class       : SpatRaster 
## dimensions  : 1473, 2441, 1  (nrow, ncol, nlyr)
## resolution  : 152.6319, 152.6319  (x, y)
## extent      : 4533838, 4906412, 1495829, 1720656  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 54.96119

Aquí se genera un histograma con la distribución de los valores de pendiente en grados en el area de estudio.El histograma de la pendiente en Putumayo muestra una distribución altamente sesgada hacia la izquierda, donde la mayoría de los valores se concentran entre 0° y 10°. Esto indica que gran parte del terreno presenta pendientes suaves o casi planas, típicas de zonas de valle, terrazas aluviales o áreas ligeramente onduladas.

A medida que aumenta el ángulo de pendiente, la frecuencia disminuye notablemente, lo cual sugiere que las pendientes pronunciadas (mayores de 20°) son escasas.

terra::hist(slope, 
     main = "Putumayo's slope", 
     xlab = "Slope (in degrees)")

## Warning: [hist] a sample of 28% of the cells was used (of which 69% was NA)

Capa 2:

Este código es similar al de la capa anterior, con la única diferencia de que aquí se selecciona la segunda capa del objeto slp_asp. Este cambio permite extraer ahora el aspecto en lugar de la pendiente.

El aspecto del terreno se refiere a la dirección hacia la cual se orientan las pendientes (en grados), medida desde el norte en sentido horario (0° a 360°). Esta variable es fundamental para estudios ambientales, agrícolas y de erosión, ya que influye en la exposición al sol, temperatura del suelo, evapotranspiración y tipo de vegetación.

(aspect =subset(slp_asp, 2))

## class       : SpatRaster 
## dimensions  : 1473, 2441, 1  (nrow, ncol, nlyr)
## resolution  : 152.6319, 152.6319  (x, y)
## extent      : 4533838, 4906412, 1495829, 1720656  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :   aspect 
## min value   :   0.0000 
## max value   : 359.9999

El histograma generado muestra una distribución relativamente uniforme de los valores de aspecto, con algunas variaciones en ciertos rangos. Esto indica que las pendientes en el territorio de Putumayo están orientadas hacia múltiples direcciones, sin una orientación predominante marcada. Las frecuencias más altas se observan entre 90° y 180°, que corresponden a pendientes orientadas al este y sureste. Sin embargo, también hay una presencia considerable en otras direcciones.

terra::hist(aspect, 
     main = "Putumayo's aspect", 
     xlab = "Aspect (in degrees)")

## Warning: [hist] a sample of 28% of the cells was used (of which 69% was NA)

#### 6.3 Conversión de pendiente de grados a porcentaje} Finalmente, se convierte la pendiente de grados a porcentaje. Una vez aplicada la función tan() al raster de pendiente (slope) y multiplicado por 100, se obtiene un nuevo raster denominado slope_perc, que representa la pendiente expresada en porcentaje. El resultado conserva las características espaciales del raster original, incluyendo dimensiones, resolución, extensión geográfica y sistema de referencia (MAGNA-SIRGAS 2018, EPSG:9377).

En cuanto a los valores, la pendiente porcentual oscila entre 0.0000% (zonas completamente planas) y 142.6091% (zonas con alta inclinación), lo que permite una interpretación más intuitiva y útil para análisis agrícolas, de conservación o planificación territorial.

(slope_perc = tan(slope*(pi/180))*100)

## class       : SpatRaster 
## dimensions  : 1473, 2441, 1  (nrow, ncol, nlyr)
## resolution  : 152.6319, 152.6319  (x, y)
## extent      : 4533838, 4906412, 1495829, 1720656  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :    slope 
## min value   :   0.0000 
## max value   : 142.6091

Aquí nuevamente se genera un histograma para visualizar valores de pendiente en porcentaje almacenados en el objeto slope_perc.Se observa una marcada concentración de valores en pendientes muy bajas (cercanas a 0%), lo cual indica que la mayor parte del territorio es plano o con muy poca inclinación.A medida que aumenta la pendiente, la frecuencia disminuye drásticamente, confirmando que las pendientes pronunciadas son poco comunes en el área analizada.

terra::hist(slope_perc, 
     main = "Pendiente del Putumayo", 
     xlab = "Pendiente (en porcentaje)")

## Warning: [hist] a sample of 28% of the cells was used (of which 69% was NA)

# 7 Computing zonal statistics (20%) En este punto se realiza un análisis espacial para calcular estadísticas zonales sobre valores de una capa raster (en este caso, el raster de pendiente), delimitadas por zonas específicas (municipios). Este tipo de análisis es útil para conocer características como el promedio de pendiente dentro de cada zona administrativa del Putumayo. La clasificación oficial del IGAC (Instituto Geográfico Agustín Codazzi) para pendientes en Colombia, esta diseñada con fines agrológicos. Esta clasificación divide las pendientes en siete rangos, de acuerdo con su porcentaje (%), desde terrenos ligeramente planos hasta muy escarpados, y asigna un color específico a cada rango (usando valores RGB) para facilitar su visualización cartográfica.

knitr::include_graphics("Z(0,0)_ (3).jpg")

#### 7.1 Reclasificación del raster de pendientes

Para aplicar esta clasificación al raster de pendientes (convertido previamente a porcentaje), se crea una matriz de reclasificación (m) que establece los intervalos y los códigos correspondientes (de 1 a 7). Luego se aplica la función classify() para reclasificar los valores de pendiente de acuerdo con esta matriz.

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)

Una vez reclasificados los datos, se utiliza la función exact_extract() del paquete exactextractr para calcular la pendiente promedio por municipio dentro del departamento del Putumayo, se genera una nueva columna en el objeto espacial munic, donde se almacena el valor promedio de pendiente (en porcentaje) para cada polígono municipal.

(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))

## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |======================================================================| 100%

##  [1] 31.958179 17.830488 13.841450  1.924543  3.003686  2.194637  1.746912
##  [8] 15.264664 37.670464  1.768489 30.336948  2.489117 19.559658

El histograma presenta la distribución de la pendiente promedio (en porcentaje) para los municipios del departamento del Putumayo, calculada a partir del raster de pendiente reclasificado.Se observa que la mayoría de los municipios presentan pendientes promedio bajas, específicamente entre 0 y 10%.

hist(munic$mean_slope, 
     main = "Pendiente promedio Putumayo", 
     xlab = "Pendiente (en percentaje)")

Se utilizó la función exact_extract() del paquete exactextractr para extraer la moda (valor más frecuente) de las clases de pendiente en cada municipio.

(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))

## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)

## 
  |                                                                            
  |                                                                      |   0%
  |                                                                            
  |=====                                                                 |   8%
  |                                                                            
  |===========                                                           |  15%
  |                                                                            
  |================                                                      |  23%
  |                                                                            
  |======================                                                |  31%
  |                                                                            
  |===========================                                           |  38%
  |                                                                            
  |================================                                      |  46%
  |                                                                            
  |======================================                                |  54%
  |                                                                            
  |===========================================                           |  62%
  |                                                                            
  |================================================                      |  69%
  |                                                                            
  |======================================================                |  77%
  |                                                                            
  |===========================================================           |  85%
  |                                                                            
  |=================================================================     |  92%
  |                                                                            
  |======================================================================| 100%

##  [1] 5 4 2 1 1 1 1 4 5 1 5 1 1

Luego, se generó un histograma para visualizar la frecuencia de municipios según la clase de pendiente dominante. El histograma muestra la distribución de los municipios del Putumayo según la clase de pendiente predominante tras la reclasificación. En el eje X se representan las categorías de pendiente (1 a 5), y en el eje Y, la frecuencia de municipios que caen en cada categoría.

La categoría 1 (ligeramente plano, 0-3%) es la más común, con 8 municipios clasificados en este rango, lo que indica que gran parte del territorio tiene relieve plano o muy suave, lo cual es favorable para actividades agrícolas mecanizadas.

Las categorías 2 y 3 (ligeramente inclinado a moderadamente inclinado) tienen frecuencias menores (2 y 1 municipios respectivamente), lo que representa zonas con pendiente suave a moderada.

Las categorías 4 y 5 (fuertemente inclinado y moderadamente escarpado) aparecen con una frecuencia de 2 municipios cada una, evidenciando áreas más pendientes que pueden tener restricciones para el uso agropecuario, especialmente mecanizado.

hist(munic$class, 
     main = "Pendiente reclasificada del Putumayo", 
     xlab = "pendiente (como categoría)")

En esta parte del análisis se realizó la reproyección de los datos raster de pendiente al sistema de coordenadas geográficas WGS 84 (EPSG:4326), utilizando la función project() del paquete terra. Este paso es necesario para asegurar que los datos puedan visualizarse correctamente sobre mapas base que usan coordenadas en grados decimales (latitud y longitud), como los mapas de Leaflet o Google Maps.

Se reproyectaron dos capas:

La capa rc, que representa la pendiente reclasificada en categorías, y se guardó como rc.geo.

(rc.geo = project(rc, "EPSG:4326"))

## class       : SpatRaster 
## dimensions  : 1482, 2437, 1  (nrow, ncol, nlyr)
## resolution  : 0.001373237, 0.001373237  (x, y)
## extent      : -77.18857, -73.84199, -0.5635595, 1.471577  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     slope 
## min value   : 0.3325361 
## max value   : 7.0000000

La capa slope_perc, que representa la pendiente en porcentaje, y se guardó como slope.geo.

Con esta transformación, ambas capas están listas para análisis geoespaciales compatibles con sistemas globales de coordenadas.

(slope.geo = project(slope_perc, "EPSG:4326"))

## class       : SpatRaster 
## dimensions  : 1482, 2437, 1  (nrow, ncol, nlyr)
## resolution  : 0.001373237, 0.001373237  (x, y)
## extent      : -77.18857, -73.84199, -0.5635595, 1.471577  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :        slope 
## min value   : 7.344723e-03 
## max value   : 1.411980e+02

7 8 Plotting terrain slope (20%)

Para visualizar la pendiente del terreno en el departamento del Putumayo, se elaboró un mapa temático interactivo utilizando leaflet. En primer lugar, se definió una paleta de colores (palredgreen) que asigna tonalidades específicas a los valores de pendiente, facilitando así la diferenciación visual entre clases. Se utilizaron colores desde el verde claro (darkseagreen3) para pendientes suaves, hasta el rojo oscuro (darkred) para pendientes más pronunciadas, pasando por tonos intermedios como amarillo (yellow2), naranja (orange) y café (brown2).

palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
  na.color = "transparent")

Luego, con la función leaflet(), se superpusieron los límites municipales (munic) sobre un mapa base (addTiles()), permitiendo una mejor contextualización geográfica. Los polígonos municipales fueron representados con un contorno gris translúcido, y se incluyó un cuadro emergente (popup) que muestra el nombre del municipio y la clase de pendiente predominante.

A través de la función addRasterImage(), se integró la capa raster de pendientes en porcentaje (slope.geo) al mapa, lo que permite observar con la variación espacial de este atributo en todo el departamento. Finalmente, se añadió una leyenda con addLegend(), titulada “Pendiente del terreno en Putumayo (%)”, que facilita la interpretación del mapa.

leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% 
    addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.4, fillOpacity = 0.10,
    popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
                          "Clase de Pendiente: ", munic$class, "<br>")) %>%
  addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8)  %>%
  addLegend(pal = palredgreen, values = values(slope.geo),
    title = "Pendiente del terreno en Putumayo (%)")

## 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'

Click on different sites for getting municipalities’ names and most frequent percentage slope values at each.

8 9 Alternative visualization of terrain slope (15%)

En esta visualización alternativa, se representa la pendiente del terreno del departamento del Putumayo en grados, lo cual permite una mejor comprensión topográfica en contextos más técnicos (por ejemplo, geológicos o de estabilidad de laderas). A diferencia del mapa anterior, este presenta etiquetas con el nombre de cada municipio acompañado de su pendiente media en grados, lo que facilita la lectura rápida del comportamiento general del relieve.

munic$mean_slope_deg <- extract(slope.geo, munic, fun = mean, na.rm = TRUE)[,2]

## Warning: [extract] transforming vector data to the CRS of the raster

centroides <- st_centroid(munic)

## Warning: st_centroid assumes attributes are constant over geometries

coords <- st_coordinates(centroides)

palgreen <- colorNumeric(c("lightgreen", "khaki", "orange", "tomato", "darkred"),
                         values(slope.geo), na.color = "transparent")

leaflet(munic) %>%
  addTiles() %>%
  setView(lng = -75.9, lat = 8.7, zoom = 9) %>%
  addPolygons(color = "gray", weight = 1, smoothFactor = 0.5,
              opacity = 0.4, fillOpacity = 0.1,
              popup = paste("Municipio: ", munic$mpio_cnmb, "<br>",
                            "Slope class: ", munic$class, "<br>")) %>%
  addRasterImage(slope.geo, colors = palgreen, opacity = 0.8) %>%
  addLegend(pal = palgreen, values = values(slope.geo),
            title = "Pendiente del terreno (°)") %>%
  addLabelOnlyMarkers(
    lng = coords[,1],
    lat = coords[,2],
    label = paste(munic$mpio_cnmb, ": ", round(munic$mean_slope_deg, 1), "°"),
    labelOptions = labelOptions(noHide = TRUE, direction = "top", textOnly = TRUE)
  )

## 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'

9 10 Visualizing terrain aspect (15%)

En este bloque se calculó el aspecto del terreno, es decir, la orientación de las pendientes en grados, usando el modelo digital de elevación (dem). Luego, se redujo la resolución del raster con aggregate() para evitar errores al visualizarlo en un mapa interactivo con leaflet. Finalmente, se utilizó una paleta circular de colores tipo brújula para representar las distintas orientaciones (norte, sur, este, oeste y puntos intermedios), junto con una leyenda que indica los valores en grados.Permite identificar hacia dónde están orientadas las laderas en el departamento del Putumayo, lo cual es clave para la gestión del agua, la erosión del suelo y el uso agrícola, ya que ciertos cultivos se desarrollan mejor en determinadas orientaciones por la cantidad de luz solar que reciben.

# Cálculo del aspecto del terreno (orientación de las pendientes)
aspect <- terrain(dem, v = "aspect", unit = "degrees")

# Reducimos la resolución del raster para que sea compatible con leaflet
aspect_small <- terra::aggregate(aspect, fact = 5)

# Definimos una paleta circular tipo brújula para representar los aspectos
pal_aspect <- colorNumeric(
  palette = c("blue", "cyan", "green", "yellow", "orange", "red", "magenta", "blue"),
  domain = values(aspect_small),
  na.color = "transparent"
)

# Mapa interactivo con leaflet para visualizar el aspecto del terreno
leaflet() %>%
  addTiles() %>%
  addRasterImage(aspect_small, colors = pal_aspect, opacity = 0.8) %>%
  addLegend(
    pal = pal_aspect,
    values = values(aspect_small),
    title = "Aspecto del terreno (°)",
    labFormat = labelFormat(suffix = "°")
  )

10 Bibliography

Lizarazo, I., 2025. Atributos del terreno geomorfométrico en R. Disponible en https://rpubs.com/ials2un/geomorphométrico

sessionInfo()

## R version 4.3.3 (2024-02-29 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.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_0.9.1  terra_1.8-29        
## [4] leaflet_2.2.2        sf_1.0-20            elevatr_0.99.0      
## 
## loaded via a namespace (and not attached):
##  [1] s2_1.1.7           utf8_1.2.4         sass_0.4.8         generics_0.1.3    
##  [5] class_7.3-22       KernSmooth_2.23-22 jpeg_0.1-10        lattice_0.22-6    
##  [9] digest_0.6.34      magrittr_2.0.3     rgl_1.3.1          evaluate_0.23     
## [13] grid_4.3.3         fastmap_1.2.0      jsonlite_1.8.8     e1071_1.7-16      
## [17] DBI_1.2.3          promises_1.3.2     purrr_1.0.2        fansi_1.0.6       
## [21] scales_1.3.0       crosstalk_1.2.1    codetools_0.2-19   jquerylib_0.1.4   
## [25] cli_3.6.2          shiny_1.10.0       rlang_1.1.3        units_0.8-7       
## [29] ellipsis_0.3.2     munsell_0.5.0      base64enc_0.1-3    cachem_1.1.0      
## [33] yaml_2.3.8         raster_3.6-32      tools_4.3.3        dplyr_1.1.4       
## [37] colorspace_2.1-0   httpuv_1.6.14      png_0.1-8          vctrs_0.6.5       
## [41] R6_2.5.1           mime_0.12          proxy_0.4-27       lifecycle_1.0.4   
## [45] classInt_0.4-11    htmlwidgets_1.6.4  pkgconfig_2.0.3    progressr_0.15.1  
## [49] bslib_0.6.1        pillar_1.9.0       later_1.3.2        glue_1.7.0        
## [53] Rcpp_1.0.12        highr_0.10         tidyselect_1.2.1   tibble_3.2.1      
## [57] xfun_0.42          rstudioapi_0.17.1  knitr_1.45         farver_2.1.1      
## [61] xtable_1.8-4       htmltools_0.5.7    rmarkdown_2.25     wk_0.9.4          
## [65] compiler_4.3.3     sp_2.2-0