Calculo de variable geomorfometricas en el departamento del Amazonas
Sara Aponte
2025.11.04
1. Introduccion:
En este cuaderno se presenta la forma en que se calculó los atributos geomorfométricos de un terreno a partir de modelos digitales de elevación en un formato DEM.
2.Configuracion
Como paso inicial de limpia la memoria:
rm(list=ls())Lo siguiente la instalación y el uso de librerías.
library(sp)## Warning: package 'sp' was built under R version 4.5.2library(sf)## Warning: package 'sf' was built under R version 4.5.2## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUElibrary(leaflet)## Warning: package 'leaflet' was built under R version 4.5.2library(terra)## Warning: package 'terra' was built under R version 4.5.2## terra 1.8.70library(MultiscaleDTM)## Warning: package 'MultiscaleDTM' was built under R version 4.5.2library(exactextractr)## Warning: package 'exactextractr' was built under R version 4.5.2library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.5.2library(tidyterra)## Warning: package 'tidyterra' was built under R version 4.5.2##
## Adjuntando el paquete: 'tidyterra'## The following object is masked from 'package:stats':
##
## filterCada librería cumple una función en el cuaderno: - Terra: Lectura y procesamiento de datos espaciales Ráster. - sf: Lectura y procesamiento de datos espaciales vectoriales. - leaflet: Permitió la creación de mapas móviles o interactivos. - MultiscaleDTM: Análisis de terreno en capas espaciales distintas aplicando DEM. - Exactextractr: Extrae valores ráster dentro de polígonos - sp: Manejo de datos espaciales en R proporciona clases y métodos para puntos, líneas, polígonos y cuadrículas. - ggplot2: La relación entre variables y el aspecto gráfico con las capas geométricas como puntos y líneas. - tidyterra: Mezcla terra con la proyección y visualización de efectos espaciales.
Introduccion a MultiscaleDTM:
Su función principal, SlpAsp(), calcula la pendiente (slope) y la exposición (aspect). Este método es una versión mejorada de los algoritmos clásicos que usan ventanas de 3×3 celdas para estimar cómo cambia la elevación entre una celda y sus vecinas Para estos cálculos se utilizan distintos métodos, según las celdas que se consideren alrededor de la celda central: Rook: usa solo las cuatro celdas que están arriba, abajo, a la izquierda y a la derecha. Es rápido y simple, pero no detecta variaciones diagonales. Queen: incluye las ocho celdas vecinas (cardinales y diagonales). Da resultados más completos y suaves, aunque tarda un poco más en procesar. Boundary: utiliza todas las celdas del borde de la ventana (por ejemplo, las 16 externas en una ventana de 5×5). Permite analizar patrones de relieve amplios, aunque puede suavizar los detalles locales.
3.Datos de entrada
En esta parte se necesitan cargar los datos ráster del departamento de Amazonas, por esto se usa la librería terra. Adicional a esto se comprueba la ruta para ubicar los archivos necesarios y se aplican y leen en el cuaderno.
getwd()## [1] "C:/Users/DELL/Desktop"A continuación, se asigna “dem” a la memoria de R y la función “rast” cumple la función para cargar el objeto de elevación.
(dem = terra::rast("Datos/COL_ref/COL_msk_alt.tif"))## class : SpatRaster
## size : 2136, 1800, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -81.8, -66.8, -4.3, 13.5 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source : COL_msk_alt.tif
## name : COL_msk_alt
## min value : -28
## max value : 5453Con la función aggregate() sirve para modificar la resolución espacial de una ráster, combinando celdas originales con celdas grandes. Esta modificación se asigna a un nuevo dem (dem2) pero con una resolución mas grande.
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Esto disminuye el total de celdas, reduciendo el tamaño del archivo para el facil procesamiento.
dem2## class : SpatRaster
## size : 1068, 900, 1 (nrow, ncol, nlyr)
## resolution : 0.01666667, 0.01666667 (x, y)
## extent : -81.8, -66.8, -4.3, 13.5 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : COL_msk_alt
## min value : -3.00
## max value : 5222.75file.exists("Datos/Amazonas Munic.gpkg")## [1] TRUEAhora con la librería sf a la variable munic se le asigna el gpkg que contiene los límites de los municipios de Amazonas:
(munic <- sf::st_read("Datos/Amazonas Munic.gpkg"))## Reading layer `cortado' from data source
## `C:\Users\DELL\Desktop\Datos\Amazonas Munic.gpkg' using driver `GPKG'
## Simple feature collection with 8 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.38985 ymin: -4.228429 xmax: -69.36835 ymax: 0.1603
## Geodetic CRS: WGS 84## Simple feature collection with 8 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.38985 ymin: -4.228429 xmax: -69.36835 ymax: 0.1603
## Geodetic CRS: WGS 84
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2
## 1 53 COL Colombia 1 Amazonas 1 El Encanto
## 2 53 COL Colombia 1 Amazonas 2 La Chorrera
## 3 53 COL Colombia 1 Amazonas 3 La Pedrera
## 4 53 COL Colombia 1 Amazonas 4 Leticia
## 5 53 COL Colombia 1 Amazonas 5 Mirití-Paraná
## 6 53 COL Colombia 1 Amazonas 6 Puerto Nariño
## 7 53 COL Colombia 1 Amazonas 7 Puerto Santander
## 8 53 COL Colombia 1 Amazonas 8 Tarapacá
## TYPE_2 ENGTYPE_2 NL_NAME_2 VARNAME_2
## 1 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 2 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 3 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 4 Municipio Municipality <NA> <NA>
## 5 Corregimiento Departamento Corregimiento Departamento <NA> Miriti Parana
## 6 Municipio Municipality <NA> <NA>
## 7 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 8 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## geom
## 1 MULTIPOLYGON (((-73.23018 -...
## 2 MULTIPOLYGON (((-71.9802 -1...
## 3 MULTIPOLYGON (((-72.61586 -...
## 4 MULTIPOLYGON (((-70.0579 -4...
## 5 MULTIPOLYGON (((-70.2188 -0...
## 6 MULTIPOLYGON (((-70.28952 -...
## 7 MULTIPOLYGON (((-71.4002 -0...
## 8 MULTIPOLYGON (((-70.1296 -3...Ahora se deben alinear los datos en coordenadas geográficas que coincidan, para proceder a eso se debe anteceder con consultar los crs.
(Amazonas_crs = crs(munic))## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"crs(dem2)## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"crs(dem2) <- Amazonas_crs(crs(dem2))## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n MEMBER[\"World Geodetic System 1984 (G2296)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"Se debe recortar los limites del Amazonas al dem2 creando la nueva capa dem3
dem3 = terra::crop(dem2,munic, mask=TRUE)4.Transformar coordenadas
Se calculan los atributos del terreno, pasar el DEM a coordenadas planas para transformar los datos de entrada al sistema Colombiano
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 264, 302, 1 (nrow, ncol, nlyr)
## resolution : 1848.502, 1848.502 (x, y)
## extent : 4846116, 5404364, 1088497, 1576501 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : COL_msk_alt
## min value : 60.75
## max value : 306.21(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 8 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4845400 ymin: 1090195 xmax: 5404114 ymax: 1575767
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2
## 1 53 COL Colombia 1 Amazonas 1 El Encanto
## 2 53 COL Colombia 1 Amazonas 2 La Chorrera
## 3 53 COL Colombia 1 Amazonas 3 La Pedrera
## 4 53 COL Colombia 1 Amazonas 4 Leticia
## 5 53 COL Colombia 1 Amazonas 5 Mirití-Paraná
## 6 53 COL Colombia 1 Amazonas 6 Puerto Nariño
## 7 53 COL Colombia 1 Amazonas 7 Puerto Santander
## 8 53 COL Colombia 1 Amazonas 8 Tarapacá
## TYPE_2 ENGTYPE_2 NL_NAME_2 VARNAME_2
## 1 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 2 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 3 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 4 Municipio Municipality <NA> <NA>
## 5 Corregimiento Departamento Corregimiento Departamento <NA> Miriti Parana
## 6 Municipio Municipality <NA> <NA>
## 7 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## 8 Corregimiento Departamento Corregimiento Departamento <NA> <NA>
## geom
## 1 MULTIPOLYGON (((4974409 136...
## 2 MULTIPOLYGON (((5113382 135...
## 3 MULTIPOLYGON (((5042723 145...
## 4 MULTIPOLYGON (((5326564 110...
## 5 MULTIPOLYGON (((5309465 150...
## 6 MULTIPOLYGON (((5300930 113...
## 7 MULTIPOLYGON (((5177945 145...
## 8 MULTIPOLYGON (((5318898 119...En este caso se cambio el sistema crs al sistema nacional de coordenadas para los datos del gpkg asignando este cambio a la variable “munic_plane”
5.Calculo atributos del terreno
Para este paso, necesitamos calcular atributos del terreno —específicamente pendiente (slope) y exposición (aspect)— usando la función SlpAsp del paquete MultiscaleDTM en R.
(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 : 264, 302, 2 (nrow, ncol, nlyr)
## resolution : 1848.502, 1848.502 (x, y)
## extent : 4846116, 5404364, 1088497, 1576501 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0008076037, 6.205885e-03
## max values : 1.9861943743, 3.599911e+02Para este bloque el numero 1 indica la primera capa (pendiente) para atribuirla a (slope)
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 264, 302, 1 (nrow, ncol, nlyr)
## resolution : 1848.502, 1848.502 (x, y)
## extent : 4846116, 5404364, 1088497, 1576501 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0008076037
## max value : 1.9861943743La función “hist()” de la libreria terra, va a generar un histograma con la distribución de valores de pendiente en el área del Amazonas.Donde el eje X va a contener la pendiente en grados:
terra::hist(slope,
main = "Pendientes de Amazonas",
xlab = "Pendiente (en grados)")
El eje horizontal (X) se representan los grados de inclinación del terreno, mientras que en el eje vertical (Y) se muestra la frecuencia, es decir, el número de celdas o áreas que presentan cada rango de pendiente.Interpretando el grafico es correcto decir y afirmar que las pendientes del territorio del AMamzoans son pedominantemente bajas.
Para este siguiente bloque el 2 indica la capa de exposicion y se asigna a una variable llamada “aspect”
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 264, 302, 1 (nrow, ncol, nlyr)
## resolution : 1848.502, 1848.502 (x, y)
## extent : 4846116, 5404364, 1088497, 1576501 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 6.205885e-03
## max value : 3.599911e+02Ahora visualizaremos la distribucion de exposion en una tabla, con la funcion hist()
terra::hist(aspect,
main = "Exposición de Amazonas",
xlab = "Exposición (en grados)")
En este grafico las concentraciones de elevacion se ubican entre el sureste y suroeste del departamento.
Adicionalmente, convirtamos los grados de pendiente a porcentaje de pendiente, por ello con la función “tan(slope(pi/180))100” Calcula la tangente del ángulo de pendiente en radianes. Matemáticamente, la tangente de un ángulo (en este caso, la pendiente) representa la razón entre la altura vertical y la distancia horizontal.
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 264, 302, 1 (nrow, ncol, nlyr)
## resolution : 1848.502, 1848.502 (x, y)
## extent : 4846116, 5404364, 1088497, 1576501 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.001409534
## max value : 3.467952406Una vez a porcentaje lo podemos transformar en un grafico:
terra::hist(slope_perc,
main = "Pendientes de Amazonas",
xlab = "Pendiente (en porcentaje)")
Esta grafica indica que el reporte de pendiente corresponte a un porcentaje minimo.
6. Calculo de estadisticas
Ahora tenemos que calcular el promedio de los valores de pendiente por municipio. Para ello, veamos la clasificación de pendientes del IGAC para fines agronómicos:
La tabla muestra la relación entre el código (alfabetico), la clase de pendiente y el porcentaje de inclinación del terreno (%), junto con los valores de color empleados para su representación. Cada código, identificado por una letra de la a a la g, corresponde a un rango específico de pendiente y a una descripción cualitativa del relieve.
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)La matriz tiene como objetivo clasificar un modelo de pendientes expresadas en porcentaje en distintas categorías o clases, según la tabla de clasificación para usos agricolas. En primer lugar, se crea un vector denominado “m” que contiene los límites de cada rango estos son inferiones o superiores de la perndiente, junto con un valor numérico que representa el código asignado a cada clase. Luego este vector se transforma en una matriz mediante la función “matrix(m, ncol=3, byrow = TRUE”
En el siguiente bloque de codigo se utiliza una nueva variable a la que se introduce el corte de esta matriz con el raster inicial.
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%## [1] 0.2329995 0.2542693 0.2813751 0.3035096 0.3104131 0.2927065 0.3867155
## [8] 0.2565996Este ultimo dato ahora es una variable de munic y replesenta la moda de pendiente en los municipios del departamento con esto se va a generar otro histograma con la representacion media del relieve
hist(munic$mean_slope,
main = "Amazonas pendiente promedio",
xlab = "Pendiente (en grados)")
El resultado de esta proyeccion refuerza los resultados anteriores la presencia de inclinaciones fuertes en el departamento y los diferentes municipios es muy bajom, las pendiente promedio estan entre 0.25 y 0.30.
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========= | 12% | |================== | 25% | |========================== | 38% | |=================================== | 50% | |============================================ | 62% | |==================================================== | 75% | |============================================================= | 88% | |======================================================================| 100%## [1] 1 1 1 1 1 1 1 1Ahora generamos el histograma que muestra las pendientes por categorias, categorais que generamos con el codigo anterior realizando la extraccion de los valores raster como como lineas y poligonos.
hist(munic$class,
main = "Pendiente reclasificada de Amazonas",
xlab = "Pendiente (como categoría)")
Transformaremos la pendiente a coordenadas geograficas. La primera instrucción crea un nuevo ráster llamado slope.geo, que contiene los valores de pendiente en porcentaje (slope_perc) pero ahora expresados dentro del sistema de coordenadas geográficas, la linea rc.geo corresponde a la version reproyectada.
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 265, 302, 1 (nrow, ncol, nlyr)
## resolution : 0.01666665, 0.01666665 (x, y)
## extent : -74.38712, -69.35379, -4.249657, 0.1670042 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1.000000
## max value : 1.873985(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 265, 302, 1 (nrow, ncol, nlyr)
## resolution : 0.01666665, 0.01666665 (x, y)
## extent : -74.38712, -69.35379, -4.249657, 0.1670042 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.008331007
## max value : 3.1985640537.Pendiente del terreno proyectada
Para dibujar la pendiente del terreno en un mapa, primero se necesita
una paleta de colores. La función
colorNumeric() del paquete leaflet crea
una escala continua de colores para los valores del ráster
slope.geo, que muestra la pendiente en grados o
porcentaje.
Los colores van del verde (zonas planas o suaves) al rojo (pendientes
fuertes o zonas escarpadas). El argumento
na.color = "transparent" hace que las áraes sin valores no
se reflejen en la proyeccion del mapa.
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")El siguiente fragmento de código en R utiliza la librería Leaflet para crear un mapa interactivo:
leaflet(munic) %>% addTiles() %>% setView(-74.0, 4.7, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$NAME_2, "<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Pendientes en el depatamento Amazonas(%)")8.Alternativa de proyeccion
La instrucción leaflet(munic) crea el mapa tomando como
base el objeto munic, que contiene los polígonos de los municipios de
Cundinamarca. Luego, addTiles() añade una capa base de
OpenStreetMap, que sirve como referencia geográfica.
El comando setView(-74.0, 4.7, 9) centra el mapa cerca
del punto medio del departamento y establece un nivel de zoom general
(9).
Después, addPolygons() dibuja los límites municipales
con bordes grises y ajusta aspectos visuales como el grosor
(weight = 1.0), el suavizado
(smoothFactor = 0.5) y la transparencia
(opacity y fillOpacity).
Finalmente, se agrega unpopup interactivo que muestra información sobre cada municipio cuando se hace clic en él.
ggplot() +
geom_spatraster(data = slope) +
scale_fill_whitebox_c(
palette = "muted",
labels = scales::label_number(suffix = "°"),
n.breaks = 12,
guide = guide_legend(reverse = TRUE)
) +
geom_sf(data = munic_plane, fill = NA, color = "black", linewidth = 0.3) +
labs(
fill = "Pendiente (°)",
title = "Pendientes en el departamento de Amazonas"
) +
theme_bw()
El siguiente código genera un mapa de pendientes del departamento del Amazonas mediante la librería ggplot2.
La función geom_spatraster() representa los valores del
ráster slope, que contiene la pendiente del terreno en
grados. Para facilitar la interpretación visual, se aplica la escala de
color scale_fill_whitebox_c() con la paleta muted, que
muestra un gradiente de tonos suaves donde los colores más claros
representan pendientes bajas y los más oscuros, pendientes más
pronunciadas.
La capa geom_sf() superpone los límites municipales
contenidos en el objeto munic_plane, sin relleno y con
bordes negros finos (linewidth = 0.3), delimitando
claramente cada municipio dentro del mapa.
Finalmente, se añaden etiquetas y título con labs(),
especificando la variable representada (“Pendiente (°)”) y el título del
mapa. El tema theme_bw() otorga un fondo limpio y uniforme,
mejorando la visualización general del ploteo.
9.Vizualizar orientacion del proyecto
Ahora se calcula el visual del terreno a partir del modelo digital de
elevación dem3.El resultado, guardado en
aspect_AMAZONAS, muestra la dirección hacia la que se
orienta el relieve en la zona.
aspect_AMAZONAS <- terra::terrain(dem3)rgb.palette <- colorRampPalette(
c("blue", "cyan", "green", "yellow", "orange", "red", "magenta", "blue")
)
spplot(
aspect_AMAZONAS,
main = "Exposición del terreno en grados - Departamento de Amazonas",
col.regions = rgb.palette(100),
sp.layout = list(
list("sp.polygons", as_Spatial(munic), col = "black", lwd = 1.2)
),
xlim = c(ext(munic)[1], ext(munic)[2]),
ylim = c(ext(munic)[3], ext(munic)[4])
)
El código crea un mapa del aspecto del terreno en el Departamento de
Amazonas. Primero, se define una paleta de colores
rgb.palette que va desde azul hasta magenta para
representar diferentes valores de orientación de la pendiente. Luego,
con spplot() se grafica aspect_AMAZONAS,
aplicando esta paleta de colores, añadiendo los límites de los
municipios (munic) y ajustando el mapa al área de interés
mediante xlim y ylim. El resultado es un mapa
que muestra la dirección de las pendientes del terreno de manera visual
y clara.
En este caso la orientacion de la pendiente esta hacia norte del departamento.
Bibliografia
Libro guia:
Lizarazo, I., 2025. Geomorphometric terrain attributes in R. Available at https://rpubs.com/ials2un/geomorphometric
Librerias:
Baston, D. (2023). exactextractr: Fast extraction from raster datasets using polygons. R package version 0.10.1. Disponible en: https://CRAN.R-project.org/package=exactextractr Bivand, R. S., Pebesma, E., & Gómez-Rubio, V. (2013). Applied spatial data analysis with R (2nd ed.). Springer. Cheng, J., Karambelkar, B., & Xie, Y. (2023). leaflet: Create interactive web maps with the JavaScript Leaflet library. R package version 2.2.2. Disponible en: https://CRAN.R-project.org/package=leaflet Hijmans, R. J. (2024). terra: Spatial data analysis. R package version 1.7-78. Disponible en: https://CRAN.R-project.org/package=terra Pebesma, E. (2018). Simple features for R: Standardized support for spatial vector data. The R Journal, 10(1), 439–446. Waga, K., Drăguţ, L., Eisank, C., & Blaschke, T. (2022). MultiscaleDTM: Multiscale analysis of digital terrain models. R package version 1.0. Disponible en: https://CRAN.R-project.org/package=MultiscaleDTM Zia, M. (2023). Digital Terrain Analysis. En Digital Soil Mapping with R. Recuperado de https://zia207.quarto.pub/digital-terrain-analysis.html
sessionInfo()## R version 4.5.1 (2025-06-13 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] tidyterra_0.7.2 ggplot2_4.0.0 exactextractr_0.10.0
## [4] MultiscaleDTM_1.0 terra_1.8-70 leaflet_2.2.3
## [7] sf_1.0-21 sp_2.2-0
##
## loaded via a namespace (and not attached):
## [1] tidyr_1.3.1 sass_0.4.10 generics_0.1.4 class_7.3-23
## [5] KernSmooth_2.23-26 lattice_0.22-7 digest_0.6.37 magrittr_2.0.4
## [9] rgl_1.3.24 RColorBrewer_1.1-3 evaluate_1.0.5 grid_4.5.1
## [13] fastmap_1.2.0 jsonlite_2.0.0 e1071_1.7-16 DBI_1.2.3
## [17] promises_1.5.0 purrr_1.1.0 scales_1.4.0 crosstalk_1.2.2
## [21] codetools_0.2-20 jquerylib_0.1.4 cli_3.6.5 shiny_1.11.1
## [25] rlang_1.1.6 units_1.0-0 withr_3.0.2 base64enc_0.1-3
## [29] cachem_1.1.0 yaml_2.3.10 otel_0.2.0 raster_3.6-32
## [33] tools_4.5.1 dplyr_1.1.4 httpuv_1.6.16 png_0.1-8
## [37] vctrs_0.6.5 R6_2.6.1 mime_0.13 proxy_0.4-27
## [41] lifecycle_1.0.4 classInt_0.4-11 htmlwidgets_1.6.4 pkgconfig_2.0.3
## [45] gtable_0.3.6 bslib_0.9.0 pillar_1.11.1 later_1.4.4
## [49] data.table_1.17.8 glue_1.8.0 Rcpp_1.1.0 xfun_0.53
## [53] tibble_3.3.0 tidyselect_1.2.1 rstudioapi_0.17.1 knitr_1.50
## [57] farver_2.1.2 xtable_1.8-4 htmltools_0.5.8.1 labeling_0.4.3
## [61] rmarkdown_2.30 compiler_4.5.1 S7_0.2.0Realizado por: SAM