Atributos Geomorfometricos en R
1) INTRODUCCION
En este cuaderno vamos a trabajar en como calcular las variables geomorfometricas del departamento del Bolivar en Colombia , se obtendran estos datos mediante el uso de modelos digital de elevación, en esta ocasión se usara el paquete MultiscaleDTM por lo tanto nuestro primer paso es limpiar la memoria que se hallaba en R
rm(list=ls())1.1)Instalación de paquetes
Se instalaran los paquetes requeridos para poder trabajar en este cuaderno esto se realiza una sola vez por lo tanto se pueden quitar los numerales y volveros a correr o correr directamente en la consola de R
#paquetes = c("MultiscaleDTM", "exactextractr")
#install.packages(paquetes)1.2) Cargar librerias
En caso de haber conflicto en las librerias revisar que se halla instalado correctamente los paquetes
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)2) Introducción a MultiscaleDTM
El paquete MultiscleDTM contiene 31 funciones en R para calcular diversas cosas siendo la pendiente una de ellas si quieres saber mas información sobre este paquete y sus funciones has click aca a en este link
3) Carga del archivo a trabajar
Para este cuardeno necesitaremos un raster DEM obtenido de un cuaderno anterior este va ser el primer paso para poder determinar la pendiente en nuestro departamento
(dem = terra::rast("./DEPARTAMENTOS/elev_bolivar_z10.tif"))## class : SpatRaster
## size : 6121, 4130, 1 (nrow, ncol, nlyr)
## resolution : 0.0006810161, 0.0006810161 (x, y)
## extent : -76.28906, -73.47647, 6.664807, 10.83331 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_bolivar_z10.tif
## name : elev_bolivar_z10
## min value : -2790
## max value : 5577verifica que tu ruta sea correcta para que se pueda leer el archivo
3.1) Reduce el tamaño del DEM
El tamaño del raster puede que sobrecargue la memoria de R y del equipo por lo tanto se procede realizar una disminución en su tamaño
dem2 = terra::aggregate(dem,2, "mean")# Se realiza terra:: para indicar que el codigo a usar es de la libreria terra ## |---------|---------|---------|---------|========================================= 3.2) Cargar y transformación departamento
Se carga el departamento trabajar este siendo un gkpg.
(munic <- sf::st_read("./DEPARTAMENTOS/bolivarcortadocorrecto.gpkg"))# Se realiza sf::para indicar que el codigo a usar es de la libreria terra ## Reading layer `bolivar' from data source
## `C:\Users\juanc\OneDrive\Documentos\GB_CristianEljadue\GB_CristianEljadue\GB_CristianEljadue\P9\DEPARTAMENTOS\Bolivarcortadocorrecto.gpkg'
## using driver `GPKG'
## Simple feature collection with 46 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.19063 ymin: 6.99916 xmax: -73.74578 ymax: 10.80147
## Geodetic CRS: MAGNA-SIRGASPosterior a cargar el departamento procedemos a unirlo con la variable munic
dem3 = terra::crop(dem2,munic, mask=TRUE)## Warning: [crop] CRS do not matchComo se puede observar las dos variables se encuentran en diferentes sistemas de coordenadas por lo tanto es necesario que queden en un mismo sistema para trabajar en el cuaderno.
4) Transformación de coordenadas
(dem_plane = project(dem3, "EPSG:9377"))# se trabaja EPSG:9377 debido que es el sistema de coordenas establecidas en colombia desde 2020 ## class : SpatRaster
## size : 2807, 1802, 1 (nrow, ncol, nlyr)
## resolution : 150.3672, 150.3672 (x, y)
## extent : 4647514, 4918475, 2331397, 2753478 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_bolivar_z10
## min value : -486.2108
## max value : 2228.3184Como se observa ya hemos dejado a dem3 en coordenadas nacionales se realiza lo mismo para el departamento trabajado
(munic_plane = sf::st_transform(munic, "EPSG:9377"))Con las dos variables en las mismas coordendas pasamos a calcular la pendiente para nuestro departamento
5) Pendiente del departamento y aspecto
Primero vamos a calcular la pendiente y el aspecto del departamento usando la funcion SlpAsp
# w esUn vector que define las dimensiones de un rectagunlo tiene que tener dos longitudes una para columnas y otra para filas siendo el primer numero las filas y el segundo las columnas como condicion deben ser numeros impares
# Method indica cuales celdas se usan en computación pueden ser solo 4 en las direcciones de arriba,abajo,y los lados u ocho en las esquinas de esas cuatro celdas con method en default
(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 : 2807, 1802, 2 (nrow, ncol, nlyr)
## resolution : 150.3672, 150.3672 (x, y)
## extent : 4647514, 4918475, 2331397, 2753478 (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 : 55.00934, 359.9998Ahora vas a separar la pendiente y el aspecto
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 2807, 1802, 1 (nrow, ncol, nlyr)
## resolution : 150.3672, 150.3672 (x, y)
## extent : 4647514, 4918475, 2331397, 2753478 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 55.00934Graficamos con un histrograma y observamos que la pendiente en grados
terra::hist(slope,
main = "pendientes del Bolivar",
xlab = "pendiente (grados)")## Warning: [hist] a sample of 20% of the cells was used (of which 77% was NA)
Donde observamos que en el departamento del bolivar la mayoria de terrenos son planos o con una inclinación de 2 grados.
Ahora observaremos el aspecto para esto primero tenemos que que separar el
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 2807, 1802, 1 (nrow, ncol, nlyr)
## resolution : 150.3672, 150.3672 (x, y)
## extent : 4647514, 4918475, 2331397, 2753478 (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.9998Con esta separación ya podemos observar como es el aspecto en el departamento
terra::hist(aspect,
main = "Aspectos del Bolivar ",
xlab = "Aspecto (grados)")## Warning: [hist] a sample of 20% of the cells was used (of which 77% was NA)
Podemos observar que el volibar predominan las orientaciones hacia el noreste y el sureste
tambien podemos observar las pendientes en porcentaje
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 2807, 1802, 1 (nrow, ncol, nlyr)
## resolution : 150.3672, 150.3672 (x, y)
## extent : 4647514, 4918475, 2331397, 2753478 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 142.8644terra::hist(slope_perc,
main = "Pendientes del Bolivar",
xlab = "pendiente (porcentaje)")## Warning: [hist] a sample of 20% of the cells was used (of which 77% was NA)
Podemos observar que las pendientes van de 0 a 60% siendo las pendientes de 0% las que mayor estan presentes en el territorio esto concuerda con el histograma en grados hecho con anterioridad
6) Calculos estadisticos zonales
Haciendo uso de la tabla de clasificacion de pendientes reoganizaremos las pendientes del departamento y las clasificaremos.
Para esto reclasificaremos las pendientes halladas anteriormente
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)Uniremos al la variable donde esta el departamento el calculo zonal
(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% | |===== | 7% | |====== | 9% | |======== | 11% | |========= | 13% | |=========== | 15% | |============ | 17% | |============== | 20% | |=============== | 22% | |================= | 24% | |================== | 26% | |==================== | 28% | |===================== | 30% | |======================= | 33% | |======================== | 35% | |========================== | 37% | |=========================== | 39% | |============================= | 41% | |============================== | 43% | |================================ | 46% | |================================= | 48% | |=================================== | 50% | |===================================== | 52% | |====================================== | 54% | |======================================== | 57% | |========================================= | 59% | |=========================================== | 61% | |============================================ | 63% | |============================================== | 65% | |=============================================== | 67% | |================================================= | 70% | |================================================== | 72% | |==================================================== | 74% | |===================================================== | 76% | |======================================================= | 78% | |======================================================== | 80% | |========================================================== | 83% | |=========================================================== | 85% | |============================================================= | 87% | |============================================================== | 89% | |================================================================ | 91% | |================================================================= | 93% | |=================================================================== | 96% | |==================================================================== | 98% | |======================================================================| 100%## [1] 2.2869468 1.8157550 5.3432398 19.8650990 1.3807722 2.7685397
## [7] 5.2977285 3.4137778 10.1638975 0.7272428 2.0966194 3.9732344
## [13] 8.2282953 5.1401849 1.5566770 1.0228293 1.0181255 3.0084543
## [19] 0.6610879 2.0796244 20.8052635 0.5899771 12.8185844 15.0571480
## [25] 0.7127678 1.1724888 9.4742975 0.8035526 2.7265577 0.5259635
## [31] 10.6892557 4.0272660 9.6671095 8.4860201 12.6133575 3.8333232
## [37] 2.2899241 18.9046459 9.5007782 0.6974851 0.5744093 12.4652700
## [43] 3.9329317 2.3503470 5.8943992 2.5175984Graficaremos un histogramas con las pendientes medias del territorio
hist(munic$mean_slope,
main = "pendientes medias del Bolivar",
xlab = "pendiente (porcentaje)")
Uniremos la reclasificacion al la variable que tenga el departamento y lo graficaremos para observar cuales son las clasificaciones obtenidas para el departamento
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))#importante verificar las coordenas## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 2% | |=== | 4% | |===== | 7% | |====== | 9% | |======== | 11% | |========= | 13% | |=========== | 15% | |============ | 17% | |============== | 20% | |=============== | 22% | |================= | 24% | |================== | 26% | |==================== | 28% | |===================== | 30% | |======================= | 33% | |======================== | 35% | |========================== | 37% | |=========================== | 39% | |============================= | 41% | |============================== | 43% | |================================ | 46% | |================================= | 48% | |=================================== | 50% | |===================================== | 52% | |====================================== | 54% | |======================================== | 57% | |========================================= | 59% | |=========================================== | 61% | |============================================ | 63% | |============================================== | 65% | |=============================================== | 67% | |================================================= | 70% | |================================================== | 72% | |==================================================== | 74% | |===================================================== | 76% | |======================================================= | 78% | |======================================================== | 80% | |========================================================== | 83% | |=========================================================== | 85% | |============================================================= | 87% | |============================================================== | 89% | |================================================================ | 91% | |================================================================= | 93% | |=================================================================== | 96% | |==================================================================== | 98% | |======================================================================| 100%## [1] 1 1 1 5 1 1 1 1 1 1 1 1 2 2 1 1 1 1 1 1 5 1 1 4 1 1 1 1 1 1 2 1 2 1 4 1 1 4
## [39] 1 1 1 1 1 1 2 1hist(munic$class,
main = "Pendientes reclasificadas para el Bolivar",
xlab = "pendiente(categoria)")
En su gran mayoria las pendientes en el departamento son de clase 1
Para poder hacer un mapa para de las pendientes es necesario que las pendientes reclasificadas y sus porcentajes trabajen en coordenas geograficas
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 2811, 1825, 1 (nrow, ncol, nlyr)
## resolution : 0.001362072, 0.001362072 (x, y)
## extent : -76.22416, -73.73838, 6.988779, 10.81756 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 2811, 1825, 1 (nrow, ncol, nlyr)
## resolution : 0.001362072, 0.001362072 (x, y)
## extent : -76.22416, -73.73838, 6.988779, 10.81756 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 127.56667) Mapa de pendiente en porcentaje del departamento
Se realiza una paleta de colores para los porcentajes de pendientes
palredgreen <- colorNumeric(c("green3","yellow2", "darkorange", "purple4", "darkred"), values(slope.geo),
na.color = "transparent") leaflet(munic) %>% addTiles() %>% setView(-75, 8.7, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"clases de pendiente: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),title = "Pendientes del bolivar (%)")## 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'