Atributos de terreno geomorfológicos en R
Brandon Alejandro Salgado Garzón
23/07/2025
- 1) Introducción 🌍
- 2) configuración 🛠️
- 3) introducción a la libreria multiscaleDTM 🧩
- 4) leyendo datos de entrada 🧠
- 5) transformar coordenadas (10%) 🔧
- 6) Cálculo de atributos del terreno (20%) 🌱
- 7) Cálculo de estadísticas zonales (20%)
- 8) Trazando la pendiente del terreno (20%)
- 9) Visualización alternativa de la pendiente del terreno 15%
- 10) Bibliografia
- Reproducibilidad
1) Introducción 🌍
En este cuaderno se sigue el procedimiento para calcular datos geomorfometricos a partir de un modelo raster. Se va a utilizar la libreria de R “MultiscaleDTM”.
Se va a medir la pendiente y la Orientación del departamento del Magdalena con ayuda de esta libreria, esta información es necesaria en la agricultura para establecer la inclinación necesaria para que fluya el agua de manera adecuada y se eviten problemas de encharcamiento del terreno, la orientación de un terreno tambien es fundamental en el mundo agricola, se usa para planificar actividades de campo, ubicacion dentro del terreno, medir áreas y obtener informacion sobre elevaciones y distancias.
2) configuración 🛠️
Antes de iniciar debemos descargar las siguientes librerias para facilitar el manejo de los datos.
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.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.8.54library(MultiscaleDTM)
library(exactextractr)La libreria “elevatr” prporciona acceso a las redes de elevación de datos y transforma datos de elevación como un objeto de entidades simples “sf” desde servicios de elevación de puntos o como un objeto raster.
la libreria “leaflet” es una herramienta usada para crear mapas interactivos. Utiliza una interfaz para la biblioteca de javascrpit que permite consttruir visualizaciones espaciales dinamicas, suele usarse para mapas geograficos interactivos gracias a su adaptabilidad.
la libreria “terra” es una herramienta de datos espaciles, que maneja tanto dats vectoriales como raster. Proporciona metodos para operaciones como intersecciones, buffers y tareas de raster especificas como operaciones locales, globales, focales y zonales.
la libreria “exactextractr” es una herramienta diseñada para resumir de manera eficiente los valores de ráster sobre areas poligonales, este proceso se llama “estadisticas zonales”. Maneja casos en los que las celdas de la cuadricula están solo parcialmente cubiertas por poligonos, garantizando resultados precisos.
3) introducción a la libreria multiscaleDTM 🧩
La librería MultiscaleDTM es un paquete de R diseñado para el
análisis avanzado de modelos digitales del terreno (DTM/DEM) mediante el
cálculo de métricas topográficas a diferentes escalas espaciales. Su
principal fortaleza radica en permitir la extracción de atributos
geomorfológicos del terreno como pendiente, curvatura, rugosidad,
orientación, índice de posición topográfica (TPI), entre otros en
ventanas móviles de distintos tamaños, lo cual resulta clave para
caracterizar la variabilidad espacial del relieve. Esta capacidad
multiescalar es particularmente útil en disciplinas como la
geomorfología, hidrología, ecología del paisaje y agricultura de
precisión, donde el relieve del terreno influye en procesos como la
acumulación de agua, la erosión del suelo, el drenaje y el rendimiento
de cultivos. MultiscaleDTM opera sobre objetos SpatRaster del paquete
terra, permitiendo integrar fácilmente estos análisis con otros flujos
de trabajo geoespaciales en R. 
Funciones principales
- SlpAsp: Es usada para calcular la orientación y pendiente considerando un subconjunto de celdas dentrp de la ventana focal, específicamente las cuatro celdas en el borde de la ventana focal directamente arriba, abajo, izquierda y derecha de la ventana focal para el metodo de “torre”, cuatro celdas de esquina adicionales para el metodo de “reina” o todas las celdas de borde para el metodo de “limite”.
DirSlp: calcula la pendiente multiescala en una dirección especificada.
Qfit: Calcula la pendiente, orientación, curvaturay las caracteristicas morfométricas ajustando una superficie cuadratica a la ventana focal utilizando minimos cuadrados ordinarios utilizando la ecuación Z=unX^2 +bY^2 + cXY + dX + eY + f, donde a-f son los parámetros de regresión, Z es la elevación/profundidad, X son las coordenadas este/oeste en la ventana focal en relación con la celda focal e Y son las coordenadas norte/sur en la ventana focal en relación con la celda focal.
- Pfit: Calcula la pendiente multiescala y el aspecto de la orientación ajustando una superficie plana a la ventana focal utilizando mínimos cuadrados ordinarios. Proporciona resultados equivalentes al uso de un ajuste cuadratico, solo que exigiendo menos al computador.
¿Que es la posición relativa?
La posición relativa representa si un área es un alto o un mínimo local en relación con una altura de referencia. Se calcula como el valor de la celda focal menos el valor de una elevación de referencia, que a menudo es la media de los valores incluidos en la ventana focal, pero también podrían ser otras funciones, como la mediana, el mínimo o el máximo de los valores incluidos. Los valores positivos indican máximos topográficos locales y los valores negativos indican mínimos. La posición relativa se puede expresar en unidades del ráster DTM de entrada o se puede estandarizar en relación con la topografía local dividiendo por la desviación estándar o el rango de los valores de elevación incluidos en la ventana focal.
¿Como se puede usar para la elaboracion de mapas agricolas?
Esta libreria permite calcular atributos geomorfologicos del terreno a multiples escalas de modelos digitales de elevación DEM. Atributos como pendiente, curvatura, amplitud y entropia, resultan fundamentales para analizar condiciones del suelo como el drenaje, riesgo de erosion y variabilidad espacial, permitiendo de esta manera zonificar areas agricolas, optimizar practicas de riego y fertilización y la toma de decisiones importantes en la agricultura de presición. Se pueden generar mapas raster que describan microformas del terreno y caracteristicas topograficas que influyen directamente en la productividad del cultivo; además, estos pueden integrarse con capas vectoriales (como parcelas) usando herramientas como exactextractr para realizar análisis espaciales más detallados y específicos.
4) leyendo datos de entrada 🧠
Para empezar lo primero que debemos hacer es cargar al directorio la capa raster del magdalena, que fue previamente descargada.
# Cargamos el archivo .tif del DEM del departamento del Magdalena usando la función rast() de la libreria terra.
# El DEM es la base que se va a utilizar para calcular la pendiente y la orientación.
(dem = terra::rast("C:\\Users\\brand\\OneDrive\\Escritorio\\GB2\\P6\\MAGDALENA\\elv_magda_z10.tif"))## class : SpatRaster
## size : 4078, 2589, 1 (nrow, ncol, nlyr)
## resolution : 0.0006789017, 0.0006789017 (x, y)
## extent : -75.23438, -73.4767, 8.754526, 11.52309 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elv_magda_z10.tif
## name : elv_magda_z10
## min value : -2439
## max value : 5676Ahora se procede a reducir la resolución del DEM para usar menos memoria y acelerar el procesamiento, especialmente si el DEM es muy detalladado como en este caso.
#Reducimos la resolucion del raster original agregando los valores de las celdas en bloques de 2x2 y calculando su promedio.
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Una vez que reducimos su tamaño, podemos leer el geopackage de el departamento con sus municipios para mas adelante recortar el DEM y limitar el analisis a una región especifica.
# usamos la libreria sf para leer la información del geopackage del dpto. del Magdalena.
(munic <- sf::st_read("C:\\Users\\brand\\OneDrive\\Escritorio\\GB2\\P6\\MAGDALENA\\magdalena.gpkg"))## Reading layer `magdalena' from data source
## `C:\Users\brand\OneDrive\Escritorio\GB2\P6\MAGDALENA\magdalena.gpkg'
## using driver `GPKG'
## Simple feature collection with 30 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -74.9466 ymin: 8.936489 xmax: -73.54184 ymax: 11.34891
## Geodetic CRS: MAGNA-SIRGASRecortamos el DEM para quedarnos solo con el area del departamento, eliminando el resto del raster.
#Recortamos el raster dem2 usando el vector munic, y con mask = TRUE convierte en NA las areas fuera de los limites municipales.
dem3 = terra::crop(dem2,munic, mask=TRUE)## Warning: [crop] CRS do not match5) transformar coordenadas (10%) 🔧
Ahora pasamos el DEM a coordenadas planas con el sistema de coordenadas nacional, que es el EPSG:9377, corrspondiente a la proyección transversal MAGNA-SIRGAS/Origen nacional, el cual fue adoptado por el IGAC desde 2020.
¿porque se cambio? Antes, Colombia usaba el sistema UTM dividido en zonas (por ejemplo, zona 18N o 19N), lo cual fragmentaba el país en sectores con distintos orígenes. Esto dificultaba análisis geoespaciales integrales y la interoperabilidad de datos.
Mayor practicidad Al trabajar con DEM, parcelas agrícolas, coberturas, suelos, capas catastrales, etc., es altamente recomendable transformar todos los datos a EPSG:9377, ya que es el sistema de referencia aceptado por entidades oficiales como IGAC y el Instituto Geográfico Militar, y está optimizado para análisis geoespaciales a escala nacional.
Ahora transformamos los datos de entrada al sistema de coordenadas nacional.
#Con la función project() cambiamos el sistema de coordenadas al sistema MAGNA-SIRGAS.
#Con este sistema calcular la pendiente y la orientación es mas facil.
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elv_magda_z10
## min value : -158.9653
## max value : 5653.6406Ahora transformamos el archivo con los poligonons de los municipios al sistema MAGNA-SIRGAS . Es muy importante que el raster y el vector estén en el mismo sistema de referencia espacial para que puedan interactuar correctamente por ejemplo, si queremos recortar el raster con el shapefile, o extraer valores por zona.
#Con la función st_transform() transformamos el sistema de coordenadas al nacional y pasamos las coordenadas a coordenadas planas en metros.
(munic_plane = sf::st_transform(munic, "EPSG:9377"))6) Cálculo de atributos del terreno (20%) 🌱
Calculamos la pendiente y el aspecto multiescalar mediante la función SlpAsp, que calcula estas variables por medio de formulaciones clásicas de las pendiente restringidas a una ventana de 3x3 a multiples escalas utilizando solo las celdas en los bordes de la ventana focal.
argumentos:
r: DTM como un SpatRaster o RasterLayer en un sistema de coordenadas proyectadas donde las unidades del mapa coinciden con las unidades de elevación y profundidad.
w: Un vector de longitud 2 que especifica las dimensiones de la ventana rectangular a utilizar, donde el primer número es el número de filas y el segundo número es el número de columnas. El tamaño de la ventana debe ser un número impar.
unit: grados o radianes
method: “tower”, “queen” (por defecto) o “boundary”. El método indica qué celdas utilizar en los cálculos. “torre” utiliza solo las 4 celdas de los bordes directamente arriba, abajo, izquierda y derecha; “reina” agrega cuatro celdas adicionales en las esquinas; “frontera” utiliza todas las celdas de los bordes
detalles: La pendiente se calcula atan(sqrt(dz.dx^2 + dz.dy^2)) y el aspecto se calcula como (-pi/2)-atan_2(dz.dy, dz.dx) y luego se restringe de 0 a 2 pi/0 a 360 grados. dz.dx es la diferencia entre la media ponderada del lado derecho de la ventana focal y la media ponderada del lado izquierdo de la ventana focal dividida por la distancia x de la ventana focal en unidades de mapa. dz.dy es la diferencia entre la media ponderada de la parte superior de la ventana focal y la media ponderada de la parte inferior de la ventana focal dividida por la distancia y de la ventana focal en unidades de mapa. Las celdas utilizadas en estos cálculos dependen del metodo elegido. Para los métodos “queen” y “limit”, las celdas de esquina tienen la mitad del peso de todas las demás celdas utilizadas en los cálculos.
#por medio de la funcion SlpAsp calculamos los atributos geomorfologicos pendiente y orientación (slope, aspect) a partir del DEM.
#La pendiente se calcula en grados
#La orientación se calcula en grados azimutales, es decir que el norte esta los 0° y el este a los 90°, y asi sucesivamente.
(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 : 1785, 1034, 2 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 0.0000
## max values : 61.32469, 359.9998ahora dividimos el objeto slp_asp en dos partes, la pendiente y la orientación.
# extraemos la primera capa del objeto Slp_asp, que corresponde a la pendiente en grados
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 61.32469Podemos observar que la pendiente del departamento varia desde los 0 hasta los 62° aproximadamente, lo cual se puede explicar debido a las grandes llanuras que presenta el departamento en la cuenca del rio Magdalena y la gran pendiente por la sierra nevada de santa marta.
#Este codigo genera un historigrama que muestra la distribucion de la pendiente del raster utilizando la funcion hist() de la libreria "terra".
terra::hist(slope,
main = "pendiente del Magdalena",
xlab = "Pendiente (en grados)")## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
Con esta grafica podemos evidenciar que los valores de la pendiente se concentran entre los 0° y 5° con una frecuencia superior a 300000 celdas, indicando que gran parte del departamento del Magdalena es plano con zonas ligeramente inclinados, probablemente en zonas costeras, la cuenca del rio magdalena y las zonas agricolas del bajo Magdalena.
Tambien podemos evidenciar que hay un menor numero de celdas entre los 100° y 30°, y aún menos hacia los 40° a 60°, por lo que estas zonas corresponden a zonas de relieve montañoso o escarpado, como la sierra nevada de santa marta o las zonas altas del sur.
#Extraemos la segunda capa del objeto slp_asp, que corresponde a la orientación de las pendientes.
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.0000
## max value : 359.9998A partir de esta información podemos tomar decisiones importantes respecto a la exposición solar que puede recibir un cultivo y ubicar las laderas con menor insolacion, pero para poder analizar esto, es necesario realizar un diagrama de barras que permita observar el comportamiento de los datos.
#Con este codigo se genera un historigrama que indica la orientación de las pendientes en el departamento del Magdalena.
terra::hist(aspect,
main = "Orientacion del magdalena",
xlab = "Orientacion (en grados)")## Warning: [hist] a sample of 54% of the cells was used (of which 45% was NA)
Mediante este grafico se puede analizar la orientación de las pendientes, donde cada grado significa lo siguiente:
0° = Norte
90° = Este
180° = Sur
270° = Oeste
360° ≈ Norte otra vez
Por lo que se puede obervar como hay mas celdas orientadas entre los 270° y 330°, es decir hacia el oeste y noroeste, por lo que estas laderas pueden recibir más radiación solar en la tarde. siendo util para cultivos que toleran mas calor. Por otro lado tambien se puede observar que hay menor frecuencia de celdas entre los 60° y los 120°, es decir menos laderas orientadas al este, por lo que son zonas conmenor insolacion, siendo utiles para cultivos que requieran más humedad.
Ahora convertiremos la pendiente en grados a pendiente en porcentaje.
#Con esta función podemos convertir los grados a porcentaje para la pendiente.
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1785, 1034, 1 (nrow, ncol, nlyr)
## resolution : 149.7512, 149.7512 (x, y)
## extent : 4786059, 4940902, 2545552, 2812858 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 182.8408La pendiente minima es del 0% por lo que hay zonas completamente planas, ttipico en muchas partes del magdalena, mientras que la pendiente maxima es del 182.84% lo que indica terrenos muy inclinados que probablemente esten inclinados en la sierra nevada de santa marta. Esta información es util en la agricultura para el ordenamiento del territorio, infraestructura agraria y ordenamiento territorial.
Pero para poder analizar mejor la información es neceario un grafico que permita visualizar los datos de manera mas precisa.
terra::hist(slope_perc,
main = "Pendiente en el Magdalena",
xlab = "pendiente (en porcentaje)")## Warning: [hist] a sample of 54% of the cells was used (of which 44% was NA)
La mayoria de los datos de pendiente se encuentran concentrados entre el 0% y el 10%, indicando que el departamento del magdalena es mayormente plano o con leves inclinaciones. A medida que la pendiente aumenta, la frecuencia disminuye drasticamente.
Una pequeña fracción del departamento refleja pendientes superioes al 50% que corresponden a la sierra nevada o regiones montañosas aisladas.
Con esta información se puede concluir que el departamento del Magdalena presenta una alta aptitud agricola gracias a la alta frecuencia de pendientes bajsa reflejando el potencial que tiene el departamento para la producción agricola. Tambien se puede concluir que las areas con una pendiente elevada deben ser analizadas para conservación, protección contra la erosión o cultivos especiales adaptados a terrenos inclinados.
7) Cálculo de estadísticas zonales (20%)
Las estadisticas zonales permiten resumir datos raster como pendiente, altitud o humedad dentro de areas definidas por zonas que pueden ser poligonos administrativos, unidades de cultivo, cuencas, etc.
Suelen usarse en agricultura de precisión para caracterizar zonas del campo segun pendiente o inclinación para ajustar la fertilización o el riego.
En este caso vamos a calcular el promedio de los valores de pendiente por municipio, por lo que vamos a examinar la clasificación de pendiente segun el IGAC para la elaboración del mapa.
con esta clasificación procederemos a reclasificar el raster de la pendiente.
#Vamos a agrupar los valores del raster slope_perc en clases definidas por rangos.
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
75, 185, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)Esta reclasificación es util para realizar analisis agricolas de eroción, planificación, etc.
A continuación podemos calcular la estadisticas zonales usando exactextractr.
## Calculamos la pendiente media dentro de cada poligono del objeto munic apartir del raster slope_perc y lo guardamos en la columna mean_slope.
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 3% | |===== | 7% | |======= | 10% | |========= | 13% | |============ | 17% | |============== | 20% | |================ | 23% | |=================== | 27% | |===================== | 30% | |======================= | 33% | |========================== | 37% | |============================ | 40% | |============================== | 43% | |================================= | 47% | |=================================== | 50% | |===================================== | 53% | |======================================== | 57% | |========================================== | 60% | |============================================ | 63% | |=============================================== | 67% | |================================================= | 70% | |=================================================== | 73% | |====================================================== | 77% | |======================================================== | 80% | |========================================================== | 83% | |============================================================= | 87% | |=============================================================== | 90% | |================================================================= | 93% | |==================================================================== | 97% | |======================================================================| 100%## [1] 34.3130493 0.9711066 32.4357758 1.3878907 0.9487980 1.5636197
## [7] 34.1772461 1.8798674 0.9540940 0.7222221 0.5690166 21.9244385
## [13] 0.9260811 1.3292797 2.5145619 1.1275929 1.0748733 1.7687278
## [19] 0.4949985 0.5313346 1.6415391 0.4778264 1.2094039 0.9828508
## [25] 1.0238775 1.5891575 0.5086798 2.3828094 1.4427539 2.6509492Una ve calculadas las estadisticas zonales podemos observar que da como resultado la pendiente media por poligono, indicando el porcentaje de la pendientre promedio de un poligono, demostrando de esta manera una gran variablidad en la pendiente media entre zonas.
hist(munic$mean_slope,
main = "la pendiente media del magdalena",
xlab = "Pendiente (en porcentaje)")
Por medio de este grafico se puede evidenciar como es que la mayoria de los municipios tienen pedientes medias bajas, que se encuentran entre el 0% y el 5%, que indica qque predominan zonas planas o con un relieve suave, tambien hay presencia de municipios con pendientes medias mayores a 10% y unos pocos con valores entre el 20% y el 35%, Por lo cual se puede concluir que el departamento del magdalena es predominantemente plano.
Una vez que hemos calculado la pendiente en porcentaje, la vamos a clasificar como una categoria.
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode')) ## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 3% | |===== | 7% | |======= | 10% | |========= | 13% | |============ | 17% | |============== | 20% | |================ | 23% | |=================== | 27% | |===================== | 30% | |======================= | 33% | |========================== | 37% | |============================ | 40% | |============================== | 43% | |================================= | 47% | |=================================== | 50% | |===================================== | 53% | |======================================== | 57% | |========================================== | 60% | |============================================ | 63% | |=============================================== | 67% | |================================================= | 70% | |=================================================== | 73% | |====================================================== | 77% | |======================================================== | 80% | |========================================================== | 83% | |============================================================= | 87% | |=============================================================== | 90% | |================================================================= | 93% | |==================================================================== | 97% | |======================================================================| 100%## [1] 5 1 5 1 1 1 5 1 1 1 1 5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1## Calculamos la clase modal de la pendiente como una categoria y la guardamos en una nueva variable munic$class.En el resultado de este comando se puede evidenciar como es que la mayoria de los municipios tienen una pendiente de calse 1, lo que indica que la pendiente mas común esta entre 0% y 3%.
Ahora procedemos a realizar un grafico para observar mejor el comportamiento de los datos.
hist(munic$class,
main = "la pendiente reclasificada del Magdalena",
xlab = "Pendiente (como una categoria)")
En este grafico se puede corroborar que mas de 25 municipios del departamento del magdalena tienen pendientes predominantemente entre 0% y 3%, confirmando que la mayoria del departamento es plano, mientras que la clase 5 aparece en 2 o 3 municipios.
Ahora procedemos a reproyectar el raaster rc desde el sistema de coordenadas MAGNA-SIRGAS hacia el sistema de coordenadas geograficas WGS 84 que utiliza grados decimales de longitud y latitud.
## reproyectamos el raster a corrdenadas geograficas por medio del comando "rc".
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1785, 1048, 1 (nrow, ncol, nlyr)
## resolution : 0.001357862, 0.001357862 (x, y)
## extent : -74.96115, -73.53811, 8.931431, 11.35521 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7Podemos evidenciar que el raster contiene 1785 filas por 1048 columnas, lo que indica que tiene un total de 1.87 millones de celdas, en donde cada pixel tiene una resolución de 0.001357862 grados, lo que es alrededor de 151 metros en el ecuador, tambien se pueden ver los valores minimo y maximo de la pendiente reclasificada, que va desde 0 hasta 167.75 grados.
Ahora transfromamos el raster al sistema de coordenadas EPSG:4326 para poder realizar el mapa mas adelante.
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1785, 1048, 1 (nrow, ncol, nlyr)
## resolution : 0.001357862, 0.001357862 (x, y)
## extent : -74.96115, -73.53811, 8.931431, 11.35521 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 167.75318) Trazando la pendiente del terreno (20%)
Para realizar el mapa interactivo de la pendiente en el departamento del magdalena, lo primero que debemos hacer es crear una paleta de colores continua para poder visualizar el raster de pendientes con un degradado que vaya desde verde claro a rojo oscuro.
##Creamos una función de color por medio del comando ColorNumeric() basada en un gradiente númerico.
palredgreen <- colorNumeric(c("darkseagreen3","yellow", "orange2", "brown3", "darkred"), values(slope.geo),
na.color = "transparent")Por medio de esta gradiente de colores se puede representar de manera clara las zonas planas en color verde, suavemente inclinadas en naranja, empinadas o motañosas en rojo, siendo ideal para el análisis agrícola, la planificación del uso del suelo o estudios de erosión.
Procedemos a crear un mapa interactivo que combine el vector de los municipios con información sobre el nombre y la clase de pendiente, el raster de la pendiente y la paleta continua de colores en donde se representen los cambios en la pendiente.
## creamos el mapa a partir del codigo leaflet(munic) como capa base
leaflet(munic) %>% addTiles() %>% setView(-74.3, 10.3, 7.9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Clase de pendiente: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "pendiente del terreno en Magdalena (%)")## 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'