Este cuaderno tiene como objetivo ilustrar el cálculo de atributos geomorfométricos del terreno a partir de Modelos Digitales de Elevación (MDE) cuadriculados, utilizando el paquete MultiscaleDTM en R. Esta herramienta permite derivar indicadores esenciales como pendiente, orientación (aspect), curvatura, posición relativa y rugosidad, aplicándolos a múltiples escalas espaciales (Ilich et al., 2023). Los atributos geomorfométricos son fundamentales en diversas disciplinas:son clave para entender la hidrología del terreno, la erosión, y la distribución espacial de propiedades del suelo (Franklin, 2020; Kokulan et al., 2018). En contextos agrícolas, la topografía—especialmente la pendiente y la orientación—afecta el flujo de agua, la exposición solar y la variabilidad de nutrientes del suelo, impactando directamente el rendimiento de los cultivos (Geopard Tech, 2020; 2022). El área de estudio será el departamento del Tolima, seleccionado por su relevancia agrícola y su diversidad de relieves, lo que lo convierte en un escenario ideal para ilustrar la relación entre geomorfometría y producción agropecuaria.
Los paquetes necesarios deben instalarse previamente desde la consola. Paquetes ncesarios = c(“MultiscaleDTM”, “exactextractr”) con la funcion install.packages(paquetes)
Se debe empezar por limpiar la memoria: Esto se hace para comenzar un script desde cero, para evitar que datos anteriores interfieran y asegurar que el código funcione con un entorno limpio (por ejemplo, antes de compartirlo o reproducirlo).
rm(list=ls())
Ahora, se cargan las bibliotecas necesarias:
library(elevatr)
## elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'. Use
## of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being
## deprecated; however, get_elev_raster continues to return a RasterLayer. This
## will be dropped in future versions, so please plan accordingly.
library(sf)
## Linking to GEOS 3.13.1, GDAL 3.10.2, PROJ 9.5.1; sf_use_s2() is TRUE
library(leaflet)
library(terra)
## terra 1.8.54
library(MultiscaleDTM)
library(exactextractr)
package_name::function.
El paquete MultiscaleDTM en R permite calcular atributos geomorfométricos multiescala—como pendiente, orientación, curvatura, posición relativa y rugosidad—utilizando modelos digitales del terreno (MDT), es decir, rásteres de elevación con cuadrícula regular, aplicando una ventana focal definida (Ilich et al., 2023). sta herramienta es particularmente útil porque permite evaluar las propiedades del terreno a diferentes escalas espaciales, una capacidad que muchos paquetes no ofrecen (Ilich et al., 2023). El presente cuaderno tiene por objetivo introducir brevemente el uso de MultiscaleDTM, orientado a estudiantes de pregrado en agronomía de la Universidad Nacional de Colombia. Se aplicará específicamente al departamento del Tolima, seleccionado por su riqueza topográfica y relevancia agrícola, lo que lo convierte en un escenario adecuado para ejemplificar cómo estas variables geomorfológicas influyen en la dinámica del agua, la exposición solar y, potencialmente, la productividad del suelo.
En lo que sigue: - Se instalará y cargará el paquete MultiscaleDTM. - Se explicará cada fragmento de código. - Se describirá e interpretará cada resultado generado. - Se evaluará la multiescalaridad de los atributos en el contexto del Tolima.
Primero, se carga un DEM ráster para el departamento en estudio(Tolima). Para luego leer el DEM con terra:
(dem = terra::rast("C:/Users/menju/OneDrive/Documentos/GB2/Rstudio2/elev_ELEVACION_z10.tif"))
## class : SpatRaster
## size : 4092, 3078, 1 (nrow, ncol, nlyr)
## resolution : 0.000685414, 0.000685414 (x, y)
## extent : -76.28906, -74.17936, 2.811272, 5.615986 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_ELEVACION_z10.tif
## name : elev_ELEVACION_z10
## min value : -427
## max value : 5372
Reduzcamos la resolución de DEM para evitar problemas de memoria: Para reducir la resolución del DEM (hacerlo más “grueso”) y disminuir el uso de memoria y procesamiento al trabajar con grandes volúmenes de datos.
dem2 = terra::aggregate(dem,2, "mean")
## |---------|---------|---------|---------|=========================================
Se leen los municipios del Tolima utilizando la biblioteca sf:
(munic<- sf::st_read("C:/Users/menju/OneDrive/Documentos/GB2/P1/DATOS/MUNI_TOL.gpkg"))
## Reading layer `MUNICIPIOS_TOL' from data source
## `C:\Users\menju\OneDrive\Documentos\GB2\P1\DATOS\MUNI_TOL.gpkg'
## using driver `GPKG'
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS: MAGNA-SIRGAS
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 73 001 73001 TOLIMA IBAGUÉ
## 2 73 024 73024 TOLIMA ALPUJARRA
## 3 73 026 73026 TOLIMA ALVARADO
## 4 73 030 73030 TOLIMA AMBALEMA
## 5 73 043 73043 TOLIMA ANZOÁTEGUI
## 6 73 055 73055 TOLIMA ARMERO
## 7 73 067 73067 TOLIMA ATACO
## 8 73 124 73124 TOLIMA CAJAMARCA
## 9 73 148 73148 TOLIMA CARMEN DE APICALÁ
## 10 73 152 73152 TOLIMA CASABIANCA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1906 MUNICIPIO 1377.2514 2023
## 2 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 501.5113 2023
## 3 1540 MUNICIPIO 344.3356 2023
## 4 1627 MUNICIPIO 238.0376 2023
## 5 Ordenanza 21 del 30 de Marzo de 1915 MUNICIPIO 469.7113 2023
## 6 Decreto 1049 de Septiembre 29 de 1908 MUNICIPIO 440.6308 2023
## 7 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 1017.6598 2023
## 8 1913 MUNICIPIO 508.4909 2023
## 9 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 190.8361 2023
## 10 Ordenanza 26 del 22 de Junio de 1896 MUNICIPIO 174.9306 2023
## shape_Leng shape_Area codigo geom
## 1 2.4800420 0.11217106 NA MULTIPOLYGON (((-75.36517 4...
## 2 1.0547137 0.04080342 NA MULTIPOLYGON (((-74.99192 3...
## 3 1.0695888 0.02805438 NA MULTIPOLYGON (((-74.97668 4...
## 4 1.0091029 0.01940220 NA MULTIPOLYGON (((-74.74744 4...
## 5 1.3101775 0.03826692 NA MULTIPOLYGON (((-75.20098 4...
## 6 1.4441004 0.03592418 NA MULTIPOLYGON (((-74.88687 5...
## 7 2.2522480 0.08276960 NA MULTIPOLYGON (((-75.3251 3....
## 8 1.0293681 0.04140402 NA MULTIPOLYGON (((-75.48219 4...
## 9 0.7339295 0.01554073 NA MULTIPOLYGON (((-74.74471 4...
## 10 1.0803200 0.01425931 NA MULTIPOLYGON (((-75.07491 5...
Se debe tener en cuenta que ambos conjuntos de datos están en coordenadas geográficas.
Ahora, se recorta el objeto dem2 a los límites de nuestro departamento:
# 1. Convertir sf a SpatVector
munic_vect <- terra::vect(munic)
# 2. Reproyectar al CRS del raster
munic_vect <- terra::project(munic_vect, dem2)
# 3. Recortar el DEM usando el objeto reproyectado
dem3 <- terra::crop(dem2, munic_vect, mask = TRUE)
El Instituto Geográfico Agustín Codazzi (IGAC) adoptó en el año 2020 un nuevo sistema de coordenadas planas de origen único nacional para todo el territorio colombiano. El código EPSG correspondiente a este sistema es el EPSG:9377, conocido como MAGNA-SIRGAS / Colombia Bogotá.
Este sistema de referencia se fundamenta en la proyección Transversa de Mercator con origen en Bogotá D.C., y pertenece al marco geodésico MAGNA-SIRGAS, el cual es parte de la Red Geodésica Nacional de Colombia. Esta adopción responde a la necesidad de contar con un sistema de coordenadas unificado que facilite la interoperabilidad geoespacial, elimine la fragmentación entre zonas UTM y mejore la precisión en el análisis de información geográfica, especialmente en modelos digitales del terreno, cartografía base y catastros.
Importancia del uso del EPSG:9377 en análisis de terreno El cálculo de atributos topográficos como pendiente, orientación (aspecto), curvatura y flujo superficial requiere que el Modelo Digital de Elevación (MDE) esté en un sistema de coordenadas proyectadas, es decir, en unidades métricas (como metros), no geográficas (grados). Esto es crucial, ya que las operaciones de derivación espacial (como gradientes o transformaciones de superficie) dependen de distancias y áreas reales.
Trabajar con EPSG:9377 garantiza resultados consistentes a lo largo del territorio colombiano sin necesidad de subdividir el país por zonas UTM, como se hacía anteriormente. Además, favorece la integración con sistemas oficiales como el Catastro Multipropósito, la Infraestructura Colombiana de Datos Espaciales (ICDE), y los datos de planificación territorial y ambiental.
# check the parameters of this terra function here
# https://rdrr.io/github/rspatial/terra/man/project.html
(dem_plane = project(dem3, "EPSG:9377"))
## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_ELEVACION_z10
## min value : -1.783039
## max value : 5277.944824
También es necesario realizar una transformación similar para el objeto vectorial:
# check the parameters of this sf function here
# https://github.com/rstudio/cheatsheets/blob/main/sf.pdf
(munic_plane = sf::st_transform(munic, "EPSG:9377"))
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4654882 ymin: 1875660 xmax: 4836359 ymax: 2146078
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 73 001 73001 TOLIMA IBAGUÉ
## 2 73 024 73024 TOLIMA ALPUJARRA
## 3 73 026 73026 TOLIMA ALVARADO
## 4 73 030 73030 TOLIMA AMBALEMA
## 5 73 043 73043 TOLIMA ANZOÁTEGUI
## 6 73 055 73055 TOLIMA ARMERO
## 7 73 067 73067 TOLIMA ATACO
## 8 73 124 73124 TOLIMA CAJAMARCA
## 9 73 148 73148 TOLIMA CARMEN DE APICALÁ
## 10 73 152 73152 TOLIMA CASABIANCA
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 1906 MUNICIPIO 1377.2514 2023
## 2 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 501.5113 2023
## 3 1540 MUNICIPIO 344.3356 2023
## 4 1627 MUNICIPIO 238.0376 2023
## 5 Ordenanza 21 del 30 de Marzo de 1915 MUNICIPIO 469.7113 2023
## 6 Decreto 1049 de Septiembre 29 de 1908 MUNICIPIO 440.6308 2023
## 7 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 1017.6598 2023
## 8 1913 MUNICIPIO 508.4909 2023
## 9 Decreto 650 del 13 de Octubre de 1887 MUNICIPIO 190.8361 2023
## 10 Ordenanza 26 del 22 de Junio de 1896 MUNICIPIO 174.9306 2023
## shape_Leng shape_Area codigo geom
## 1 2.4800420 0.11217106 NA MULTIPOLYGON (((4737717 207...
## 2 1.0547137 0.04080342 NA MULTIPOLYGON (((4778805 194...
## 3 1.0695888 0.02805438 NA MULTIPOLYGON (((4780825 207...
## 4 1.0091029 0.01940220 NA MULTIPOLYGON (((4806322 210...
## 5 1.3101775 0.03826692 NA MULTIPOLYGON (((4755947 207...
## 6 1.4441004 0.03592418 NA MULTIPOLYGON (((4790926 212...
## 7 2.2522480 0.08276960 NA MULTIPOLYGON (((4741831 196...
## 8 1.0293681 0.04140402 NA MULTIPOLYGON (((4724687 206...
## 9 0.7339295 0.01554073 NA MULTIPOLYGON (((4806435 202...
## 10 1.0803200 0.01425931 NA MULTIPOLYGON (((4770073 212...
Se procede a calcular la pendiente y la orientación.
# Explain what is w
# Explain what is method
# Change if needed
(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 : 1787, 1197, 2 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0000, 7.917609e-04
## max values : 69.1466, 3.599997e+02
Ahora, dividamos el objeto slp_asp en dos partes.
Capa 1:
(slope = subset(slp_asp, 1))
## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 69.1466
¿Cuál es la distribución de valores de pendiente?
terra::hist(slope,
main = "Tolima's slope",
xlab = "Slope (in degrees)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
El histograma mostrado representa la distribución de la pendiente del
terreno (en grados) en el departamento del Tolima, Colombia. Se observa
una clara asimetría hacia la derecha, lo que indica que la mayoría del
área analizada presenta pendientes suaves. Específicamente, la mayor
concentración de valores se encuentra entre los 0° y 10°, con una
frecuencia que supera los 120,000 píxeles, lo cual sugiere que una gran
proporción del terreno es relativamente plano o ligeramente
inclinado.
A medida que aumenta el valor de la pendiente, la frecuencia disminuye de forma progresiva. Pendientes mayores a 30° son considerablemente menos frecuentes, y los valores extremos (mayores a 50°) son casi inexistentes, lo que confirma que las áreas abruptas son escasas en la región.
Este patrón es típico en territorios con una combinación de zonas planas agrícolas (como el valle del Magdalena) y zonas montañosas en la cordillera central, reflejando la heterogeneidad del relieve tolimense. Este tipo de análisis es crucial en estudios de aptitud agrícola, riesgo de erosión y ordenamiento territorial, ya que la pendiente influye directamente en la capacidad de uso del suelo y en procesos como escorrentía o retención de agua.
(aspect =subset(slp_asp, 2))
## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 7.917609e-04
## max value : 3.599997e+02
terra::hist(aspect,
main = "Tolima aspect",
xlab = "Aspect (in degrees)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
El histograma presentado muestra la distribución del aspecto del terreno
(orientación de las pendientes) en el departamento del Tolima, expresado
en grados desde 0° (Norte) hasta 360° (también Norte, cerrado el
círculo). Este atributo geomorfométrico es crucial para comprender la
exposición solar, la circulación de aire, y el microclima local, todos
factores determinantes en procesos ecológicos, agrícolas y de
erosión.
En la gráfica se observa una distribución relativamente equilibrada, aunque con una mayor frecuencia de aspectos entre los 90° y 150°, lo cual corresponde a pendientes orientadas hacia el Este y Sureste. Este patrón puede deberse a la influencia de la cordillera Central y la topografía del valle del río Magdalena, donde las formas del terreno tienden a inclinarse en esa dirección.
Se destaca una ligera disminución de frecuencias entre los 180° y 270°, es decir, en las orientaciones Sur a Oeste, posiblemente porque esas laderas son menos predominantes en las unidades de relieve analizadas. Por otra parte, existe una recuperación hacia el rango de 270° a 360°, sugiriendo una mayor presencia de laderas con exposición Noroeste.
Este tipo de análisis es útil para:
Planificación agrícola, ya que las orientaciones determinan la radiación solar disponible.
Gestión de cuencas hidrográficas, pues el aspecto influye en el derretimiento de nieve (en zonas de alta montaña) o en la dirección del escurrimiento superficial.
Ecología del paisaje, al afectar la distribución de la vegetación, humedad del suelo y condiciones microclimáticas.
Además, se convierten los grados de pendiente en porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)
## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 262.5146
terra::hist(slope_perc,
main = "Tolima's slope",
xlab = "Slope (in percentage)")
## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
El histograma presentado muestra la distribución de la pendiente del
terreno en porcentaje (%) para el departamento del Tolima, calculada a
partir de un Modelo Digital de Elevación (MDE). Esta variable es clave
para analizar procesos erosivos, escorrentía, accesibilidad del terreno
y aptitud agrícola.
Se observa una distribución fuertemente sesgada hacia la izquierda, con la mayor parte de las pendientes por debajo del 30%. Las clases de pendiente con valores bajos (0–20%) son las más frecuentes, superando ampliamente las 120.000 observaciones en los primeros intervalos. Esto indica que una proporción significativa del territorio del Tolima es relativamente plano o con pendiente suave, lo cual es favorable para actividades agrícolas y de infraestructura.
A medida que aumentan los valores de pendiente, la frecuencia disminuye rápidamente. Las pendientes mayores al 100% (equivalente a un ángulo de inclinación de 45°) son escasas, y representan zonas de topografía abrupta, posiblemente correspondientes a áreas de montaña, como las zonas cercanas al Nevado del Tolima o la cordillera Central. Esta distribución sugiere que el uso del suelo y la ocupación humana en el Tolima están favorecidos por las amplias zonas de baja pendiente.mSin embargo, las zonas de alta pendiente, aunque menos frecuentes, deben ser cuidadosamente gestionadas para evitar riesgos de erosión, deslizamientos y para garantizar prácticas agrícolas sostenibles.
Una operación de estadísticas zonales calcula estadísticas sobre los valores de celda de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En este ncaso, nos interesa calcular el promedio de los valores de pendiente por municipio.
Primero, se conoce la clasificación de pendientes del IGAC para fines agrológicos:
(Reclacificar un Ráster)
# Reclassify the slope raster
#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
75, 160, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)
Ahora, se calculan 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% | |=== | 4% | |==== | 6% | |====== | 9% | |======= | 11% | |========= | 13% | |========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================== | 26% | |=================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |=============================== | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |======================================= | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |=================================================== | 72% | |==================================================== | 74% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 85% | |============================================================= | 87% | |=============================================================== | 89% | |================================================================ | 91% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%
## [1] 28.327206 22.369404 15.429619 3.669801 33.040310 9.219158 27.107445
## [8] 36.053734 15.664359 35.969097 28.613630 18.607035 7.297594 18.536375
## [15] 20.129711 1.252090 18.143921 1.507208 24.973742 2.159666 34.832108
## [22] 16.668184 22.115511 11.497709 32.145004 11.752955 16.978456 27.755978
## [29] 16.803055 18.587324 20.755070 8.394407 34.529945 15.837931 10.465549
## [36] 41.639145 27.062590 30.451038 1.194180 36.182438 11.381176 31.772017
## [43] 13.304363 17.442604 11.510288 35.072704 18.109522
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
hist(munic$mean_slope,
main = "Tolima's mean slope",
xlab = "Slope (in percentage)")
El histograma muestra la pendiente promedio (en porcentaje) de los
municipios del Tolima. Se observa que la mayoría de los municipios
tienen pendientes medias entre 15% y 25%, con un pico en el intervalo de
15–20%. Solo unos pocos municipios presentan pendientes promedio
superiores al 35%, lo que sugiere que en general el relieve municipal es
moderadamente inclinado.
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
## | | | 0% | |= | 2% | |=== | 4% | |==== | 6% | |====== | 9% | |======= | 11% | |========= | 13% | |========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================== | 26% | |=================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |=============================== | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |======================================= | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |=================================================== | 72% | |==================================================== | 74% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 85% | |============================================================= | 87% | |=============================================================== | 89% | |================================================================ | 91% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%
## [1] 5 4 1 1 5 1 5 5 1 5 5 4 1 4 4 1 4 1 5 1 5 4 4 1 5 4 4 5 1 5 4 1 5 4 1 5 4 5
## [39] 1 5 1 5 1 4 1 5 4
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)
hist(munic$class,
main = "Tolima's reclassified slope",
xlab = "Slope (as a category)")
El histograma muestra la pendiente reclasificada de los municipios del
Tolima en cinco categorías. Las clases 1 y
5 (pendientes muy bajas y muy altas, respectivamente)
son las más frecuentes, cada una con alrededor de 16
municipios. Las clases intermedias presentan menor frecuencia,
destacándose la clase 2, que no tiene representación.
Esto evidencia una distribución bimodal del relieve en el
departamento.
Transformar nuevamente la pendiente a coordenadas geográficas:
# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))
## class : SpatRaster
## size : 1793, 1199, 1 (nrow, ncol, nlyr)
## resolution : 0.001370805, 0.001370805 (x, y)
## extent : -76.11484, -74.47124, 2.868096, 5.325948 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.2429543
## max value : 259.9135742
#This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))
## class : SpatRaster
## size : 1793, 1199, 1 (nrow, ncol, nlyr)
## resolution : 0.001370805, 0.001370805 (x, y)
## extent : -76.11484, -74.47124, 2.868096, 5.325948 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1.254625e-04
## max value : 2.599136e+02
Primero, se seleciona una paleta de colores para la pendiente:
# expand the number of colors to improve your visualization
# https://r-graph-gallery.com/42-colors-names.html
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")
Ahora es el momento de trazar el gráfico. Asegúrate de cambiar varias opciones en el código si es necesario:
leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>%
addPolygons(color = "black", 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 Tolima (%)")
## 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'