** 📚 Tabla de contenido **

1. 🧭 Introducción

2. ⚙️ Configuración

3. 📥Introducción a MultiscaleDTM

4. 🧾Lectura de datos de entrada

5. 🧭Transformación de coordenadas

6. 📈Cálculo de atributos del terreno

7. 📈Cálculo de estadísticas zonales

8. 📐Trazado de la pendiente del terreno

9. ✅ Visualización alternativa de la pendiente del terreno

10. ✅ Visualización de la orientación del terreno

11. 📚Bibliografía

1. Introducción 🧭

Este cuaderno tiene como objetivo ilustrar el cálculo de atributos geomorfométricos del terreno a partir de Modelos Digitales de Elevación (MDE) cuadriculados, utilizando el paquete MultiscaleDTM en R. Esta herramienta permite derivar indicadores esenciales como pendiente, orientación (aspect), curvatura, posición relativa y rugosidad, aplicándolos a múltiples escalas espaciales (Ilich et al., 2023). Los atributos geomorfométricos son fundamentales en diversas disciplinas:son clave para entender la hidrología del terreno, la erosión, y la distribución espacial de propiedades del suelo (Franklin, 2020; Kokulan et al., 2018). En contextos agrícolas, la topografía—especialmente la pendiente y la orientación—afecta el flujo de agua, la exposición solar y la variabilidad de nutrientes del suelo, impactando directamente el rendimiento de los cultivos (Geopard Tech, 2020; 2022). El área de estudio será el departamento del Tolima, seleccionado por su relevancia agrícola y su diversidad de relieves, lo que lo convierte en un escenario ideal para ilustrar la relación entre geomorfometría y producción agropecuaria.

2. Configuración ⚙️

Los paquetes necesarios deben instalarse previamente desde la consola. Paquetes ncesarios = c(“MultiscaleDTM”, “exactextractr”) con la funcion install.packages(paquetes)

Se debe empezar por limpiar la memoria: Esto se hace para comenzar un script desde cero, para evitar que datos anteriores interfieran y asegurar que el código funcione con un entorno limpio (por ejemplo, antes de compartirlo o reproducirlo).

rm(list=ls())

Ahora, se cargan las bibliotecas necesarias:

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.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE
library(leaflet)
library(terra)
## terra 1.8.54
library(MultiscaleDTM)
library(exactextractr)

package_name::function.

3. Introducción a MultiscaleDTM 📥

El paquete MultiscaleDTM en R permite calcular atributos geomorfométricos multiescala—como pendiente, orientación, curvatura, posición relativa y rugosidad—utilizando modelos digitales del terreno (MDT), es decir, rásteres de elevación con cuadrícula regular, aplicando una ventana focal definida (Ilich et al., 2023). sta herramienta es particularmente útil porque permite evaluar las propiedades del terreno a diferentes escalas espaciales, una capacidad que muchos paquetes no ofrecen (Ilich et al., 2023). El presente cuaderno tiene por objetivo introducir brevemente el uso de MultiscaleDTM, orientado a estudiantes de pregrado en agronomía de la Universidad Nacional de Colombia. Se aplicará específicamente al departamento del Tolima, seleccionado por su riqueza topográfica y relevancia agrícola, lo que lo convierte en un escenario adecuado para ejemplificar cómo estas variables geomorfológicas influyen en la dinámica del agua, la exposición solar y, potencialmente, la productividad del suelo.

En lo que sigue: - Se instalará y cargará el paquete MultiscaleDTM. - Se explicará cada fragmento de código. - Se describirá e interpretará cada resultado generado. - Se evaluará la multiescalaridad de los atributos en el contexto del Tolima.

4. Lectura de datos de entrada 🧾

Primero, se carga un DEM ráster para el departamento en estudio(Tolima). Para luego leer el DEM con terra:

(dem = terra::rast("C:/Users/menju/OneDrive/Documentos/GB2/Rstudio2/elev_ELEVACION_z10.tif"))
## class       : SpatRaster 
## size        : 4092, 3078, 1  (nrow, ncol, nlyr)
## resolution  : 0.000685414, 0.000685414  (x, y)
## extent      : -76.28906, -74.17936, 2.811272, 5.615986  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs 
## source      : elev_ELEVACION_z10.tif 
## name        : elev_ELEVACION_z10 
## min value   :               -427 
## max value   :               5372

Reduzcamos la resolución de DEM para evitar problemas de memoria: Para reducir la resolución del DEM (hacerlo más “grueso”) y disminuir el uso de memoria y procesamiento al trabajar con grandes volúmenes de datos.

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

Se leen los municipios del Tolima utilizando la biblioteca sf:

(munic<-  sf::st_read("C:/Users/menju/OneDrive/Documentos/GB2/P1/DATOS/MUNI_TOL.gpkg"))
## Reading layer `MUNICIPIOS_TOL' from data source 
##   `C:\Users\menju\OneDrive\Documentos\GB2\P1\DATOS\MUNI_TOL.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS:  MAGNA-SIRGAS
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS:  MAGNA-SIRGAS
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr        mpio_cnmbr
## 1          73        001      73001     TOLIMA            IBAGUÉ
## 2          73        024      73024     TOLIMA         ALPUJARRA
## 3          73        026      73026     TOLIMA          ALVARADO
## 4          73        030      73030     TOLIMA          AMBALEMA
## 5          73        043      73043     TOLIMA        ANZOÁTEGUI
## 6          73        055      73055     TOLIMA            ARMERO
## 7          73        067      73067     TOLIMA             ATACO
## 8          73        124      73124     TOLIMA         CAJAMARCA
## 9          73        148      73148     TOLIMA CARMEN DE APICALÁ
## 10         73        152      73152     TOLIMA        CASABIANCA
##                               mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1                                   1906 MUNICIPIO  1377.2514      2023
## 2  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO   501.5113      2023
## 3                                   1540 MUNICIPIO   344.3356      2023
## 4                                   1627 MUNICIPIO   238.0376      2023
## 5   Ordenanza 21 del 30 de Marzo de 1915 MUNICIPIO   469.7113      2023
## 6  Decreto 1049 de Septiembre 29 de 1908 MUNICIPIO   440.6308      2023
## 7  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO  1017.6598      2023
## 8                                   1913 MUNICIPIO   508.4909      2023
## 9  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO   190.8361      2023
## 10  Ordenanza 26 del 22 de Junio de 1896 MUNICIPIO   174.9306      2023
##    shape_Leng shape_Area codigo                           geom
## 1   2.4800420 0.11217106     NA MULTIPOLYGON (((-75.36517 4...
## 2   1.0547137 0.04080342     NA MULTIPOLYGON (((-74.99192 3...
## 3   1.0695888 0.02805438     NA MULTIPOLYGON (((-74.97668 4...
## 4   1.0091029 0.01940220     NA MULTIPOLYGON (((-74.74744 4...
## 5   1.3101775 0.03826692     NA MULTIPOLYGON (((-75.20098 4...
## 6   1.4441004 0.03592418     NA MULTIPOLYGON (((-74.88687 5...
## 7   2.2522480 0.08276960     NA MULTIPOLYGON (((-75.3251 3....
## 8   1.0293681 0.04140402     NA MULTIPOLYGON (((-75.48219 4...
## 9   0.7339295 0.01554073     NA MULTIPOLYGON (((-74.74471 4...
## 10  1.0803200 0.01425931     NA MULTIPOLYGON (((-75.07491 5...

Se debe tener en cuenta que ambos conjuntos de datos están en coordenadas geográficas.

Ahora, se recorta el objeto dem2 a los límites de nuestro departamento:

# 1. Convertir sf a SpatVector
munic_vect <- terra::vect(munic)
# 2. Reproyectar al CRS del raster
munic_vect <- terra::project(munic_vect, dem2)
# 3. Recortar el DEM usando el objeto reproyectado
dem3 <- terra::crop(dem2, munic_vect, mask = TRUE)

5. Transformación de coordenadas 🧭

El Instituto Geográfico Agustín Codazzi (IGAC) adoptó en el año 2020 un nuevo sistema de coordenadas planas de origen único nacional para todo el territorio colombiano. El código EPSG correspondiente a este sistema es el EPSG:9377, conocido como MAGNA-SIRGAS / Colombia Bogotá.

Este sistema de referencia se fundamenta en la proyección Transversa de Mercator con origen en Bogotá D.C., y pertenece al marco geodésico MAGNA-SIRGAS, el cual es parte de la Red Geodésica Nacional de Colombia. Esta adopción responde a la necesidad de contar con un sistema de coordenadas unificado que facilite la interoperabilidad geoespacial, elimine la fragmentación entre zonas UTM y mejore la precisión en el análisis de información geográfica, especialmente en modelos digitales del terreno, cartografía base y catastros.

Importancia del uso del EPSG:9377 en análisis de terreno El cálculo de atributos topográficos como pendiente, orientación (aspecto), curvatura y flujo superficial requiere que el Modelo Digital de Elevación (MDE) esté en un sistema de coordenadas proyectadas, es decir, en unidades métricas (como metros), no geográficas (grados). Esto es crucial, ya que las operaciones de derivación espacial (como gradientes o transformaciones de superficie) dependen de distancias y áreas reales.

Trabajar con EPSG:9377 garantiza resultados consistentes a lo largo del territorio colombiano sin necesidad de subdividir el país por zonas UTM, como se hacía anteriormente. Además, favorece la integración con sistemas oficiales como el Catastro Multipropósito, la Infraestructura Colombiana de Datos Espaciales (ICDE), y los datos de planificación territorial y ambiental.

# check the parameters of this terra function here
# https://rdrr.io/github/rspatial/terra/man/project.html
(dem_plane = project(dem3, "EPSG:9377"))
## class       : SpatRaster 
## size        : 1787, 1197, 1  (nrow, ncol, nlyr)
## resolution  : 151.7803, 151.7803  (x, y)
## extent      : 4654850, 4836531, 1875472, 2146704  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        : elev_ELEVACION_z10 
## min value   :          -1.783039 
## max value   :        5277.944824

También es necesario realizar una transformación similar para el objeto vectorial:

# check the parameters of this sf function here
# https://github.com/rstudio/cheatsheets/blob/main/sf.pdf
(munic_plane = sf::st_transform(munic, "EPSG:9377"))
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4654882 ymin: 1875660 xmax: 4836359 ymax: 2146078
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr        mpio_cnmbr
## 1          73        001      73001     TOLIMA            IBAGUÉ
## 2          73        024      73024     TOLIMA         ALPUJARRA
## 3          73        026      73026     TOLIMA          ALVARADO
## 4          73        030      73030     TOLIMA          AMBALEMA
## 5          73        043      73043     TOLIMA        ANZOÁTEGUI
## 6          73        055      73055     TOLIMA            ARMERO
## 7          73        067      73067     TOLIMA             ATACO
## 8          73        124      73124     TOLIMA         CAJAMARCA
## 9          73        148      73148     TOLIMA CARMEN DE APICALÁ
## 10         73        152      73152     TOLIMA        CASABIANCA
##                               mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1                                   1906 MUNICIPIO  1377.2514      2023
## 2  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO   501.5113      2023
## 3                                   1540 MUNICIPIO   344.3356      2023
## 4                                   1627 MUNICIPIO   238.0376      2023
## 5   Ordenanza 21 del 30 de Marzo de 1915 MUNICIPIO   469.7113      2023
## 6  Decreto 1049 de Septiembre 29 de 1908 MUNICIPIO   440.6308      2023
## 7  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO  1017.6598      2023
## 8                                   1913 MUNICIPIO   508.4909      2023
## 9  Decreto 650 del 13 de Octubre de 1887 MUNICIPIO   190.8361      2023
## 10  Ordenanza 26 del 22 de Junio de 1896 MUNICIPIO   174.9306      2023
##    shape_Leng shape_Area codigo                           geom
## 1   2.4800420 0.11217106     NA MULTIPOLYGON (((4737717 207...
## 2   1.0547137 0.04080342     NA MULTIPOLYGON (((4778805 194...
## 3   1.0695888 0.02805438     NA MULTIPOLYGON (((4780825 207...
## 4   1.0091029 0.01940220     NA MULTIPOLYGON (((4806322 210...
## 5   1.3101775 0.03826692     NA MULTIPOLYGON (((4755947 207...
## 6   1.4441004 0.03592418     NA MULTIPOLYGON (((4790926 212...
## 7   2.2522480 0.08276960     NA MULTIPOLYGON (((4741831 196...
## 8   1.0293681 0.04140402     NA MULTIPOLYGON (((4724687 206...
## 9   0.7339295 0.01554073     NA MULTIPOLYGON (((4806435 202...
## 10  1.0803200 0.01425931     NA MULTIPOLYGON (((4770073 212...

6. Cálculo de atributos del terreno 📈

Se procede a calcular la pendiente y la orientación.

# Explain what is w 
# Explain what is method
# Change if needed
(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        : 1787, 1197, 2  (nrow, ncol, nlyr)
## resolution  : 151.7803, 151.7803  (x, y)
## extent      : 4654850, 4836531, 1875472, 2146704  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## names       :   slope,       aspect 
## min values  :  0.0000, 7.917609e-04 
## max values  : 69.1466, 3.599997e+02

Ahora, dividamos el objeto slp_asp en dos partes.

Capa 1:

(slope = subset(slp_asp, 1))
## class       : SpatRaster 
## size        : 1787, 1197, 1  (nrow, ncol, nlyr)
## resolution  : 151.7803, 151.7803  (x, y)
## extent      : 4654850, 4836531, 1875472, 2146704  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :   slope 
## min value   :  0.0000 
## max value   : 69.1466

¿Cuál es la distribución de valores de pendiente?

terra::hist(slope, 
     main = "Tolima's slope", 
     xlab = "Slope (in degrees)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)

El histograma mostrado representa la distribución de la pendiente del terreno (en grados) en el departamento del Tolima, Colombia. Se observa una clara asimetría hacia la derecha, lo que indica que la mayoría del área analizada presenta pendientes suaves. Específicamente, la mayor concentración de valores se encuentra entre los 0° y 10°, con una frecuencia que supera los 120,000 píxeles, lo cual sugiere que una gran proporción del terreno es relativamente plano o ligeramente inclinado.

A medida que aumenta el valor de la pendiente, la frecuencia disminuye de forma progresiva. Pendientes mayores a 30° son considerablemente menos frecuentes, y los valores extremos (mayores a 50°) son casi inexistentes, lo que confirma que las áreas abruptas son escasas en la región.

Este patrón es típico en territorios con una combinación de zonas planas agrícolas (como el valle del Magdalena) y zonas montañosas en la cordillera central, reflejando la heterogeneidad del relieve tolimense. Este tipo de análisis es crucial en estudios de aptitud agrícola, riesgo de erosión y ordenamiento territorial, ya que la pendiente influye directamente en la capacidad de uso del suelo y en procesos como escorrentía o retención de agua.

Capa 2:

(aspect =subset(slp_asp, 2))
## class       : SpatRaster 
## size        : 1787, 1197, 1  (nrow, ncol, nlyr)
## resolution  : 151.7803, 151.7803  (x, y)
## extent      : 4654850, 4836531, 1875472, 2146704  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :       aspect 
## min value   : 7.917609e-04 
## max value   : 3.599997e+02
terra::hist(aspect, 
     main = "Tolima aspect", 
     xlab = "Aspect (in degrees)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)

El histograma presentado muestra la distribución del aspecto del terreno (orientación de las pendientes) en el departamento del Tolima, expresado en grados desde 0° (Norte) hasta 360° (también Norte, cerrado el círculo). Este atributo geomorfométrico es crucial para comprender la exposición solar, la circulación de aire, y el microclima local, todos factores determinantes en procesos ecológicos, agrícolas y de erosión.

En la gráfica se observa una distribución relativamente equilibrada, aunque con una mayor frecuencia de aspectos entre los 90° y 150°, lo cual corresponde a pendientes orientadas hacia el Este y Sureste. Este patrón puede deberse a la influencia de la cordillera Central y la topografía del valle del río Magdalena, donde las formas del terreno tienden a inclinarse en esa dirección.

Se destaca una ligera disminución de frecuencias entre los 180° y 270°, es decir, en las orientaciones Sur a Oeste, posiblemente porque esas laderas son menos predominantes en las unidades de relieve analizadas. Por otra parte, existe una recuperación hacia el rango de 270° a 360°, sugiriendo una mayor presencia de laderas con exposición Noroeste.

Este tipo de análisis es útil para:

  • Planificación agrícola, ya que las orientaciones determinan la radiación solar disponible.

  • Gestión de cuencas hidrográficas, pues el aspecto influye en el derretimiento de nieve (en zonas de alta montaña) o en la dirección del escurrimiento superficial.

  • Ecología del paisaje, al afectar la distribución de la vegetación, humedad del suelo y condiciones microclimáticas.

Además, se convierten los grados de pendiente en porcentaje de pendiente:

(slope_perc = tan(slope*(pi/180))*100)
## class       : SpatRaster 
## size        : 1787, 1197, 1  (nrow, ncol, nlyr)
## resolution  : 151.7803, 151.7803  (x, y)
## extent      : 4654850, 4836531, 1875472, 2146704  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377) 
## source(s)   : memory
## name        :    slope 
## min value   :   0.0000 
## max value   : 262.5146
terra::hist(slope_perc, 
     main = "Tolima's slope", 
     xlab = "Slope (in percentage)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)

El histograma presentado muestra la distribución de la pendiente del terreno en porcentaje (%) para el departamento del Tolima, calculada a partir de un Modelo Digital de Elevación (MDE). Esta variable es clave para analizar procesos erosivos, escorrentía, accesibilidad del terreno y aptitud agrícola.

Se observa una distribución fuertemente sesgada hacia la izquierda, con la mayor parte de las pendientes por debajo del 30%. Las clases de pendiente con valores bajos (0–20%) son las más frecuentes, superando ampliamente las 120.000 observaciones en los primeros intervalos. Esto indica que una proporción significativa del territorio del Tolima es relativamente plano o con pendiente suave, lo cual es favorable para actividades agrícolas y de infraestructura.

A medida que aumentan los valores de pendiente, la frecuencia disminuye rápidamente. Las pendientes mayores al 100% (equivalente a un ángulo de inclinación de 45°) son escasas, y representan zonas de topografía abrupta, posiblemente correspondientes a áreas de montaña, como las zonas cercanas al Nevado del Tolima o la cordillera Central. Esta distribución sugiere que el uso del suelo y la ocupación humana en el Tolima están favorecidos por las amplias zonas de baja pendiente.mSin embargo, las zonas de alta pendiente, aunque menos frecuentes, deben ser cuidadosamente gestionadas para evitar riesgos de erosión, deslizamientos y para garantizar prácticas agrícolas sostenibles.

7. Cálculo de estadísticas zonales 📈

Una operación de estadísticas zonales calcula estadísticas sobre los valores de celda de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En este ncaso, nos interesa calcular el promedio de los valores de pendiente por municipio.

Primero, se conoce la clasificación de pendientes del IGAC para fines agrológicos:

(Reclacificar un Ráster)

# Reclassify the slope 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)

Ahora, se calculan 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%  |                                                                              |=                                                                     |   2%  |                                                                              |===                                                                   |   4%  |                                                                              |====                                                                  |   6%  |                                                                              |======                                                                |   9%  |                                                                              |=======                                                               |  11%  |                                                                              |=========                                                             |  13%  |                                                                              |==========                                                            |  15%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  21%  |                                                                              |================                                                      |  23%  |                                                                              |==================                                                    |  26%  |                                                                              |===================                                                   |  28%  |                                                                              |=====================                                                 |  30%  |                                                                              |======================                                                |  32%  |                                                                              |========================                                              |  34%  |                                                                              |=========================                                             |  36%  |                                                                              |===========================                                           |  38%  |                                                                              |============================                                          |  40%  |                                                                              |==============================                                        |  43%  |                                                                              |===============================                                       |  45%  |                                                                              |=================================                                     |  47%  |                                                                              |==================================                                    |  49%  |                                                                              |====================================                                  |  51%  |                                                                              |=====================================                                 |  53%  |                                                                              |=======================================                               |  55%  |                                                                              |========================================                              |  57%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  62%  |                                                                              |=============================================                         |  64%  |                                                                              |==============================================                        |  66%  |                                                                              |================================================                      |  68%  |                                                                              |=================================================                     |  70%  |                                                                              |===================================================                   |  72%  |                                                                              |====================================================                  |  74%  |                                                                              |======================================================                |  77%  |                                                                              |=======================================================               |  79%  |                                                                              |=========================================================             |  81%  |                                                                              |==========================================================            |  83%  |                                                                              |============================================================          |  85%  |                                                                              |=============================================================         |  87%  |                                                                              |===============================================================       |  89%  |                                                                              |================================================================      |  91%  |                                                                              |==================================================================    |  94%  |                                                                              |===================================================================   |  96%  |                                                                              |===================================================================== |  98%  |                                                                              |======================================================================| 100%
##  [1] 28.327206 22.369404 15.429619  3.669801 33.040310  9.219158 27.107445
##  [8] 36.053734 15.664359 35.969097 28.613630 18.607035  7.297594 18.536375
## [15] 20.129711  1.252090 18.143921  1.507208 24.973742  2.159666 34.832108
## [22] 16.668184 22.115511 11.497709 32.145004 11.752955 16.978456 27.755978
## [29] 16.803055 18.587324 20.755070  8.394407 34.529945 15.837931 10.465549
## [36] 41.639145 27.062590 30.451038  1.194180 36.182438 11.381176 31.772017
## [43] 13.304363 17.442604 11.510288 35.072704 18.109522
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
hist(munic$mean_slope, 
     main = "Tolima's mean slope", 
     xlab = "Slope (in percentage)")

El histograma muestra la pendiente promedio (en porcentaje) de los municipios del Tolima. Se observa que la mayoría de los municipios tienen pendientes medias entre 15% y 25%, con un pico en el intervalo de 15–20%. Solo unos pocos municipios presentan pendientes promedio superiores al 35%, lo que sugiere que en general el relieve municipal es moderadamente inclinado.

(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
##   |                                                                              |                                                                      |   0%  |                                                                              |=                                                                     |   2%  |                                                                              |===                                                                   |   4%  |                                                                              |====                                                                  |   6%  |                                                                              |======                                                                |   9%  |                                                                              |=======                                                               |  11%  |                                                                              |=========                                                             |  13%  |                                                                              |==========                                                            |  15%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  21%  |                                                                              |================                                                      |  23%  |                                                                              |==================                                                    |  26%  |                                                                              |===================                                                   |  28%  |                                                                              |=====================                                                 |  30%  |                                                                              |======================                                                |  32%  |                                                                              |========================                                              |  34%  |                                                                              |=========================                                             |  36%  |                                                                              |===========================                                           |  38%  |                                                                              |============================                                          |  40%  |                                                                              |==============================                                        |  43%  |                                                                              |===============================                                       |  45%  |                                                                              |=================================                                     |  47%  |                                                                              |==================================                                    |  49%  |                                                                              |====================================                                  |  51%  |                                                                              |=====================================                                 |  53%  |                                                                              |=======================================                               |  55%  |                                                                              |========================================                              |  57%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  62%  |                                                                              |=============================================                         |  64%  |                                                                              |==============================================                        |  66%  |                                                                              |================================================                      |  68%  |                                                                              |=================================================                     |  70%  |                                                                              |===================================================                   |  72%  |                                                                              |====================================================                  |  74%  |                                                                              |======================================================                |  77%  |                                                                              |=======================================================               |  79%  |                                                                              |=========================================================             |  81%  |                                                                              |==========================================================            |  83%  |                                                                              |============================================================          |  85%  |                                                                              |=============================================================         |  87%  |                                                                              |===============================================================       |  89%  |                                                                              |================================================================      |  91%  |                                                                              |==================================================================    |  94%  |                                                                              |===================================================================   |  96%  |                                                                              |===================================================================== |  98%  |                                                                              |======================================================================| 100%
##  [1] 5 4 1 1 5 1 5 5 1 5 5 4 1 4 4 1 4 1 5 1 5 4 4 1 5 4 4 5 1 5 4 1 5 4 1 5 4 5
## [39] 1 5 1 5 1 4 1 5 4
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
hist(munic$class, 
     main = "Tolima's reclassified slope", 
     xlab = "Slope (as a category)")

El histograma muestra la pendiente reclasificada de los municipios del Tolima en cinco categorías. Las clases 1 y 5 (pendientes muy bajas y muy altas, respectivamente) son las más frecuentes, cada una con alrededor de 16 municipios. Las clases intermedias presentan menor frecuencia, destacándose la clase 2, que no tiene representación. Esto evidencia una distribución bimodal del relieve en el departamento.

Transformar nuevamente la pendiente a coordenadas geográficas:

# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))
## class       : SpatRaster 
## size        : 1793, 1199, 1  (nrow, ncol, nlyr)
## resolution  : 0.001370805, 0.001370805  (x, y)
## extent      : -76.11484, -74.47124, 2.868096, 5.325948  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :       slope 
## min value   :   0.2429543 
## max value   : 259.9135742
#This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))
## class       : SpatRaster 
## size        : 1793, 1199, 1  (nrow, ncol, nlyr)
## resolution  : 0.001370805, 0.001370805  (x, y)
## extent      : -76.11484, -74.47124, 2.868096, 5.325948  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :        slope 
## min value   : 1.254625e-04 
## max value   : 2.599136e+02

8. Trazado de la pendiente del terreno 📐

Primero, se seleciona una paleta de colores para la pendiente:

# expand the number of colors to improve your visualization
# https://r-graph-gallery.com/42-colors-names.html
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
  na.color = "transparent")

Ahora es el momento de trazar el gráfico. Asegúrate de cambiar varias opciones en el código si es necesario:

leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>% 
    addPolygons(color = "black", 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 = "Terrain slope in Tolima (%)")
## 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'

Es posible hacer clic en diferentes sitios para obtener los nombres de los municipios y los valores porcentuales de pendiente más frecuentes en cada uno.

9. Visualización alternativa de la pendiente del terreno

Visualización de la pendiente del terreno en grados respecto a los nombres de los municipios junto con sus valores de pendiente media correspondientes.

summary(values(slope.geo))  # o
##      slope        
##  Min.   :  0.00   
##  1st Qu.:  9.23   
##  Median : 21.46   
##  Mean   : 23.96   
##  3rd Qu.: 36.01   
##  Max.   :259.91   
##  NA's   :1103290
hist(values(slope.geo), main = "Distribucion de pendiente (grados)")

# Paleta para pendiente en grados
palSlopeDeg <- colorNumeric(
  palette = c("lightblue", "yellow", "orange", "red", "darkred"),
  domain = values(slope), na.color = "transparent"
)


leaflet(munic) %>%
  addTiles() %>%
  addRasterImage(slope.geo, colors = palSlopeDeg, opacity = 0.8) %>%
  addPolygons(color = "black", weight = 1.0, fill = FALSE,
              label = ~paste0(mpio_cnmbr, ": ", round(slope.geo, 1), "%"),
              popup = ~paste("<strong>Municipio: </strong>", mpio_cnmbr, "<br>",
                             "<strong>Pendiente media: </strong>", round(slope.geo, 1), "%")) %>%
  addLegend(pal = palSlopeDeg, values = values(slope.geo),
            title = "Slope (degrees)", position = "bottomright")
## 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'

El mapa muestra la pendiente del terreno en el Tolima en grados. Se observa que las zonas con pendientes más altas (tonos rojos intensos) se concentran en el occidente y sur del departamento, asociadas con áreas montañosas de la cordillera Central. En contraste, las regiones del centro-este presentan pendientes más suaves (tonos amarillos y claros), lo cual sugiere un relieve más plano, típicamente asociado con zonas agrícolas y de mayor accesibilidad. Esta variabilidad topográfica es clave para planificar el uso del suelo y actividades agropecuarias.

10. Visualización de la orientación del terreno

Visualiza la orientación del terreno con una paleta similar a esta.

# Guardar en un nombre diferente para evitar el conflicto
writeRaster(slope.geo, "elev_ELEVACION_z10.tif", overwrite = TRUE)

# Volver a cargar el raster limpio
aspect_raster <- rast("elev_ELEVACION_z10.tif")
names(aspect_raster) <- "aspect"
Sys.setlocale("LC_ALL", "es_ES.UTF-8")
## [1] "LC_COLLATE=es_ES.UTF-8;LC_CTYPE=es_ES.UTF-8;LC_MONETARY=es_ES.UTF-8;LC_NUMERIC=C;LC_TIME=es_ES.UTF-8"

Clasificar el ráster aspect en clases cardinales

Acontinuacion se muestra el paso a paso para obtener el mapa de orientacion de terreno, se representara en forma de texto e imagen “pantallazo” debido a que el mapa representa una ocupación de memoria extremadamente alta impidiendo la publicación en Rpubs, se tomó esta desición despues de multiples intentos fallidos.

#library(leaflet)
#library(terra)
#library(sf)
#library(dplyr)

# ---------------------------------------------------
# 1. CARGAR Y PREPARAR DATOS RASTER (ASPECTO)
# ---------------------------------------------------

# Cargar raster de aspecto
#aspect_raster <- rast("elev_ELEVACION_z10.tif")

# Verificar y reproyectar si es necesario
#if (!grepl("4326", crs(aspect_raster))) {
#aspect_raster <- project(aspect_raster, "EPSG:4326")
#}

# Clasificar el aspecto en 8 direcciones cardinales
#breaks <- c(0, 22.5, 67.5, 112.5, 157.5, 202.5, 247.5, 292.5, 337.5, 360)
#aspect_class <- classify(aspect_raster, breaks, include.lowest = TRUE)

# ---------------------------------------------------
# 2. CARGAR DATOS VECTORIALES (MUNICIPIOS Y CIUDADES)
# ---------------------------------------------------

# Opción A: Cargar desde archivos shapefile (reemplaza con tus rutas)
#municipios <- st_read("C:/Users/menju/OneDrive/Documentos/GB2/Rstudio4-Informe/MGN_ADM_MPIO_GRAFICO.shp")

# Filtrar solo los municipios del Tolima (código DANE 73)
#tolima <- municipios %>% 
 # filter(mpio_cnmbr == "Tolima") %>%  # 73 es el código del Tolima
 # st_transform(4326)  # Asegurar CRS WGS84

# Crear datos de ciudades principales del Tolima
#ciudades_tolima <- data.frame(
#  nombre = c("Ibagué", "Espinal", "Honda", "Mariquita", "Líbano"),
  #lon = c(-75.232, -74.884, -74.736, -74.893, -75.062),
  #lat = c(4.445, 4.149, 5.208, 5.199, 4.921)
#) %>% 
 # st_as_sf(coords = c("lon", "lat"), crs = 4326)


# ---------------------------------------------------
# 3. DEFINIR PALETA Y ETIQUETAS
# ---------------------------------------------------

#etiquetas <- c("Norte", "Noreste", "Este", "Sureste", 
 #              "Sur", "Suroeste", "Oeste", "Noroeste")

#paleta_aspecto <- colorFactor(
  #palette = c("#FF0000", "#FFA500", "#FFFF00", "#00FF00",
  #           "#00FFFF", "#0000FF", "#4B0082", "#800080"),
 # domain = 1:8,
  #na.color = "transparent"
#)

# ---------------------------------------------------
# 4. CONSTRUIR EL MAPA INTERACTIVO
# ---------------------------------------------------

#mapa_final <- leaflet() %>%
  # Capa base
  #addTiles(group = "OpenStreetMap") %>%
  #addProviderTiles(providers$CartoDB.Positron, group = "CartoDB") %>%
  
  # Centrar en la región de interés
  #setView(lng = -75.9, lat = 4.8, zoom = 8) %>%
  
  # Capa de aspecto (orientación)
  #addRasterImage(
  #  aspect_class,
  #  colors = paleta_aspecto,
  #  opacity = 0.8,
  #  group = "Orientación",
  #  project = FALSE
  #) %>%
  
  # Leyenda para la orientación
  #addLegend(
  #  pal = paleta_aspecto,
  #  values = 1:8,
  #  title = "Orientación del Terreno",
  #  labFormat = labelFormat(
  #    suffix = "",
  #    transform = function(x) etiquetas[as.numeric(x)]
  #  ),
  #  position = "bottomright",
  #  group = "Leyenda"
  #) %>%
  
  # Barra de escala
  #addScaleBar(position = "bottomleft") %>%
  
  # Capa de municipios (si existe)
 # {
 #   if(exists("municipios")) {
 #     (.) %>% addPolygons(
 #       data = municipios,
 #       color = "#444444",
 #       weight = 1,
 #       fill = FALSE,
 #       group = "MGN_ADM_MPIO_GRAFICO",
 #     
 #     )
 #   } else {
 #     (.) %>% addControl(
 #       "Datos municipales no disponibles",
 #       position = "topright"
 #     )
 #   }
 # } %>%
  
  # Capa de ciudades (si existe)
 # {
 #   if(exists("ciudades")) {
 #     (.) %>% addCircleMarkers(
 #       data = ciudades,
 #       radius = 5,
 #       color = "black",
 #       fillColor = "white",
 #       fillOpacity = 1,
 #       stroke = TRUE,
 #       weight = 1,
 #       group = "Ciudades",
 #       label = ~nombre,
 #       labelOptions = labelOptions(
 #         noHide = FALSE,
 #         direction = "auto",
 #         textOnly = FALSE
 #       )
 #     )
 #   } else {
 #     (.)
 #   }
 # } %>%
  
  # Control de capas
  #addLayersControl(
  #  baseGroups = c("OpenStreetMap", "CartoDB"),
  #  overlayGroups = c("Orientación", "Municipios", "Ciudades"),
  #  options = layersControlOptions(collapsed = TRUE)
  #)                                                                       

# Mostrar el mapa
#mapa_final

El mapa muestra la orientación del terreno en el Tolima. Predominan las áreas con orientación hacia el norte (color rojo), especialmente en la zona central y occidental del departamento. Esta orientación puede tener implicaciones sobre la radiación solar recibida, temperatura del suelo y condiciones de humedad, factores clave en la planificación agrícola y ambiental (Burrough & McDonnell, 1998). Las demás orientaciones son mucho menos frecuentes, con áreas reducidas hacia el noreste, este y sureste.

Bibliografía 📚

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

Instituto Geográfico Agustín Codazzi – IGAC. (2020). Normas Técnicas para la Adopción del Marco Geodésico de Referencia MAGNA-SIRGAS con Origen Único Nacional. Bogotá, Colombia.

Instituto Geográfico Agustín Codazzi – IGAC. (2022). Boletín Cartográfico No. 27 – Sistema de Coordenadas Origen Único Nacional EPSG:9377. Bogotá, Colombia.

Hijmans, R. J. (2023). terra: Spatial Data Analysis. R package version 1.7-65. https://cran.r-project.org/package=terra

Snyder, J. P. (1987). Map Projections – A Working Manual. U.S. Geological Survey Professional Paper 1395. Washington: United States Government Printing Office.

Franklin, S. E. (2020). Remote Sensing for Biodiversity and Wildlife Management: Synthesis and Applications. McGill-Queen’s Press.

Geopard Tech. (2020). Using terrain attributes to enhance precision agriculture. Geopard Blog.

Geopard Tech. (2022). Topographic factors and their influence on crop yield. Geopard Blog.

Ilich, A., Drăguț, L., & Csillik, O. (2023). MultiscaleDTM: An R package for multiscale terrain analysis using local geomorphometric variables. Computers & Geosciences, 171, 105202. https://doi.org/10.1016/j.cageo.2022.105202

Kokulan, V., Minasny, B., Malone, B. P., & McBratney, A. B. (2018). Mapping soil drainage classes using terrain attributes and Random Forest. Geoderma, 311, 120–129. https://doi.org/10.1016/j.geoderma.2017.09.018

Burrough, P. A., & McDonnell, R. A. (1998). Principles of geographical information systems. Oxford University Press.