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.
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)## |---------|---------|---------|---------|=========================================
## 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
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.
## |---------|---------|---------|---------|=========================================
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.
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.
## 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
## 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...
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
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.
## 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:
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.
## 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
Expreseremos la pendiente en porcentaje mediante la siguiente conversion:
## 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
Dividimos la pendiente continua en 7 clases agronómicas estándar para clasificar la aptitud del terreno.
Extraemos métricas para caracterizar el relieve promedio de cada municipio del departamento del Cauca.
## | | | 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
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.
## | | | 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
# 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).
## 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
## 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
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 (%)")