Atributos geomorfométricos del departamento de Arauca
Daniel Castiblanco
23/07/2025
Atributos de terreno geomorfológicos en R
- Introducción 🌍
- Configuración 🛠️
- Introducción a la librería
multiscaleDTM 🧪
- Leyendo datos de entrada
🧠
- Transformar coordenadas
(10%) 🔧
- Cálculo de atributos del terreno
(20%) 🌿
- Cálculo de estadísticas
zonales (20%) 📊
- Trazando la pendiente del terreno (20%)
🧭
- Bibliografía 📚
🌍1) Introducción
Este cuaderno tiene como objetivo aplicar herramientas de R para calcular variables geomorfométricas a partir de un modelo de elevación en formato ráster. En particular, se empleará la librería MultiscaleDTM para calcular dos propiedades clave del terreno en el departamento del Magdalena: la pendiente y la orientación.
Ambas variables son esenciales en agricultura:
La pendiente permite identificar cómo debe fluir el agua para evitar acumulaciones y mejorar el drenaje.
La orientación del terreno facilita la planeación agrícola, la ubicación de cultivos, el diseño de drenajes, y permite interpretar condiciones de iluminación, temperatura y viento, aspectos claves en la productividad.
🛠️ 2)Configuracion
Preparación del entorno Antes de comenzar, es necesario cargar varias librerías de R que facilitarán el procesamiento espacial:
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 TRUElibrary(leaflet)
library(terra)## terra 1.8.54library(MultiscaleDTM)
library(exactextractr)(elevatr ) permite descargar datos de elevación en forma de objetos sf o raster desde servicios en línea. Aunque está migrando hacia terra, sigue ofreciendo funciones compatibles con raster.
(sf) facilita el manejo de datos vectoriales basados en geometrías simples.
(leaflet) se utiliza para crear mapas web interactivos, integrando datos geoespaciales con visualizaciones dinámicas.
(terra) trabaja con datos espaciales raster y vectoriales, y proporciona herramientas para realizar análisis espaciales complejos.
(exactextractr) permite calcular estadísticas de valores raster dentro de polígonos, incluso si las celdas están parcialmente cubiertas.
🧩 3) Introducción a la librería MultiscaleDTM
MultiscaleDTM es un paquete especializado en el análisis multiescalar de modelos digitales del terreno. Su principal utilidad es extraer atributos topográficos como pendiente, curvatura, rugosidad, orientación, TPI (índice de posición topográfica), entre otros, usando ventanas móviles de distintos tamaños. Esta capacidad multiescalar es muy útil para estudios de geomorfología, hidrología y agricultura de precisión.
La librería trabaja con objetos tipo SpatRaster del paquete terra, facilitando su integración con flujos de trabajo espaciales en R.
knitr::include_graphics("C:\\Users\\Pc\\Pictures\\Screenshots\\Captura_de_pantalla.png")
##Funciones principales
SlpAsp: Esta función se utiliza para calcular la pendiente y la orientación del terreno, considerando un subconjunto de celdas dentro de una ventana focal. Dependiendo del método seleccionado, se consideran diferentes celdas:
1). En el método de “torre”, se usan las cuatro celdas directamente adyacentes (superior, inferior, izquierda y derecha).
2). En el método de “reina”, además de las anteriores, se incluyen las cuatro celdas en las esquinas.
3). Con el método de “límite”, se consideran todas las celdas que rodean el borde de la ventana focal.
knitr::include_graphics("C:\\Users\\Pc\\Downloads\\descarga.png")
DirSlp: Esta función estima la pendiente a múltiples escalas en una dirección previamente definida, lo que permite analizar variaciones direccionales del relieve dentro de una superficie.
Qfit: Esta herramienta calcula la pendiente, orientación, curvatura y otras propiedades morfométricas del terreno ajustando una superficie cuadrática sobre la ventana focal. El ajuste se realiza mediante mínimos cuadrados ordinarios a partir de la ecuación:
knitr::include_graphics("C:\\Users\\Pc\\Pictures\\Screenshots\\Captura de pantalla 2025-07-21 185309.png")
donde los coeficientes 𝑎 a a 𝑓 f son parámetros obtenidos por regresión, 𝑍 Z representa la elevación o profundidad, y 𝑋 X e 𝑌 Y corresponden a las coordenadas relativas (este/oeste y norte/sur, respectivamente) dentro de la ventana focal respecto a la celda central.
¿Deseas que integre estos fragmentos en un texto más largo, como parte de un informe o presentación?
knitr::include_graphics("C:\\Users\\Pc\\Downloads\\descarga (1).png")
Pfit calcula la pendiente y el aspecto a partir de un ajuste plano en la ventana focal, usando mínimos cuadrados ordinarios. Es menos exigente computacionalmente que un ajuste cuadrático, pero ofrece resultados similares.
Posición relativa indica si una celda es un punto alto o bajo en relación con su entorno. Se obtiene restando el valor de la celda focal a una elevación de referencia (media, mediana, mínimo, etc.). Puede expresarse en unidades del DTM o estandarizarse según la topografía local.
Aplicación en mapas agrícolas:
Esta librería permite calcular atributos topográficos como pendiente, curvatura y entropía desde modelos digitales de elevación. Estos son clave para evaluar drenaje, riesgo de erosión y variabilidad del suelo, facilitando la zonificación agrícola, el manejo del riego y fertilización, y la toma de decisiones en agricultura de precisión. Los resultados pueden combinarse con capas vectoriales para análisis espaciales más detallados.
knitr::include_graphics("C:\\Users\\Pc\\Downloads\\descarga (2).png")
4) Lectura de datos de entrada 🧠
El primer paso consiste en cargar al entorno de trabajo la capa ráster correspondiente al departamento de Arauca, la cual debe haber sido previamente descargada. Esta capa, que representa el modelo digital de elevación (DEM) de la zona, servirá como base para la extracción de atributos geomorfológicos que permitirán caracterizar el relieve y apoyar el análisis territorial en el contexto agrícola del departamento.
# 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\\Pc\\Desktop\\COGUA\\COGUA\\P1\\ALTURA\\ARAUCA.tif")#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")# Instalar librerías si es necesario:
# install.packages(c("sf", "dplyr", "geodata"))
library(sf)
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(geodata) # cómodo para bajar GADM## Warning: package 'geodata' was built under R version 4.5.1# === OPCIÓN A: bajar GADM en una sola línea ===
gadm_lvl2 <- geodata::gadm("COL", level = 2, path = tempdir()) # SpatVector
muni_col <- st_as_sf(gadm_lvl2) # a 'sf'
# Filtrar solo Arauca (NAME_1 == "Arauca")
muni_ara <- muni_col %>% filter(NAME_1 == "Arauca")
# Guardar como GeoPackage (una sola capa llamada "municipios_arauca")
st_write(muni_ara,
"arauca_municipios.gpkg",
layer = "municipios_arauca",
driver = "GPKG",
delete_dsn = TRUE) # sobreescribe si ya existe## Deleting source `arauca_municipios.gpkg' using driver `GPKG'
## Writing layer `municipios_arauca' to data source
## `arauca_municipios.gpkg' using driver `GPKG'
## Writing 7 features with 13 fields and geometry type Polygon.munic <- sf::st_read("arauca_municipios.gpkg")## Reading layer `municipios_arauca' from data source
## `C:\Users\Pc\Desktop\GB2\TRABAJO INFORME #2\arauca_municipios.gpkg'
## using driver `GPKG'
## Simple feature collection with 7 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36546 ymin: 6.033162 xmax: -69.42688 ymax: 7.104162
## Geodetic CRS: WGS 84plot(munic["NAME_2"]) # si el campo con nombres municipales es NAME_2
normalizePath("arauca_municipios.gpkg")## [1] "C:\\Users\\Pc\\Desktop\\GB2\\TRABAJO INFORME #2\\arauca_municipios.gpkg"(munic <- sf::st_read("C:\\Users\\Pc\\Desktop\\GB2\\TRABAJO INFORME #2\\arauca_municipios.gpkg"))## Reading layer `municipios_arauca' from data source
## `C:\Users\Pc\Desktop\GB2\TRABAJO INFORME #2\arauca_municipios.gpkg'
## using driver `GPKG'
## Simple feature collection with 7 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36546 ymin: 6.033162 xmax: -69.42688 ymax: 7.104162
## Geodetic CRS: WGS 84## Simple feature collection with 7 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -72.36546 ymin: 6.033162 xmax: -69.42688 ymax: 7.104162
## Geodetic CRS: WGS 84
## GID_2 GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1 NAME_2 VARNAME_2
## 1 COL.3.1_2 COL Colombia COL.3_2 Arauca <NA> Arauca <NA>
## 2 COL.3.2_2 COL Colombia COL.3_2 Arauca <NA> Arauquita <NA>
## 3 COL.3.3_2 COL Colombia COL.3_2 Arauca <NA> Cravo Norte <NA>
## 4 COL.3.4_2 COL Colombia COL.3_2 Arauca <NA> Fortul <NA>
## 5 COL.3.5_2 COL Colombia COL.3_2 Arauca <NA> Puerto Rondón <NA>
## 6 COL.3.6_2 COL Colombia COL.3_2 Arauca <NA> Saravena <NA>
## 7 COL.3.7_2 COL Colombia COL.3_2 Arauca <NA> Tame <NA>
## NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2 geom
## 1 <NA> Municipio Municipality <NA> CO.AR.AC POLYGON ((-70.65508 6.48394...
## 2 <NA> Municipio Municipality <NA> CO.AR.AQ POLYGON ((-71.09048 6.98354...
## 3 <NA> Municipio Municipality <NA> CO.AR.CN POLYGON ((-70.65508 6.48394...
## 4 <NA> Municipio Municipality <NA> CO.AR.FO POLYGON ((-71.69981 6.86233...
## 5 <NA> Municipio Municipality <NA> CO.AR.PR POLYGON ((-70.88526 6.62282...
## 6 <NA> Municipio Municipality <NA> CO.AR.SA POLYGON ((-71.69448 7.03584...
## 7 <NA> Municipio Municipality <NA> CO.AR.TA POLYGON ((-71.47402 6.67631...#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)#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 : 65, 176, 1 (nrow, ncol, nlyr)
## resolution : 1843.165, 1843.165 (x, y)
## extent : 5070831, 5395228, 2224257, 2344063 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : ARAUCA
## min value : 84.46235
## max value : 4371.05371####Convertimos el archivo que contiene los polígonos de los municipios al sistema de referencia espacial MAGNA-SIRGAS. Es fundamental que tanto el raster como el vector compartan el mismo sistema de coordenadas, ya que esto permite que interactúen adecuadamente. Por ejemplo, al recortar el raster usando el shapefile o al extraer valores por zonas, se requiere esa compatibilidad.
#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"))## Simple feature collection with 7 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 5070157 ymin: 2225666 xmax: 5395455 ymax: 2343637
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## GID_2 GID_0 COUNTRY GID_1 NAME_1 NL_NAME_1 NAME_2 VARNAME_2
## 1 COL.3.1_2 COL Colombia COL.3_2 Arauca <NA> Arauca <NA>
## 2 COL.3.2_2 COL Colombia COL.3_2 Arauca <NA> Arauquita <NA>
## 3 COL.3.3_2 COL Colombia COL.3_2 Arauca <NA> Cravo Norte <NA>
## 4 COL.3.4_2 COL Colombia COL.3_2 Arauca <NA> Fortul <NA>
## 5 COL.3.5_2 COL Colombia COL.3_2 Arauca <NA> Puerto Rondón <NA>
## 6 COL.3.6_2 COL Colombia COL.3_2 Arauca <NA> Saravena <NA>
## 7 COL.3.7_2 COL Colombia COL.3_2 Arauca <NA> Tame <NA>
## NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2 geom
## 1 <NA> Municipio Municipality <NA> CO.AR.AC POLYGON ((5259240 2275063, ...
## 2 <NA> Municipio Municipality <NA> CO.AR.AQ POLYGON ((5210870 2330097, ...
## 3 <NA> Municipio Municipality <NA> CO.AR.CN POLYGON ((5259240 2275063, ...
## 4 <NA> Municipio Municipality <NA> CO.AR.FO POLYGON ((5143603 2316471, ...
## 5 <NA> Municipio Municipality <NA> CO.AR.PR POLYGON ((5233716 2290308, ...
## 6 <NA> Municipio Municipality <NA> CO.AR.SA POLYGON ((5144139 2335651, ...
## 7 <NA> Municipio Municipality <NA> CO.AR.TA POLYGON ((5168611 2295982, ...6) Cálculo de atributos del terreno 🌱
Se estiman la pendiente y el aspecto del terreno a múltiples escalas usando la función SlpAsp. Esta herramienta aplica fórmulas tradicionales para el cálculo de pendientes dentro de una ventana de 3x3, considerando únicamente las celdas ubicadas en los bordes de dicha ventana, y permite operar a diferentes escalas espaciales.
Parámetros principales:
r: Modelo Digital del Terreno (DTM) en formato SpatRaster o RasterLayer, el cual debe estar en un sistema de coordenadas proyectadas donde las unidades del mapa coincidan con las de elevación.
w: Vector de dos elementos que define el tamaño de la ventana de análisis. El primer valor representa el número de filas y el segundo, el de columnas. Ambos deben ser valores impares.
unit: Unidad para expresar los resultados, ya sea en grados o radianes.
method: Método para definir qué celdas se incluirán en los cálculos. Puede ser:
“tower”: solo usa las celdas de los bordes vertical y horizontal (arriba, abajo, izquierda y derecha).
“queen” (valor por defecto): incluye también las esquinas.
“boundary”: emplea todas las celdas de los bordes.
Cálculo detallado: La pendiente se obtiene mediante atan(sqrt(dz.dx^2 + dz.dy^2)), y el aspecto se determina con (-π/2) - atan2(dz.dy, dz.dx), restringiéndolo a un rango entre 0 y 2π (o de 0 a 360° si se usan grados).
dz.dx representa la diferencia entre los promedios ponderados de las celdas del lado derecho e izquierdo, dividida por la distancia horizontal.
dz.dy corresponde a la diferencia entre los promedios ponderados de las celdas superiores e inferiores, dividida por la distancia vertical.
El conjunto de celdas consideradas en los cálculos depende del método elegido. En los métodos “queen” y “boundary”, las celdas ubicadas en las esquinas reciben un peso equivalente a la mitad del asignado a las celdas restantes.
#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 : 65, 176, 2 (nrow, ncol, nlyr)
## resolution : 1843.165, 1843.165 (x, y)
## extent : 5070831, 5395228, 2224257, 2344063 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 6.831753e-04, 0.1579205
## max values : 1.508657e+01, 359.8435459ahora 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 : 65, 176, 1 (nrow, ncol, nlyr)
## resolution : 1843.165, 1843.165 (x, y)
## extent : 5070831, 5395228, 2224257, 2344063 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 6.831753e-04
## max value : 1.508657e+01Podemos observar que la pendiente del departamento de Arauca varía aproximadamente entre 0 y 15°, lo cual se explica por la presencia de extensas llanuras en la región de la Orinoquía, caracterizadas por un relieve predominantemente plano. Las leves pendientes se concentran principalmente en las zonas de transición hacia la cordillera oriental, al suroeste del departamento, donde comienzan a aparecer formaciones más elevadas.
#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 de Arauca",
xlab = "Pendiente (en grados)")
Con esta gráfica podemos evidenciar que los valores de la pendiente en
Arauca se concentran mayoritariamente entre 0° y 5°, con una frecuencia
muy alta que supera las 5.000 celdas, lo que indica que gran parte del
territorio tiene un relieve plano o suavemente inclinado. Esto es
coherente con la geografía del departamento, que forma parte de la
región de la Orinoquía colombiana, caracterizada por extensas llanuras
de inundación y sabanas.
También se observa que a medida que la pendiente aumenta, la frecuencia disminuye drásticamente. Las pendientes entre 5° y 15° están presentes en una proporción mucho menor, lo que sugiere que estas zonas corresponden a áreas de transición topográfica hacia la cordillera Oriental, localizadas principalmente en el suroeste del departamento.
#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 : 65, 176, 1 (nrow, ncol, nlyr)
## resolution : 1843.165, 1843.165 (x, y)
## extent : 5070831, 5395228, 2224257, 2344063 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.1579205
## max value : 359.8435459A partir de esta información, es posible tomar decisiones clave sobre la exposición solar que pueden recibir los cultivos, identificando las laderas con menor insolación y, por ende, con condiciones más favorables o desfavorables según el tipo de cultivo. Sin embargo, para poder analizar adecuadamente este factor, es necesario elaborar un diagrama de barras que permita visualizar el comportamiento del aspecto del terreno, es decir, la orientación de las pendientes respecto al sol, lo cual facilita una interpretación más precisa en términos agroecológicos y de planificación territorial.
#Con el siguiente código se genera un histograma que muestra la orientación de las pendientes en el departamento de Arauca:
terra::hist(aspect,
main = "Orientación del terreno en Arauca",
xlab = "Orientación (en grados)")
Mediante este gráfico se puede analizar la orientación de las pendientes
en Arauca. Cada valor en grados representa una dirección cardinal
aproximada:
0° = Norte
90° = Este
180° = Sur
270° = Oeste
360° ≈ Norte nuevamente
A partir del histograma, se observa que hay una mayor concentración de celdas orientadas entre los 270° y 330°, es decir, hacia el oeste y noroeste. Estas laderas tienden a recibir mayor radiación solar en horas de la tarde, lo que podría ser favorable para cultivos que toleran altas temperaturas y requieren más luz solar.
Por otro lado, se nota una menor frecuencia de celdas entre los 60° y 120°, que corresponden a laderas orientadas al este. Estas zonas reciben el sol principalmente en la mañana y pueden mantener mayores niveles de humedad, lo cual es beneficioso para cultivos que prefieren condiciones más frescas y húmedas.
#Con esta función podemos convertir los grados a porcentaje para la pendiente.
(slope_perc = tan(slope*(pi/180))*100) ## class : SpatRaster
## size : 65, 176, 1 (nrow, ncol, nlyr)
## resolution : 1843.165, 1843.165 (x, y)
## extent : 5070831, 5395228, 2224257, 2344063 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.001192366
## max value : 26.956926780La pendiente mínima del terreno en el departamento de Arauca es cercana al 0%, lo que indica la presencia de zonas completamente planas, típicas de las llanuras de la Orinoquía. La pendiente máxima alcanza aproximadamente 27°, lo que equivale a un 182.84% al convertirla a porcentaje mediante la fórmula: pendiente (%) = tan(°) × 100. Este valor sugiere la existencia de terrenos muy inclinados, localizados probablemente en las zonas de transición hacia la cordillera Oriental, en el suroeste del departamento. Esta información es fundamental para la planificación agrícola, el diseño de infraestructura rural y el ordenamiento territorial. No obstante, para interpretar de forma más precisa estos datos, es recomendable generar un gráfico que permita visualizar la distribución de la pendiente en porcentaje de manera clara y detallada.
terra::hist(slope_perc,
main = "Pendiente en el Arauca",
xlab = "pendiente (en porcentaje)")
La mayoría de los valores de pendiente en el departamento de Arauca se concentran entre el 0% y el 10%, con una frecuencia muy elevada, lo que confirma que la región es mayoritariamente plana o con leves inclinaciones. Esta característica es coherente con la geografía de Arauca, ubicada en la Orinoquía colombiana, una zona compuesta principalmente por llanuras aluviales y sabanas. A medida que aumenta el valor de la pendiente, la frecuencia disminuye drásticamente, indicando que las zonas con pendientes más pronunciadas son escasas y se concentran en áreas de transición hacia la cordillera Oriental, en el suroeste del departamento. Aunque se identifican pendientes que alcanzan hasta el 25%, estas representan una fracción muy pequeña del territorio. Con base en esta información, se puede concluir que Arauca posee una alta aptitud agrícola, gracias a la predominancia de terrenos con baja pendiente, lo que favorece actividades como el cultivo extensivo y la mecanización agrícola. No obstante, las zonas con pendientes más altas deben ser consideradas estratégicamente para procesos de conservación, protección contra la erosión o para la implementación de cultivos adaptados a condiciones inclinadas.
7) Cálculo de estadísticas zonales
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.
knitr::include_graphics("C:\\Users\\Pc\\Downloads\\descarga (3).png")
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.
knitr::include_graphics("C:\\Users\\Pc\\Downloads\\descarga (4).png")
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.
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========== | 14% | |==================== | 29% | |============================== | 43% | |======================================== | 57% | |================================================== | 71% | |============================================================ | 86% | |======================================================================| 100%## [1] 0.04428873 0.10874008 0.04845228 7.07922697 0.06909953 2.44810247 5.51830435hist(munic$mean_slope,
main = "la pendiente media Arauca",
xlab = "Pendiente (en porcentaje)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |========== | 14% | |==================== | 29% | |============================== | 43% | |======================================== | 57% | |================================================== | 71% | |============================================================ | 86% | |======================================================================| 100%## [1] 1 1 1 1 1 1 1hist(munic$class,
main = "la pendiente reclasificada Arauca",
xlab = "Pendiente (como una categoria)")
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 66, 176, 1 (nrow, ncol, nlyr)
## resolution : 0.01666668, 0.01666668 (x, y)
## extent : -72.3597, -69.42637, 6.013347, 7.113348 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 5(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 66, 176, 1 (nrow, ncol, nlyr)
## resolution : 0.01666668, 0.01666668 (x, y)
## extent : -72.3597, -69.42637, 6.013347, 7.113348 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.003513035
## max value : 26.5141601568) Trazando la pendiente del terreno
palredgreen <- colorNumeric(c("darkseagreen3","yellow", "orange2", "brown3", "darkred"), values(slope.geo),
na.color = "transparent")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 (%)")9) Bibliografía
Beers, T., Dress, P. E., & Wensel, L. C. (1966). Aspect transformation in site productivity research. Journal of Forestry, 64(10), 691–692. (Fundamento sobre el uso de orientación en grados y su interpretación ecológica)
Dunnington, D. (2023). exactextractr: Fast Extraction from Raster Datasets Using Polygons (versión 0.9.2) [Paquete R]. https://github.com/isciences/exactextractr (Usado para calcular estadísticas zonales precisas)
Hijmans, R. J. (2023). terra: Spatial Data Analysis (versión 1.7-65) [Paquete R]. https://rspatial.org/terra/
Instituto Geográfico Agustín Codazzi (IGAC). (2020). Implementación del nuevo sistema de coordenadas planas con origen nacional EPSG:9377 (MAGNA-SIRGAS 2018). https://geoportal.igac.gov.co/
Lizarazo, I. (2024). Geomorphometric terrain attributes using R. RPubs. https://rpubs.com/ials2un/geomorphometric (Tutorial base del que se usaron los fragmentos de código)
Pebesma, E. (2018). Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal, 10(1), 439–446. https://doi.org/10.32614/RJ-2018-009 (Base teórica del paquete sf, usado para datos vectoriales)
R Core Team. (2023). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/