Variables Geomorfometricas del Departamento del Huila
Cristian Felipe Hincapie Romero
15/06/2025
1. Introducción
En este cuaderno aprenderemos a calcular atributos geomorfometricos del terreno a partir de modelos digitales de elevación cuadriculados. El objetivo es aprender conceptos y herramientas básicas de geoinformartica para si mismo y para quien le interese.
Usaremos la libreria MultiscaleDTM para calcular atributos del terreno.
2. Parte Práctica
En primer lugar, cargaremos las librerias que vamos a utilizar (¡importante!, recuerda limpiar la consola antes de empezar con el codigo rm(list=ls())), si no tenemos instalada alguna de las librerias podemos instalarla en la consola con install.packages(“nombre_de_la_libreria”)
podemos eliminar los mensajes innecesarios o los mensajes de alerta con message = FALSE o warning = FALSE respectivamente.
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)Si nos aparece algún error, es importante revisar lo que nos indica dicho error, muchas veces puede ser algún conflicto por funciones de librerias distintas que tienen el mismo nombre
3. MultiscaleDTM
Este paquete de R en especifíco permite calcular atributos geomorfométricos multiescala a partir de Modelos Digitales del Terreno (DTM) o de Batimetría (MDB). Utiliza una ventana móvil definida por el usuario para procesar rásters de elevación con celdas regulares (cuadriculados), generando métricas del terreno adaptables a diferentes escalas espaciales.
SlpAsp calcula la pendiente y el aspecto multiescala según Misiuk et al. (2021), que es una modificación de los algoritmos tradicionales de pendiente y aspecto 3 x 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 en los bordes (arriba, abajo, izquierda y derecha) de la celda focal para el método “rook”, cuatro esquinas adicionales para el método “queen”, o todas las celdas del borde para el método “boundary”.
4. Lectura de los Datos
Vamos a cargar el raster DEM que descargamos para nuestro departamento, para leer el DEM vamos a hacer uso de la libreria terra.
(dem = terra::rast("C:/Users/TEMP/Documents/GB2/P1/DATOS/ELEVACION/HUILA_ELEVATION.tif"))## class : SpatRaster
## size : 3582, 3586, 1 (nrow, ncol, nlyr)
## resolution : 0.0006862561, 0.0006862561 (x, y)
## extent : -76.64062, -74.17971, 1.406085, 3.864255 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : HUILA_ELEVATION.tif
## name : HUILA_ELEVATION
## min value : -747
## max value : 5372Reducimos la resolución para que no hallan problemas con nuestro computador al costarle trabajar con los datos.
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Vamos a leer nuestros municipios apoyandonos de la liberia sf
(munic <- sf::st_read("C:/Users/TEMP/Documents/GB2/P1/DATOS/DEPARTAMENTOS/HUILA.gpkg"))## Reading layer `HUILA' from data source
## `C:\Users\TEMP\Documents\GB2\P1\DATOS\DEPARTAMENTOS\HUILA.gpkg'
## using driver `GPKG'
## Simple feature collection with 37 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS: MAGNA-SIRGASEl conjunto de datos esta en coordenadas geográficas.
Cortamos el objeto dem2 acorde a los límites de nuestro de departamento.
dem3 = terra::crop(dem2, munic, mask= TRUE)## Warning: [crop] CRS do not match5. Transformando Coordenadas
Para calcular correctamente los atributos del terreno, el modelo digital de elevación (DEM) debe encontrarse en un sistema de coordenadas planas. Como se recordará, desde 2020 el IGAC adoptó un sistema con un único origen nacional. ¿Cuál es el código EPSG correspondiente a este sistema?
Ajustamos los dos fragmentos de código que se muestran a continuación para transformar los datos de entrada al sistema de coordenadas planas nacionales vigente en Colombia.
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : HUILA_ELEVATION
## min value : 315.8539
## max value : 5280.1484También es necesario aplicar una transformación equivalente al objeto vectorial para que coincida con el sistema de coordenadas planas nacionales.
(munic_plane = sf::st_transform(munic, "EPSG:9377"))6. Cálculo de Atributos del Terreno
Vamos a calcular la pendiente y el aspecto del terreno. Para ello, dirígete a la consola y escribe ?SlpAsp. Esto abrirá en la ventana de ayuda (a la derecha) la documentación relacionada con esta función.
Este método calcula la pendiente y orientación del terreno utilizando un enfoque multiescala, apoyado en el algoritmo desarrollado por Misiuk y colaboradores en 2021. A diferencia de los métodos tradicionales que se limitan a una ventana de análisis de 3x3 celdas, este permite trabajar con diferentes tamaños de ventana. Además, el algoritmo fue adaptado para aceptar ventanas de forma rectangular, no solamente cuadradas.
Ahora llamemos a la función SlpAsp.
(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 : 1669, 1619, 2 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 1.145656e-04
## max values : 78.29077, 3.599998e+02Ahora, dividamos el objeto slp_asp en dos partes.
Capa 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 78.29077¿Cuál es la distribución de los valores de la pendiente?.
terra::hist(slope,
main = "Pendiente en el Huila",
xlab = "Pendiente (en grados)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
Capa 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.145656e-04
## max value : 3.599998e+02terra::hist(aspect,
main = "Huila Orientacion",
xlab = "Orientacion (en grados)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
De paso, convirtamos los grados de pendiente en porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 482.4902terra::hist(slope_perc,
main = "Pendiente en el Huila",
xlab = "Pendiente (En porcentaje)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
7. Cálculo de Estadísticas Zonales
Este método calcula la pendiente y orientación del terreno utilizando un enfoque multiescala, apoyado en el algoritmo desarrollado por Misiuk y colaboradores en 2021. A diferencia de los métodos tradicionales que se limitan a una ventana de análisis de 3x3 celdas, este permite trabajar con diferentes tamaños de ventana. Además, el algoritmo fue adaptado para aceptar ventanas de forma rectangular, no solamente cuadradas.
Primero, conozcamos la clasificación de pendientes del IGAC para fines agronómicos:
¡Manos a la obra!.
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)Calculemos las estadísticas zonales utilizando exactextractr:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 23.069536 22.971500 13.887896 18.557335 32.051628 14.012997 23.509218
## [8] 18.159473 27.879963 26.955832 23.962423 23.030651 26.672739 18.523106
## [15] 31.574409 15.578115 22.617243 23.450455 27.235809 21.257706 20.056879
## [22] 18.243835 19.343857 20.165983 19.859833 22.496557 21.036348 25.655258
## [29] 32.813728 24.513773 21.133305 17.930473 19.953043 30.260111 21.859444
## [36] 7.383247 12.466011hist(munic$mean_slope,
main = "Pendiente media del Huila",
xlab = "Pendiente (en porcentaje)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4hist(munic$class,
main = "Pendiente reclasificada del Huila",
xlab = "Pendiente (como categoria)")
Ahora, volvamos a transformar la pendiente en coordenadas geográficas:
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 452.9977(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 452.99778. Ploteando la Pend. del Terreno
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ahora podemos plotear, puedes manipular en el código para ajustalo a tus gustos.
leaflet(munic) %>% addTiles() %>% setView(-75.5, 2.6, 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 = "Pendiente del terreno en Huila (%)")9. Visualización Alter. de la Pend. del Terreno
Vamos a visualizar las pendientes del departamento de Huila, pero a diferencia del ejemplo anterior, lo haremos en grados y con las etiquetas que señalan los municipios de nuestro departamento. Utilizamos la función jitter para que los nombres de los municipios no queden sobrepuestos, y para poder plotear hacemos un proceso similar al que hicimos con su pendiente en porcentaje.
Convertimos nuestro objeto en un raster, anterior a la aplicación de esta función era un objeto vectorial.
(munic$mean_slope_deg <- exactextractr::exact_extract(slope, munic, 'mean'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 12.644789 12.797388 7.769434 10.266978 17.480240 7.859426 12.981675
## [8] 9.856176 15.294997 14.747098 13.205628 12.670143 14.697580 10.204793
## [15] 17.106588 8.650787 12.526482 12.855496 15.001105 11.832603 11.193542
## [22] 10.081303 10.848270 11.223944 11.026142 11.898920 11.653651 14.111075
## [29] 17.829418 13.594316 11.653077 9.940928 10.950447 16.399023 12.131503
## [36] 4.158818 7.014504No resulta necesario, pero podemos volvere a visualizar el histograma de la pendiente.
hist(munic$mean_slope_deg,
main = "Pendiente media del Huila",
xlab = "Pendiente (en grados)")
Nuevamente, transformamos un objeto vectorial en un raster.
(munic$class_deg <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4Reprojectamos nuestro objeto Spat.Raster para que tenga las coordenadas con las que estamos trabajando y no halla errores por un disparidad de las mismas.
(rc.geo_deg = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 452.9977De nuevo, repetimos el proceso con el objeto SpatRaster, slope.geo_dg.
(slope.geo_dg = project(slope, "EPSG:4326"))## class : SpatRaster
## size : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 77.37261Para poder colocar las etiquetas de los municipios y que no se sobrepongan unas con otras a causa de situarlas en las mismas coordenadas, tenemos que colocar centroides en cada una de las coordenadas individuales de cada municipio.
munic_centroids <- st_centroid(munic)coords <- st_coordinates(munic_centroids)Usamos set.seed(123) para que cada vez que ejecutemos el codigo obtendremos los mismos valores en cada coordenada, jitter hace referencia a un “ruido aleatorio” que hace que hayan unos minimos valores al azar al rededor de las coordenadas para que las etiquetas no se sobrepongan.
set.seed(123)
coords_jitter <- data.frame(
lng = jitter(coords[, 1], amount = 0.05),
lat = jitter(coords[, 2], amount = 0.05),
label = munic$mpio_cnmbr
)Usamos la paleta de colores anteriormente utilizada para este caso.
palredgreen_dg <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Y ploteamos para poder observar las pendientes del departamento del Huila en grados y las etiquetas de cada municipio para podernos situar mejor.
leaflet(munic) %>% addTiles() %>% setView(-75.5, 2.6, 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_dg, "<br>")) %>%
addRasterImage(slope.geo_dg, colors = palredgreen_dg, opacity = 0.8) %>%
addLabelOnlyMarkers(
data = coords_jitter,
lng = ~lng, lat = ~lat,
label = ~label,
labelOptions = labelOptions(
noHide = TRUE,
direction = "auto",
textsize = "11px",
opacity = 0.8
)
) %>%
addLegend(pal = palredgreen, values = values(slope.geo_dg),
title = "Pendiente del Terreno en Huila (grados)")10. visualización de la orientación del terreno
Ahora visualizaremos la orientación en grados, entendamos orientación como la dirección en la que apunta una pendiente, el proceso es relativamente similar al del punto anterior, ya que la orientación se mide en grados
(munic$mean_aspect_dg <- exactextractr::exact_extract(aspect, munic, 'mean'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 197.7818 181.6022 151.4429 156.6729 205.3990 187.4814 200.9467 218.4144
## [9] 222.9548 176.1548 231.3761 222.6475 212.0125 224.6192 172.8979 165.3855
## [17] 178.7280 164.5075 234.4746 165.4130 177.1659 158.6518 172.8113 150.1872
## [25] 182.9868 268.0864 158.3165 160.7825 158.3753 188.8365 148.3418 161.0296
## [33] 213.5841 170.9169 202.2813 215.3844 150.4501Obsevamos el histograma de la orientación del departamento deL Huila
hist(munic$mean_aspect_dg,
main = "Orientacion media del Huila",
xlab = "Orientacion (en grados)")
Transformamos los poligonos en un raster.
(munic$class_aspdeg <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4Reprojectamos la capa para que este en nuestro sistema de coordenadas de referencia.
(slope.geo_aspdg = project(aspect, "EPSG:4326"))## class : SpatRaster
## size : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : aspect
## min value : 0.4234968
## max value : 359.7618713la paleta de colores se realizo emulando una paleta observada en ArcGis que esta citada en las referencias. Para terminos prácticos se le coloco anotherpal
anotherpal <- colorNumeric(c("white","brown", "yellow2", "green3", "lightyellow", "lightgreen"), values(slope.geo),
na.color = "transparent")Como resultado obtenemos un mapa del departamento del Huila con sus respectivas orientaciones en grados
leaflet(munic) %>% addTiles() %>% setView(-75.5, 2.6, 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_dg, "<br>")) %>%
addRasterImage(slope.geo_aspdg, colors = anotherpal, opacity = 0.8) %>%
addLabelOnlyMarkers(
data = coords_jitter,
lng = ~lng, lat = ~lat,
label = ~label,
labelOptions = labelOptions(
noHide = TRUE,
direction = "auto",
textsize = "11px",
opacity = 0.8
)
) %>%
addLegend(pal = anotherpal, values = values(slope.geo_dg),
title = "Orientacion del Terreno en Huila (grados)")Referencias
Lizarazo, I., 2025. Geomorphometric terrain attributes in R. Available at https://rpubs.com/ials2un/geomorphometric
Esri. (n.d.). How Aspect works. ArcGIS Pro. https://pro.arcgis.com/en/pro-app/latest/tool-reference/spatial-analyst/how-aspect-works.htm
sessionInfo()## R version 4.5.0 (2025-04-11 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26100)
##
## 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