1 Introduccion:

El análisis de atributos geomorfométricos (como la pendiente y la orientación de laderas) es fundamental para la planificación territorial, la conservación de suelos y la evaluación de riesgos de remoción en masa. Para este estudio destacamos el uso de las siguientes librerías clave:

MultiscaleDTM: Herramienta especializada para el cálculo avanzado de geomorfometría a partir de modelos digitales de elevación. * terra y sf: Motor principal para el procesamiento de datos ráster y vectoriales. * exactextractr: Permite la extracción estadística rápida y precisa de valores ráster dentro de polígonos municipales. * leaflet: Plataforma para la creación de cartografía interactiva.

rm(list=ls())
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)
library(sp)

2 Carga y Preprocesamiento del Modelo Digital de Elevación (DEM):

Cargamos el DEM del departamento del Cauca y aplicamos un filtro de seguridad para evitar datos erróneos o valores negativos sobre la superficie terrestre.

ruta = terra::rast("C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno5\\data\\elev_cauca_z10.tif")
## Se establece que el DEM el minimo valor que pueda tomar sea 0
dem <- terra::clamp(ruta, lower = 0, values = TRUE)
## |---------|---------|---------|---------|=========================================                                          
## Vizualozaremos nuestra variable
dem
## class       : SpatRaster
## size        : 4095, 3586, 1  (nrow, ncol, nlyr)
## resolution  : 0.0006863432, 0.0006863432  (x, y)
## extent      : -78.04687, -75.58565, 0.7028454, 3.513421  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## varname     : elev_cauca_z10
## name        : elev_cauca_z10
## min value   :              0
## max value   :           5373

3 Reducción de Resolución y Delimitación Territorial:

Para optimizar el uso de la memoria RAM del equipo y evitar saturaciones durante los cálculos matriciales, agregamos píxeles y recortamos el mapa ajustándolo exactamente a la frontera del departamento.

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

Leeremos los municipios del departamento del cauca para tener una clara delimitacion de las dimensiones del dapartamento

munic <- sf::st_read("C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno5\\data\\cauca_munic.gpkg")
## Reading layer `cauca_munic' from data source 
##   `C:\Users\pc\OneDrive\Documentos\GB2\RSTUDIO\cuaderno5\data\cauca_munic.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
## Geodetic CRS:  WGS 84

Tanto el DEM como los municipos estan en coordenadas geograficas WGS 84, mediante la funcion crop presente en la libreria terra recortaremos dem2 a nuestro objeto geoespacial que corresponde a la region de intres.

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

4 Reproyección del Ráster:

Las fórmulas de pendiente requieren que tanto las distancias horizontales (X, Y) como la altitud (Z) estén expresadas en las mismas unidades (metros). Por ello, transformamos las coordenadas geográficas (grados) al sistema plano oficial de Colombia.

(dem_plane = project(dem3, "EPSG:9377"))
## class       : SpatRaster
## size        : 1727, 1598, 1  (nrow, ncol, nlyr)
## resolution  : 152.4318, 152.4318  (x, y)
## extent      : 4451159, 4694745, 1664036, 1927286  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## name        : elev_cauca_z10
## min value   :              0
## max value   :    5359.402832

5 Reproyectando los Objetos Vectoriales

(munic_plane = sf::st_transform(munic, "EPSG:9377"))
## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4451719 ymin: 1664078 xmax: 4694583 ymax: 1926529
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
##    dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr   mpio_cnmbr
## 1          19        001      19001      CAUCA      POPAYÁN
## 2          19        022      19022      CAUCA     ALMAGUER
## 3          19        050      19050      CAUCA      ARGELIA
## 4          19        075      19075      CAUCA       BALBOA
## 5          19        100      19100      CAUCA      BOLÍVAR
## 6          19        110      19110      CAUCA BUENOS AIRES
## 7          19        130      19130      CAUCA      CAJIBÍO
## 8          19        137      19137      CAUCA      CALDONO
## 9          19        142      19142      CAUCA       CALOTO
## 10         19        212      19212      CAUCA      CORINTO
##                            mpio_crslc mpio_tipo mpio_narea mpio_nano shape_Leng
## 1                                1537 MUNICIPIO   480.3814      2024   1.374118
## 2                                1799 MUNICIPIO   236.8516      2024   0.710253
## 3  Ordenanza 2 de Noviembre 8 de 1967 MUNICIPIO   776.4555      2024   1.224153
## 4   Ordenanza 1 de Octubre 20 de 1967 MUNICIPIO   413.5523      2024   1.016327
## 5                                1793 MUNICIPIO   798.3494      2024   1.802787
## 6                                1851 MUNICIPIO   435.1190      2024   1.632520
## 7                                1824 MUNICIPIO   552.5514      2024   1.381852
## 8                                1746 MUNICIPIO   354.8557      2024   1.102961
## 9                                1543 MUNICIPIO   264.8714      2024   1.315453
## 10                               1868 MUNICIPIO   325.9754      2024   0.828674
##    shape_Area                           geom
## 1  0.03897017 MULTIPOLYGON (((4583956 184...
## 2  0.01919657 MULTIPOLYGON (((4576303 177...
## 3  0.06288063 MULTIPOLYGON (((4520521 183...
## 4  0.03349049 MULTIPOLYGON (((4539082 179...
## 5  0.06468337 MULTIPOLYGON (((4562832 179...
## 6  0.03530827 MULTIPOLYGON (((4581114 191...
## 7  0.04481975 MULTIPOLYGON (((4573051 186...
## 8  0.02880277 MULTIPOLYGON (((4626166 187...
## 9  0.02150781 MULTIPOLYGON (((4619907 191...
## 10 0.02648009 MULTIPOLYGON (((4641747 189...

6 Calculo de Atributos Geomorfométricos:

Haciendo uso de la librería MultiscaleDTM, calculamos la pendiente y el aspecto utilizando una ventana móvil de 3x3 celdas, aplicando el método QUEEN y corrigiendo celdas sin inclinación.

(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        : 1727, 1598, 2  (nrow, ncol, nlyr)
## resolution  : 152.4318, 152.4318  (x, y)
## extent      : 4451159, 4694745, 1664036, 1927286  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## names       :    slope,     aspect
## min values  :        0,   0.000018
## max values  : 68.76312, 359.999962

7 Extracción de la Capa de pendiente (slope):

El objeto slp_asp generado por la librería MultiscaleDTM es un ráster multicapa que contiene tanto la pendiente como el aspecto. La función subset() nos permite desglosar este objeto y extraer capas individuales. Al utilizar el índice 1, estamos aislando específicamente la matriz de datos de la pendiente (la inclinación del terreno) para realizar su análisis estadístico y espacial de forma independiente.

(slope = subset(slp_asp, 1))
## class       : SpatRaster
## size        : 1727, 1598, 1  (nrow, ncol, nlyr)
## resolution  : 152.4318, 152.4318  (x, y)
## extent      : 4451159, 4694745, 1664036, 1927286  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## name        :    slope
## min value   :        0
## max value   : 68.76312

A continuacion representaremos mediante histogramas los datos y analizaremos su distribución:

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

8 Nota: Extracción de la Capa de Orientación (Aspecto):

El objeto slp_asp generado por la librería MultiscaleDTM es un ráster multicapa que contiene tanto la pendiente como el aspecto. La función subset() nos permite desglosar este objeto y extraer capas individuales. Al utilizar el índice 2, estamos aislando específicamente la matriz de datos del aspecto (la orientación de las laderas) para realizar su análisis estadístico y espacial de forma independiente.

(aspect = subset(slp_asp, 2))
## class       : SpatRaster
## size        : 1727, 1598, 1  (nrow, ncol, nlyr)
## resolution  : 152.4318, 152.4318  (x, y)
## extent      : 4451159, 4694745, 1664036, 1927286  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## name        :     aspect
## min value   :   0.000018
## max value   : 359.999962
terra::hist(aspect, 
     main = "Cauca's aspect", 
     xlab = "Aspect (in degrees)")

Expreseremos la pendiente en porcentaje mediante la siguiente conversion:

(slope_perc = tan(slope*(pi/180))*100)
## class       : SpatRaster
## size        : 1727, 1598, 1  (nrow, ncol, nlyr)
## resolution  : 152.4318, 152.4318  (x, y)
## extent      : 4451159, 4694745, 1664036, 1927286  (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s)   : memory
## name        :      slope
## min value   :          0
## max value   : 257.323996
terra::hist(slope_perc, 
     main = "Cauca's slope", 
     xlab = "Slope (in percentage)")

9 Reclasificación Categorizada de la Pendiente:

Dividimos la pendiente continua en 7 clases agronómicas estándar para clasificar la aptitud del terreno.

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

10 Análisis Zonal por Municipios:

Extraemos métricas para caracterizar el relieve promedio de cada municipio del departamento del Cauca.

(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))
##   |                                                                              |                                                                      |   0%  |                                                                              |==                                                                    |   2%  |                                                                              |===                                                                   |   5%  |                                                                              |=====                                                                 |   7%  |                                                                              |=======                                                               |  10%  |                                                                              |========                                                              |  12%  |                                                                              |==========                                                            |  14%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  21%  |                                                                              |=================                                                     |  24%  |                                                                              |==================                                                    |  26%  |                                                                              |====================                                                  |  29%  |                                                                              |======================                                                |  31%  |                                                                              |=======================                                               |  33%  |                                                                              |=========================                                             |  36%  |                                                                              |===========================                                           |  38%  |                                                                              |============================                                          |  40%  |                                                                              |==============================                                        |  43%  |                                                                              |================================                                      |  45%  |                                                                              |=================================                                     |  48%  |                                                                              |===================================                                   |  50%  |                                                                              |=====================================                                 |  52%  |                                                                              |======================================                                |  55%  |                                                                              |========================================                              |  57%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  62%  |                                                                              |=============================================                         |  64%  |                                                                              |===============================================                       |  67%  |                                                                              |================================================                      |  69%  |                                                                              |==================================================                    |  71%  |                                                                              |====================================================                  |  74%  |                                                                              |=====================================================                 |  76%  |                                                                              |=======================================================               |  79%  |                                                                              |=========================================================             |  81%  |                                                                              |==========================================================            |  83%  |                                                                              |============================================================          |  86%  |                                                                              |==============================================================        |  88%  |                                                                              |===============================================================       |  90%  |                                                                              |=================================================================     |  93%  |                                                                              |===================================================================   |  95%  |                                                                              |====================================================================  |  98%  |                                                                              |======================================================================| 100%
##  [1] 12.7062006 34.4478607 42.5884094 34.2758980 27.2228680 21.2370110
##  [7] 15.2377043 21.4560223 18.5817146 29.1378345 31.0206146 24.0970097
## [13]  1.0455663 12.4084435 26.5436401 33.1763306 23.1620789 32.6481247
## [19] 20.0335407 18.0960045 23.0177364 26.4700069  0.6239001 35.3344116
## [25] 18.5916405 15.1034260 13.3238077  0.5157506 19.4977303 24.2509060
## [31] 24.4898586 11.2066650 30.5933552 23.7175999 23.4066830 28.3654919
## [37] 28.6634159 12.7001238 13.5686045 36.8871307 17.3526478  0.6462926
hist(munic$mean_slope, 
     main = "Cauca's mean slope", 
     xlab = "Slope (in percentage)")

extraemos métricas categóricas utilizando la función de ‘moda’ (mode). El algoritmo analiza el ráster de la pendiente previamente reclasificada y determina cuál es la categoría dominante dentro de los límites de cada municipio. Esto nos permite clasificar al municipio según su paisaje más representativo, almacenando este valor en la nueva columna class.

(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))
##   |                                                                              |                                                                      |   0%  |                                                                              |==                                                                    |   2%  |                                                                              |===                                                                   |   5%  |                                                                              |=====                                                                 |   7%  |                                                                              |=======                                                               |  10%  |                                                                              |========                                                              |  12%  |                                                                              |==========                                                            |  14%  |                                                                              |============                                                          |  17%  |                                                                              |=============                                                         |  19%  |                                                                              |===============                                                       |  21%  |                                                                              |=================                                                     |  24%  |                                                                              |==================                                                    |  26%  |                                                                              |====================                                                  |  29%  |                                                                              |======================                                                |  31%  |                                                                              |=======================                                               |  33%  |                                                                              |=========================                                             |  36%  |                                                                              |===========================                                           |  38%  |                                                                              |============================                                          |  40%  |                                                                              |==============================                                        |  43%  |                                                                              |================================                                      |  45%  |                                                                              |=================================                                     |  48%  |                                                                              |===================================                                   |  50%  |                                                                              |=====================================                                 |  52%  |                                                                              |======================================                                |  55%  |                                                                              |========================================                              |  57%  |                                                                              |==========================================                            |  60%  |                                                                              |===========================================                           |  62%  |                                                                              |=============================================                         |  64%  |                                                                              |===============================================                       |  67%  |                                                                              |================================================                      |  69%  |                                                                              |==================================================                    |  71%  |                                                                              |====================================================                  |  74%  |                                                                              |=====================================================                 |  76%  |                                                                              |=======================================================               |  79%  |                                                                              |=========================================================             |  81%  |                                                                              |==========================================================            |  83%  |                                                                              |============================================================          |  86%  |                                                                              |==============================================================        |  88%  |                                                                              |===============================================================       |  90%  |                                                                              |=================================================================     |  93%  |                                                                              |===================================================================   |  95%  |                                                                              |====================================================================  |  98%  |                                                                              |======================================================================| 100%
##  [1] 4 5 5 5 5 5 2 5 5 5 5 5 1 2 5 5 4 5 5 4 1 5 1 5 5 1 2 1 4 4 5 1 5 4 4 5 5 2
## [39] 2 5 4 1
hist(munic$class, 
     main = "Cauca's reclassified slope", 
     xlab = "Slope (as a category)")

# Nota: Para visualizar espacialmente los resultados en un mapa dinámico con mapas base reales, utilizaremos la plataforma leaflet. Sin embargo, esta librería requiere estrictamente que las capas espaciales estén en coordenadas geográficas (Latitud/Longitud). Por consiguiente, procedemos a reproyectar temporalmente nuestros rásteres de pendiente y aspecto desde el sistema MAGNA-SIRGAS hacia el estándar WGS 84 (EPSG:4326).

(rc.geo = project(rc, "EPSG:4326"))
## class       : SpatRaster
## size        : 1736, 1597, 1  (nrow, ncol, nlyr)
## resolution  : 0.00137266, 0.00137266  (x, y)
## extent      : -77.93644, -75.7443, 0.955096, 3.338034  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## name        :     slope
## min value   :         0
## max value   : 216.17717
(slope.geo = project(slope_perc, "EPSG:4326"))
## class       : SpatRaster
## size        : 1736, 1597, 1  (nrow, ncol, nlyr)
## resolution  : 0.00137266, 0.00137266  (x, y)
## extent      : -77.93644, -75.7443, 0.955096, 3.338034  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s)   : memory
## name        :      slope
## min value   :          0
## max value   : 235.120575

11 Mapa Interactivo de la pendiente del Terreno:

palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
  na.color = "transparent")
leaflet(munic) %>% addTiles() %>% setView(-76.8259,2.7049, 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 = "Terrain slope in Cauca (%)")

12 Mapa Interactivo del Aspecto del Terreno

# 1. Reproyectar el ráster de aspecto a coordenadas geográficas (necesario para leaflet)
aspect.geo <- terra::project(aspect, "EPSG:4326")

# 2. Crear una paleta de colores 
# Nota: Empezamos y terminamos en "red" para que el Norte (0 y 360) se vea igual
pal_aspect <- colorNumeric(
  palette = c("red", "orange", "yellow", "green", "cyan", "blue", "purple", "red"), 
  domain = c(0, 360),
  na.color = "transparent"
)

# 3. Generar el mapa interactivo con Leaflet
leaflet(munic) %>% 
  addTiles() %>% 
  setView(-76.8259, 2.7049, 9) %>% 
  # Agregar la capa del ráster de Aspecto
  addRasterImage(
    aspect.geo, 
    colors = pal_aspect, 
    opacity = 0.7, 
    group = "Aspecto"
  ) %>%
  # Agregar los límites de los municipios encima para referencia
  addPolygons(
    color = "black", 
    weight = 1.0, 
    smoothFactor = 0.5,
    opacity = 0.6, 
    fill = FALSE
  ) %>%
  # Agregar la leyenda explicativa
  addLegend(
    pal = pal_aspect, 
    values = c(0, 360), # Forzamos el rango completo de la brújula
    title = "Terrain aspect in Cauca",
    position = "bottomright",
    labFormat = labelFormat(suffix = "°")
  )