Quinto cuaderno en R: Atributos geomorfométricos del terreno en R
Valentina Lasso Rodriguez
26/06/2025
- 1. °| 🧩 Introducción |°
- 2. °| 🛠️ Preparación de la herramienta |°
- 3. °| 🔍 Introducción a MultiscaleDTM |°
- 4. °| 📑 Leer datos de entrada |°
- 5. °| 🔄 Transformando las coordenadas (10%) |°
- 6. °| 📐 Calculando atributos del terreno (20%) |°
- 7. °| 📊 Calculando estadisticas zonales (20%) |°
- 8. °| 🗺️ Graficando la pendiente del terreno (20%) |°
- 9. °| 🧭 Visualizacion alternativa de la pendiente del terreno (15%) |°
- 10. °| 🧮 Visualizando aspectos del terreno (15%) |°
- °💻 ENTORNO INFORMATICO°
1. °| 🧩 Introducción |°
Este Notebook de R ilustra como calcular atributos geomorfométricos de un terreno desde modelos de elevación digital (DEMs). De estas manera para calcular los atributos del terreno usaremos el MultiscaleDTM R package.
2. °| 🛠️ Preparación de la herramienta |°
Antes de iniciar vamos a asegurarnos que los paquetes y librerías que necesitamos están instaladas y cargadas, para darle un buen manejo a esto se recomienda limpiar la memoria antes de iniciar:
# paquetes = c("MultiscaleDTM", "exactextractr")
# install.packages(paquetes)
# Para limpiar la memoria:
rm(list=ls())
# Para cargar las librerías:
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)3. °| 🔍 Introducción a MultiscaleDTM |°
Este paquete de R calcula atributos geomorfológicos del terreno a múltiples escalas a partir de modelos digitales del terreno (DTM) en formato de cuadrícula regular (es decir, mapas raster de elevación o batimetría), utilizando una ventana de análisis especificada.
Funciones principales 🧠
Pendiente, orientación y curvatura
“SlpAsp” calcula la pendiente y orientación a múltiples escalas, según Misiuk et al. (2021), con base en una modificación de los algoritmos tradicionales de pendiente y orientación en ventanas 3 × 3 (Fleming y Hoffer, 1979; Horn et al., 1981; Ritter, 1987). Este algoritmo solo considera un subconjunto de celdas dentro de la ventana de análisis. Específicamente:
Las cuatro celdas ubicadas directamente arriba, abajo, a la izquierda y a la derecha de la celda central (focal) para el método “rook”.
Esas mismas cuatro más las celdas diagonales para el método “queen”.
Todas las celdas ubicadas en el borde de la ventana para el método “boundary”.
4. °| 📑 Leer datos de entrada |°
Primero vamos a subir el raster DEM de nuestro departamento. Se lee usando la librería “terra”
# Cambiar la ruta y el nombre del DEM como se necesite
(dem = terra::rast("C:/Users/stell/Downloads/Valentina/GB2/Notebook 4/elev_valle_z10.tif"))## class : SpatRaster
## size : 3580, 3077, 1 (nrow, ncol, nlyr)
## resolution : 0.0006856326, 0.0006856326 (x, y)
## extent : -77.69531, -75.58562, 2.811443, 5.266008 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_valle_z10.tif
## name : elev_valle_z10
## min value : -2050
## max value : 5373Podemos reducir la resolución del DEM para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Ahora leemos el “gpkg” de municipios de nuestro departamento:
(munic <- sf::st_read("C:/Users/stell/Downloads/Valentina/valle_munic.gpkg"))## Reading layer `valle_munic' from data source
## `C:\Users\stell\Downloads\Valentina\valle_munic.gpkg' using driver `GPKG'
## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
## Geodetic CRS: MAGNA-SIRGAS## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.54977 ymin: 3.091239 xmax: -75.70724 ymax: 5.047394
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 76 001 76001 VALLE DEL CAUCA CALI
## 2 76 020 76020 VALLE DEL CAUCA ALCALÁ
## 3 76 036 76036 VALLE DEL CAUCA ANDALUCÍA
## 4 76 041 76041 VALLE DEL CAUCA ANSERMANUEVO
## 5 76 054 76054 VALLE DEL CAUCA ARGELIA
## 6 76 100 76100 VALLE DEL CAUCA BOLÍVAR
## 7 76 109 76109 VALLE DEL CAUCA BUENAVENTURA
## 8 76 111 76111 VALLE DEL CAUCA GUADALAJARA DE BUGA
## 9 76 113 76113 VALLE DEL CAUCA BUGALAGRANDE
## 10 76 122 76122 VALLE DEL CAUCA CAICEDONIA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1536 MUNICIPIO 563.04092 2023
## 2 Ordenanza 12 de Marzo 31 de 1919 MUNICIPIO 62.35695 2023
## 3 Ordenanza 38 de Abril 25 de 1921 MUNICIPIO 110.46042 2023
## 4 Ordenanza 29 de 1925 MUNICIPIO 305.45167 2023
## 5 Ordenanza 15 de 1956 MUNICIPIO 90.79604 2023
## 6 Ordenanza 2 de 1884 MUNICIPIO 743.57207 2023
## 7 1872 MUNICIPIO 6292.04790 2023
## 8 1555 MUNICIPIO 825.87759 2023
## 9 1791 MUNICIPIO 396.78132 2023
## 10 Ordenanza 21 del 20 de Abril de 1923 MUNICIPIO 166.98369 2023
## shape_Leng shape_Area geom
## 1 1.1644109 0.045732983 MULTIPOLYGON (((-76.59175 3...
## 2 0.4895551 0.005077916 MULTIPOLYGON (((-75.83262 4...
## 3 0.6651004 0.008984851 MULTIPOLYGON (((-76.22406 4...
## 4 0.9459396 0.024871121 MULTIPOLYGON (((-76.01558 4...
## 5 0.4538261 0.007391017 MULTIPOLYGON (((-76.14316 4...
## 6 1.5372119 0.060482648 MULTIPOLYGON (((-76.27584 4...
## 7 6.6181260 0.510748888 MULTIPOLYGON (((-77.2381 4....
## 8 2.0065360 0.067158351 MULTIPOLYGON (((-76.31608 3...
## 9 1.0336063 0.032279294 MULTIPOLYGON (((-76.15131 4...
## 10 0.6465900 0.013590500 MULTIPOLYGON (((-75.8539 4....Ahora unimos el objeto que creamos antes con los limites de nuestro departamento:
# Re-proyectar el DEM al mismo CRS que los municipios
dem2 <- terra::project(dem2, munic)
# Luego se puede recortar sin problemas
dem3 = terra::crop(dem2,munic, mask=TRUE)5. °| 🔄 Transformando las coordenadas (10%) |°
Para computar los atributos del terreno, el DEM debe estar en coordenadas planas, por lo que se modificaran los datos de entrada:
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1430, 1353, 1 (nrow, ncol, nlyr)
## resolution : 152.1184, 152.1184 (x, y)
## extent : 4494105, 4699921, 1900030, 2117560 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_valle_z10
## min value : -252.9045
## max value : 4201.2246De igual manera lo haremos con el vector de los municipios:
(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: 4494194 ymin: 1900273 xmax: 4699750 ymax: 2116536
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 76 001 76001 VALLE DEL CAUCA CALI
## 2 76 020 76020 VALLE DEL CAUCA ALCALÁ
## 3 76 036 76036 VALLE DEL CAUCA ANDALUCÍA
## 4 76 041 76041 VALLE DEL CAUCA ANSERMANUEVO
## 5 76 054 76054 VALLE DEL CAUCA ARGELIA
## 6 76 100 76100 VALLE DEL CAUCA BOLÍVAR
## 7 76 109 76109 VALLE DEL CAUCA BUENAVENTURA
## 8 76 111 76111 VALLE DEL CAUCA GUADALAJARA DE BUGA
## 9 76 113 76113 VALLE DEL CAUCA BUGALAGRANDE
## 10 76 122 76122 VALLE DEL CAUCA CAICEDONIA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1536 MUNICIPIO 563.04092 2023
## 2 Ordenanza 12 de Marzo 31 de 1919 MUNICIPIO 62.35695 2023
## 3 Ordenanza 38 de Abril 25 de 1921 MUNICIPIO 110.46042 2023
## 4 Ordenanza 29 de 1925 MUNICIPIO 305.45167 2023
## 5 Ordenanza 15 de 1956 MUNICIPIO 90.79604 2023
## 6 Ordenanza 2 de 1884 MUNICIPIO 743.57207 2023
## 7 1872 MUNICIPIO 6292.04790 2023
## 8 1555 MUNICIPIO 825.87759 2023
## 9 1791 MUNICIPIO 396.78132 2023
## 10 Ordenanza 21 del 20 de Abril de 1923 MUNICIPIO 166.98369 2023
## shape_Leng shape_Area geom
## 1 1.1644109 0.045732983 MULTIPOLYGON (((4600987 195...
## 2 0.4895551 0.005077916 MULTIPOLYGON (((4685860 208...
## 3 0.6651004 0.008984851 MULTIPOLYGON (((4642160 202...
## 4 0.9459396 0.024871121 MULTIPOLYGON (((4665634 209...
## 5 0.4538261 0.007391017 MULTIPOLYGON (((4651417 208...
## 6 1.5372119 0.060482648 MULTIPOLYGON (((4636542 205...
## 7 6.6181260 0.510748888 MULTIPOLYGON (((4529330 200...
## 8 2.0065360 0.067158351 MULTIPOLYGON (((4631835 199...
## 9 1.0336063 0.032279294 MULTIPOLYGON (((4650286 203...
## 10 0.6465900 0.013590500 MULTIPOLYGON (((4683364 204...6. °| 📐 Calculando atributos del terreno (20%) |°
Vamos a aprender cómo calcular pendiente (slope) y orientación (aspect) del terreno:
## ?SlpAspNota: Este código nos muestra ayuda correspondiente con la función en una ventana a la derecha
Ahora ya sabiendo la información podemos usar la funció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 : 1430, 1353, 2 (nrow, ncol, nlyr)
## resolution : 152.1184, 152.1184 (x, y)
## extent : 4494105, 4699921, 1900030, 2117560 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 1.119667e-04
## max values : 74.02856, 3.599985e+02Donde “w” es el argumento que define la ventana usada para calcular los atributos geomorfologicos, por lo tanto, indica una ventana de 3 filas por 3 columnas, es decir, una ventana cuadrada de 3x3 píxeles centrada en cada celda del raster. Cabe recalcar que se puede usar valores mayores como 5x5, pero se pierde un poco el detalle.
Para el argumento “method” define la conectividad entre pixeles, en este caso se uso la distribución “queen” que haciendo una similitud al ajedrez, por lo que viéndolo aplicado aquí, un pixel central esta conectado a sus 8 pixeles adyacentes, tanto en los verticales, horizontales y diagonales.Esta distribución ayuda a captar mejor las posibles variaciones en el terreno.
Ahora, vamos a dividir el objeto en dos partes:
Primera capa:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1430, 1353, 1 (nrow, ncol, nlyr)
## resolution : 152.1184, 152.1184 (x, y)
## extent : 4494105, 4699921, 1900030, 2117560 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 74.02856terra::hist(slope,
main = "Valle del Cauca's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 52% of the cells was used (of which 54% was NA)
Segunda capa:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1430, 1353, 1 (nrow, ncol, nlyr)
## resolution : 152.1184, 152.1184 (x, y)
## extent : 4494105, 4699921, 1900030, 2117560 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.119667e-04
## max value : 3.599985e+02terra::hist(aspect,
main = "Valle del Cauca's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 52% of the cells was used (of which 54% was NA)
Adicionalmente vamos a convertir la pendiente de grados a porcentaje
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1430, 1353, 1 (nrow, ncol, nlyr)
## resolution : 152.1184, 152.1184 (x, y)
## extent : 4494105, 4699921, 1900030, 2117560 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 349.3986terra::hist(slope_perc,
main = "Valle del Cauca's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 52% of the cells was used (of which 54% was NA)
7. °| 📊 Calculando estadisticas zonales (20%) |°
Una operación de estadísticas zonales calcula las estadísticas sobre los valores de las celdas de un raster, estas zonas definidas por un conjunto de datos. En este caso, estamos interesados en calcular los valores promedio de la pendiente por municipio.
Primero, podemos conocer la clasificación de pendientes según el IGAC, se usa para fines agrologicos:
| Dirección Cardinal | Rango (°) | Color típico |
|---|---|---|
| Plano (sin pendiente) | -1 | Gris |
| Norte | 0–22.5 y 337.5–360 | Rojo |
| Noreste | 22.5–67.5 | Naranja |
| Este | 67.5–112.5 | Amarillo |
| Sureste | 112.5–157.5 | Verde claro |
| Sur | 157.5–202.5 | Verde oscuro |
| Suroeste | 202.5–247.5 | Azul claro |
| Oeste | 247.5–292.5 | Azul |
| Noroeste | 292.5–337.5 | Púrpura |
Ahora si, primero podemos empezar:
##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)Podemos calcular 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% | |=== | 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] 19.1932716 7.0536819 6.5234451 24.4195595 29.2239342 33.5480080
## [7] 14.4227972 29.3312092 10.4425707 19.4222984 29.9164867 0.5738842
## [13] 8.4911633 30.0080242 39.9681816 34.8505516 21.2526875 32.7035942
## [19] 30.7923889 27.2616692 10.1506596 19.0136471 19.5517540 14.8743944
## [25] 9.2575541 10.9017820 19.2761707 26.3611240 15.6934271 23.0123882
## [31] 16.8672638 15.2212105 24.7992954 20.7394142 26.5074520 26.5728569
## [37] 7.9285340 31.6442127 18.7146034 17.6516132 19.1481228 5.5441055hist(munic$mean_slope,
main = "Valle del Cauca's mean slope",
xlab = "Slope (in percentage)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 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] 1 2 1 5 5 5 1 5 1 5 5 1 1 5 5 5 1 5 5 5 1 1 4 1 1 1 1 1 4 5 5 1 5 5 5 5 2 5
## [39] 4 4 5 1hist(munic$class,
main = "Valle del Cauca's reclassified slope",
xlab = "Slope (as a category)")
Transformemos la pendiente otra vez a coordenadas geográficas
## La pendiente reclasificada
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1439, 1357, 1 (nrow, ncol, nlyr)
## resolution : 0.001371234, 0.001371234 (x, y)
## extent : -77.56098, -75.70021, 3.085094, 5.058299 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.000
## max value : 292.879## La pendiente en porcentaje
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1439, 1357, 1 (nrow, ncol, nlyr)
## resolution : 0.001371234, 0.001371234 (x, y)
## extent : -77.56098, -75.70021, 3.085094, 5.058299 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 297.06778. °| 🗺️ Graficando la pendiente del terreno (20%) |°
Primero vamos a necesitar una paleta de colores para la pendiente:
### Primero expandir el numero de colores para mejorar la visualización
## Encontrar una guía en esta pagina:
# https://r-graph-gallery.com/42-colors-names.html
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ahora vamos a graficar:
leaflet(munic) %>%
addTiles() %>%
setView(lng = -76.6, lat = 4.1, zoom = 8) %>%
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 Valle del Cauca (%)")## 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'Nota: Se puede dar click y moverlo por distintos municipios para ver el nombre y el valor de pendiente en porcentaje más frecuente en cada uno.
9. °| 🧭 Visualizacion alternativa de la pendiente del terreno (15%) |°
Para ver la pendiente del terreno en grados podemos hacerlo con el siguiente código:
# Primero vas o ver que el crs esté en WGS84
slope <- project(slope, "EPSG:4326")
munic <- st_transform(munic, 4326)
# Paleta de colores para pendiente en grados
pal_degrees <- colorNumeric(palette = "YlOrRd", domain = values(slope), na.color = "transparent")
# Para organizar los centroides
munic_centroids <- sf::st_centroid(munic)## Warning: st_centroid assumes attributes are constant over geometriesmunic_centroids$lon <- sf::st_coordinates(munic_centroids)[, 1]
munic_centroids$lat <- sf::st_coordinates(munic_centroids)[, 2]
munic_centroids$label <- paste0(
munic_centroids$mpio_cnmbr, ": ",
round(munic_centroids$mean_slope, 1), "°"
)
# El mapa
leaflet(munic) %>%
addTiles() %>%
setView(lng = -76.6, lat = 4.1, zoom = 8) %>%
addRasterImage(slope, colors = pal_degrees, opacity = 0.8) %>%
addPolygons(
color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = ~paste0(
"<b>Municipio:</b> ", mpio_cnmbr, "<br>",
"<b>Clase de pendiente:</b> ", class
)
) %>%
addLegend(pal = palredgreen, values = values(slope),
title = "Terrain slope in Valle del Cauca (°)")10. °| 🧮 Visualizando aspectos del terreno (15%) |°
Vamos a ver las características del terreno según las especificaciones y usando la paleta de colores convencional basada en la simbologia sugerida por ESRI / ArcGis, de esta manera:
| Dirección Cardinal | Rango (°) | Color típico |
|---|---|---|
| Plano (sin pendiente) | -1 | Gris |
| Norte | 0–22.5 y 337.5–360 | Rojo |
| Noreste | 22.5–67.5 | Naranja |
| Este | 67.5–112.5 | Amarillo |
| Sureste | 112.5–157.5 | Verde claro |
| Sur | 157.5–202.5 | Verde oscuro |
| Suroeste | 202.5–247.5 | Azul claro |
| Oeste | 247.5–292.5 | Azul |
| Noroeste | 292.5–337.5 | Púrpura |
Ahora:
# Clasificación manual según los rangos
m_aspect <- matrix(c(
-1, 0, 0, # plano
0, 22.5, 1, # norte
22.5, 67.5, 2, # noreste
67.5, 112.5, 3, # este
112.5, 157.5, 4, # sureste
157.5, 202.5, 5, # sur
202.5, 247.5, 6, # suroeste
247.5, 292.5, 7, # oeste
292.5, 337.5, 8, # noroeste
337.5, 360, 1 # norte (de nuevo)
), ncol = 3, byrow = TRUE)
aspect_class <- classify(aspect, m_aspect, right = FALSE)
# Paleta de colores sugerida por ArcGIS
aspect_colors <- c(
"gray80", # 0: plano
"red", # 1: norte
"orange", # 2: noreste
"yellow", # 3: este
"lightgreen", # 4: sureste
"darkgreen", # 5: sur
"lightblue", # 6: suroeste
"blue", # 7: oeste
"purple" # 8: noroeste
)
# Reproyectar si es necesario
aspect_class_geo <- project(aspect_class, "EPSG:4326")
# El mapa
pal_aspect <- colorFactor(palette = aspect_colors, domain = 0:8)
leaflet() %>%
addTiles() %>%
setView(lng = -76.6, lat = 4.1, zoom = 8) %>%
addRasterImage(aspect_class, colors = pal_aspect, opacity = 0.8) %>%
addLegend(
pal = pal_aspect,
values = 0:8,
title = "Aspecto del terreno",
position = "topright",
labFormat = function(type, cuts, p) {
c("Plano", "Norte", "NE", "Este", "SE", "Sur", "SO", "Oeste", "NO")
}
)## Warning in colors(.): Some values were outside the color scale and will be
## treated as NAEste mapa muestra hacia qué dirección cardinal se orienta la pendiente del terreno.
Por ejemplo:
-Un pixel rojo indica que la ladera mira hacia el norte. -Un pixel azul indica que la ladera está orientada al oeste. -Zonas grises son áreas planas sin orientación dominante.
Esto es útil en estudios como:
-Erosión, porque las laderas sur en zonas tropicales suelen ser más secas. -Uso del suelo (cultivos, conservación, planeación urbana). -Modelado de microclimas.
°💻 ENTORNO INFORMATICO°
Bibliografia: Lizarazo, I., 2025. Geomorphometric terrain attributes in R. Available at https://rpubs.com/ials2un/geomorphometric
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8 LC_CTYPE=Spanish_Colombia.utf8
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C
## [5] LC_TIME=Spanish_Colombia.utf8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] exactextractr_0.10.0 MultiscaleDTM_1.0 terra_1.8-54
## [4] leaflet_2.2.2 sf_1.0-21 elevatr_0.99.0
##
## loaded via a namespace (and not attached):
## [1] s2_1.1.9 sass_0.4.10 generics_0.1.4 class_7.3-23
## [5] KernSmooth_2.23-26 lattice_0.22-6 digest_0.6.37 magrittr_2.0.3
## [9] rgl_1.3.18 RColorBrewer_1.1-3 evaluate_1.0.3 grid_4.5.0
## [13] fastmap_1.2.0 jsonlite_2.0.0 e1071_1.7-16 DBI_1.2.3
## [17] promises_1.3.3 purrr_1.0.4 scales_1.4.0 crosstalk_1.2.1
## [21] codetools_0.2-20 jquerylib_0.1.4 cli_3.6.5 shiny_1.10.0
## [25] rlang_1.1.6 units_0.8-7 base64enc_0.1-3 cachem_1.1.0
## [29] yaml_2.3.10 raster_3.6-32 tools_4.5.0 dplyr_1.1.4
## [33] httpuv_1.6.16 png_0.1-8 vctrs_0.6.5 R6_2.6.1
## [37] mime_0.13 proxy_0.4-27 lifecycle_1.0.4 classInt_0.4-11
## [41] htmlwidgets_1.6.4 pkgconfig_2.0.3 progressr_0.15.1 bslib_0.9.0
## [45] pillar_1.10.2 later_1.4.2 glue_1.8.0 Rcpp_1.0.14
## [49] tidyselect_1.2.1 xfun_0.52 tibble_3.2.1 rstudioapi_0.17.1
## [53] knitr_1.50 farver_2.1.2 xtable_1.8-4 htmltools_0.5.8.1
## [57] rmarkdown_2.29 wk_0.9.4 compiler_4.5.0 sp_2.2-0