Atributos geomorfometricos
Thomas Espina
2026-06-13
rm(list=ls())###Calculo de atributos geomorfometricos:
##Destacaremos las librerias elevatr, MultiscaleDTM, exactextractr fundamentales para el analisis de datos geomorfometricos, tenemos que destacar la libreria MultiscaleDTM y sus funciones incluidas ya que nos permitio calcular los atributos buscados. MultiscaleDTM, ES UTILIZADA PARA cálculos de geomorfometría avanzada.
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)
library(sp)##Cargaremos un DEM ráster de nuestro departamento de interes (Cauca). En este cuaderno, utilizaré un DEM obtenido previamente en clase. Este paso resulta funcamental ya que sin este DEM no podria realizar el analisis:
ruta = terra::rast("C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno5\\data\\elev_cauca_z10.tif")
## Se establece que el DEM el minimo valor que pueda tomar sea 0
dem <- terra::clamp(ruta, lower = 0, values = TRUE)## |---------|---------|---------|---------|========================================= ## Vizualozaremos nuestra variable
dem## class : SpatRaster
## size : 4095, 3586, 1 (nrow, ncol, nlyr)
## resolution : 0.0006863432, 0.0006863432 (x, y)
## extent : -78.04687, -75.58565, 0.7028454, 3.513421 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : spat_5aa460503d6a_23204_o8up7yuDdJ6XHBX.tif
## varname : elev_cauca_z10
## name : elev_cauca_z10
## min value : 0
## max value : 5373##En el sigueinte paso se reducira la resolucion del DEM para evitar problemas con la memoria mediante la funcion aggragate que es un funcion que incluye la libreria terra
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= ##Leeremos los municipios del departamento del cauca para tener una clara delimitacion de las dimensiones del dapartamento
munic <- sf::st_read("C:\\Users\\pc\\OneDrive\\Documentos\\GB2\\RSTUDIO\\cuaderno5\\data\\cauca_munic.gpkg")## Reading layer `cauca_munic' from data source
## `C:\Users\pc\OneDrive\Documentos\GB2\RSTUDIO\cuaderno5\data\cauca_munic.gpkg'
## using driver `GPKG'
## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.92834 ymin: 0.9580285 xmax: -75.74782 ymax: 3.328941
## Geodetic CRS: WGS 84##Tanto el DEM como los municipos estan en coordenadas geograficas WGS 84, mediante la funcion crop presente en la libreria terra recortaremos dem2 a nuestro objeto geoespacial que corresponde a la REGION DE INTERES
dem3 = terra::crop(dem2,munic, mask=TRUE)##Transformando las coordenadas: A coordenadas planas debido a la normativa vigente en el pais, utilizaremos el EPSG 9377 MAGNA-SIRGAS – Origen Nacional esto con el objetivo de poder hacer la reproyeccion y calcular los atributos correspondiente
# 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 : 1727, 1598, 1 (nrow, ncol, nlyr)
## resolution : 152.4318, 152.4318 (x, y)
## extent : 4451159, 4694745, 1664036, 1927286 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_cauca_z10
## min value : 0
## max value : 5359.402832Reproyectando los objetos vectoriales ¿Que funcion tiene esta paso dentro del codigo?
(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 42 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4451719 ymin: 1664078 xmax: 4694583 ymax: 1926529
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 19 001 19001 CAUCA POPAYÁN
## 2 19 022 19022 CAUCA ALMAGUER
## 3 19 050 19050 CAUCA ARGELIA
## 4 19 075 19075 CAUCA BALBOA
## 5 19 100 19100 CAUCA BOLÍVAR
## 6 19 110 19110 CAUCA BUENOS AIRES
## 7 19 130 19130 CAUCA CAJIBÍO
## 8 19 137 19137 CAUCA CALDONO
## 9 19 142 19142 CAUCA CALOTO
## 10 19 212 19212 CAUCA CORINTO
## mpio_crslc mpio_tipo mpio_narea mpio_nano shape_Leng
## 1 1537 MUNICIPIO 480.3814 2024 1.374118
## 2 1799 MUNICIPIO 236.8516 2024 0.710253
## 3 Ordenanza 2 de Noviembre 8 de 1967 MUNICIPIO 776.4555 2024 1.224153
## 4 Ordenanza 1 de Octubre 20 de 1967 MUNICIPIO 413.5523 2024 1.016327
## 5 1793 MUNICIPIO 798.3494 2024 1.802787
## 6 1851 MUNICIPIO 435.1190 2024 1.632520
## 7 1824 MUNICIPIO 552.5514 2024 1.381852
## 8 1746 MUNICIPIO 354.8557 2024 1.102961
## 9 1543 MUNICIPIO 264.8714 2024 1.315453
## 10 1868 MUNICIPIO 325.9754 2024 0.828674
## shape_Area geom
## 1 0.03897017 MULTIPOLYGON (((4583956 184...
## 2 0.01919657 MULTIPOLYGON (((4576303 177...
## 3 0.06288063 MULTIPOLYGON (((4520521 183...
## 4 0.03349049 MULTIPOLYGON (((4539082 179...
## 5 0.06468337 MULTIPOLYGON (((4562832 179...
## 6 0.03530827 MULTIPOLYGON (((4581114 191...
## 7 0.04481975 MULTIPOLYGON (((4573051 186...
## 8 0.02880277 MULTIPOLYGON (((4626166 187...
## 9 0.02150781 MULTIPOLYGON (((4619907 191...
## 10 0.02648009 MULTIPOLYGON (((4641747 189...Aqui haremos uso de la libreria MultiscaleDTM con el objetivo de calcular los atributos correspondientes, usaremos una venta 3x3 aplicaremos el metodo QUEEN y correguiremos los valores que tenga una pendiante = 0
(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 : 1727, 1598, 2 (nrow, ncol, nlyr)
## resolution : 152.4318, 152.4318 (x, y)
## extent : 4451159, 4694745, 1664036, 1927286 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0, 0.000018
## max values : 68.76312, 359.999962¿Que significa que los valores obtenidos de aspect y slope den codigo?
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1727, 1598, 1 (nrow, ncol, nlyr)
## resolution : 152.4318, 152.4318 (x, y)
## extent : 4451159, 4694745, 1664036, 1927286 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 68.76312## ¿Que funcion tiene este fragmento de codigo? terra::hist(slope,
main = "Cauca's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 36% of the cells was used (of which 51% was NA)
## ¿Que interpretacion puede tener el siguiente histograma?
## Posible interpretacion: El cauca tiene regiones mas planas que con alta pendiente segun el histograma (aspect =subset(slp_asp,2 ))## class : SpatRaster
## size : 1727, 1598, 1 (nrow, ncol, nlyr)
## resolution : 152.4318, 152.4318 (x, y)
## extent : 4451159, 4694745, 1664036, 1927286 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.000018
## max value : 359.999962## ¿Que funcion tiene terra::hist(aspect,
main = "Cauca's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 36% of the cells was used (of which 51% was NA)
Como interpretar el aspecto de esta grafica?
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1727, 1598, 1 (nrow, ncol, nlyr)
## resolution : 152.4318, 152.4318 (x, y)
## extent : 4451159, 4694745, 1664036, 1927286 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 257.323996terra::hist(slope_perc,
main = "Cauca's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 36% of the cells was used (of which 51% was NA)
# 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)(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG: 9377) .
## To avoid on-the-fly coordinate transformation, ensure that st_crs() returns the
## same result for both raster and polygon inputs.## | | | 0% | |== | 2% | |=== | 5% | |===== | 7% | |======= | 10% | |======== | 12% | |========== | 14% | |============ | 17% | |============= | 19% | |=============== | 21% | |================= | 24% | |================== | 26% | |==================== | 29% | |====================== | 31% | |======================= | 33% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |================================ | 45% | |================================= | 48% | |=================================== | 50% | |===================================== | 52% | |====================================== | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |=============================================== | 67% | |================================================ | 69% | |================================================== | 71% | |==================================================== | 74% | |===================================================== | 76% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 86% | |============================================================== | 88% | |=============================================================== | 90% | |================================================================= | 93% | |=================================================================== | 95% | |==================================================================== | 98% | |======================================================================| 100%## [1] 12.7062006 34.4478607 42.5884094 34.2758980 27.2228680 21.2370110
## [7] 15.2377043 21.4560223 18.5817146 29.1378345 31.0206146 24.0970097
## [13] 1.0455663 12.4084435 26.5436401 33.1763306 23.1620789 32.6481247
## [19] 20.0335407 18.0960045 23.0177364 26.4700069 0.6239001 35.3344116
## [25] 18.5916405 15.1034260 13.3238077 0.5157506 19.4977303 24.2509060
## [31] 24.4898586 11.2066650 30.5933552 23.7175999 23.4066830 28.3654919
## [37] 28.6634159 12.7001238 13.5686045 36.8871307 17.3526478 0.6462926hist(munic$mean_slope,
main = "Cordoba's mean slope",
xlab = "Slope (in percentage)")
Calcularemos el valor promedio de la pendiente
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG: 9377) .
## To avoid on-the-fly coordinate transformation, ensure that st_crs() returns the
## same result for both raster and polygon inputs.## | | | 0% | |== | 2% | |=== | 5% | |===== | 7% | |======= | 10% | |======== | 12% | |========== | 14% | |============ | 17% | |============= | 19% | |=============== | 21% | |================= | 24% | |================== | 26% | |==================== | 29% | |====================== | 31% | |======================= | 33% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |================================ | 45% | |================================= | 48% | |=================================== | 50% | |===================================== | 52% | |====================================== | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |=============================================== | 67% | |================================================ | 69% | |================================================== | 71% | |==================================================== | 74% | |===================================================== | 76% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 86% | |============================================================== | 88% | |=============================================================== | 90% | |================================================================= | 93% | |=================================================================== | 95% | |==================================================================== | 98% | |======================================================================| 100%## [1] 4 5 5 5 5 5 2 5 5 5 5 5 1 2 5 5 4 5 5 4 1 5 1 5 5 1 2 1 4 4 5 1 5 4 4 5 5 2
## [39] 2 5 4 1hist(munic$class,
main = "Cauca's reclassified slope",
xlab = "Slope (as a category)")
Reclasifficaion de la pendiente del cauca.
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1736, 1597, 1 (nrow, ncol, nlyr)
## resolution : 0.00137266, 0.00137266 (x, y)
## extent : -77.93644, -75.7443, 0.955096, 3.338034 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 216.17717(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1736, 1597, 1 (nrow, ncol, nlyr)
## resolution : 0.00137266, 0.00137266 (x, y)
## extent : -77.93644, -75.7443, 0.955096, 3.338034 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 235.120575# 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")leaflet(munic) %>% addTiles() %>% setView(-76.8259,2.7049, 9) %>%
addPolygons(color = "gray", 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 Cordoba (%)")