Atributos geomorfométricos del terreno
Laura Cortez
10 de junio
##Introducción En este cuaderno con el cálculo previo de un modelo de elevación (DEM) y la ayuda de librerías como MultiscaleDTM y elevatr, se logrará calcular la pendiente y la orientación del departamento de Bolívar.
##Configuración
rm(list = ls())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.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.9.27library(MultiscaleDTM)
library(exactextractr)##Lectura de datos de entrada
getwd()## [1] "C:/Users/imthe/Documents"list.files()## [1] "ArcGISEarth"
## [2] "Bolívar.gpkg"
## [3] "Bolivar_elevacion.tif"
## [4] "Bolivar_líneas de agua_.gpkg"
## [5] "Certificación laboral Pedro Tangua Conductor SCT477.pdf"
## [6] "Cuaderno-de-atributos-geomorfométricos.Rmd"
## [7] "Cuaderno-de-atributos-geomorfométricos_files"
## [8] "Cuaderno-de-Soil-Grids.knit.md"
## [9] "Cuaderno-de-Soil-Grids.Rmd"
## [10] "Cuaderno-Sentinel-2.html"
## [11] "Cuaderno 1 _ Informe de geomática.Rmd"
## [12] "Cuaderno de atributos geomorfométricos.Rmd"
## [13] "Cuaderno de Cálculos con Dem.Rmd"
## [14] "Cuaderno de Interpolación.nb.html"
## [15] "Cuaderno de Interpolación.Rmd"
## [16] "Cuaderno de Soil Grids.nb.html"
## [17] "Cuaderno de Soil Grids.Rmd"
## [18] "Cuaderno Sentinel-2.nb.html"
## [19] "Cuaderno Sentinel-2.Rmd"
## [20] "Cuaderno sobre centroides.nb.html"
## [21] "Descargas - Acceso directo.lnk"
## [22] "desktop.ini"
## [23] "Formato Responsabilidad IndividualSalAcad- Laura Cortez.pdf"
## [24] "Formato Salidas Académicas-Laura Cortez.pdf"
## [25] "Formato Salidas Académicas-Laura Cortez.xlsx"
## [26] "GB2"
## [27] "InterpolaciónSoilGrids.html"
## [28] "InterpolaciónSoilGrids.Rmd"
## [29] "InterpolaciónSoilGrids_files"
## [30] "map.geojson"
## [31] "MGN2025_MPIO_GRAFICO"
## [32] "Mi música"
## [33] "Miprimercuaderno.html"
## [34] "Miprimercuaderno.nb.html"
## [35] "Miprimercuaderno.Rmd"
## [36] "Mis imágenes"
## [37] "Mis vídeos"
## [38] "Plantillas personalizadas de Office"
## [39] "Presentación Caracterización Finca Normandía versión 2.pptx"
## [40] "rsconnect"
## [41] "Santa_Rosa_del_Sur"
## [42] "Santa_Rosa_delSur"
## [43] "World _Cities"
## [44] "Zona de interés- Santa Rosa del Sur.geojson"list.files(pattern = "\\.tif$|\\.tiff$")## [1] "Bolivar_elevacion.tif"(dem = terra::rast ("Bolivar_elevacion.tif"))## class : SpatRaster
## size : 6125, 3616, 1 (nrow, ncol, nlyr)
## resolution : 0.0006805766, 0.0006805766 (x, y)
## extent : -75.9375, -73.47654, 6.664774, 10.83331 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : Bolivar_elevacion.tif
## name : file70349d6601e
## min value : -2029
## max value : 5578dem2= terra::aggregate(dem, 2, "mean")## |---------|---------|---------|---------|========================================= (munic <- sf :: st_read("Bolívar.gpkg"))## Reading layer `Bolívar' from data source `C:\Users\imthe\Documents\Bolívar.gpkg' using driver `GPKG'
## Simple feature collection with 6 features and 90 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -75.79552 ymin: 6.979202 xmax: -73.769 ymax: 10.80139
## Geodetic CRS: MAGNA-SIRGASdem3 = terra::crop(dem2, munic, mask=TRUE)## Warning: [crop] CRS do not match##Transformación de coordenadas
dem_plane = project(dem3, "EPSG:9377")munic_plane = sf::st_transform(munic, "EPSG:9377")##Computación de atributos de terreno
(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 : 2819, 1495, 2 (nrow, ncol, nlyr)
## resolution : 150.3149, 150.3149 (x, y)
## extent : 4691176, 4915897, 2329319, 2753057 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0, 0
## max values : 58.195922, 359.999709Primera capa
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 2819, 1495, 1 (nrow, ncol, nlyr)
## resolution : 150.3149, 150.3149 (x, y)
## extent : 4691176, 4915897, 2329319, 2753057 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 58.195922terra::hist(slope,
main = "Pendiente de Bolívar",
xlab = "Pendiente (en grados")## Warning: [hist] a sample of 24% of the cells was used (of which 73% was NA)
Segunda Capa
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 2819, 1495, 1 (nrow, ncol, nlyr)
## resolution : 150.3149, 150.3149 (x, y)
## extent : 4691176, 4915897, 2329319, 2753057 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0
## max value : 359.999709terra::hist(aspect,
main = "Orientación de Bolívar",
xlab = "Orientación (en grados)")## Warning: [hist] a sample of 24% of the cells was used (of which 73% was NA)
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 2819, 1495, 1 (nrow, ncol, nlyr)
## resolution : 150.3149, 150.3149 (x, y)
## extent : 4691176, 4915897, 2329319, 2753057 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0
## max value : 161.257861terra::hist(slope_perc,
main = "Pendiente de Bolívar",
xlab = "Pendiente (en porcentaje)")## Warning: [hist] a sample of 24% of the cells was used (of which 73% was NA)
##Computación de Estadísticas Zonales
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% | |============ | 17% | |======================= | 33% | |=================================== | 50% | |=============================================== | 67% | |========================================================== | 83% | |======================================================================| 100%## [1] 1.1860836 0.6501005 9.1903152 1.2698051 7.9200978 12.3590889hist(munic$mean_slope,
main = "Pendiente de Bolívar",
xlab = "Pendiente (en porcentajes)")
(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% | |============ | 17% | |======================= | 33% | |=================================== | 50% | |=============================================== | 67% | |========================================================== | 83% | |======================================================================| 100%## [1] 1 1 1 1 3 4hist(munic$class,
main = "Pendiente reclasificada de Bolívar",
xlab = "Pendiente (como una categoría)")
terra :: hist (slope_perc)## Warning: [hist] a sample of 24% of the cells was used (of which 73% was NA)
(rc.geo = project (rc, "EPSG:4326"))## class : SpatRaster
## size : 2822, 1516, 1 (nrow, ncol, nlyr)
## resolution : 0.001361191, 0.001361191 (x, y)
## extent : -75.82511, -73.76154, 6.972412, 10.81369 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 77.514595(slope.geo = project(slope_perc, "EPSG: 4326"))## class : SpatRaster
## size : 2822, 1516, 1 (nrow, ncol, nlyr)
## resolution : 0.001361191, 0.001361191 (x, y)
## extent : -75.82511, -73.76154, 6.972412, 10.81369 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 140.157745##Trazado de la pendiente del terreno (20%)
palredgreen <- colorNumeric (c("darkseagreen3", "yellow2","orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")slope.deg <- terra::terrain(dem, v = "slope", unit = "degrees")
slope.pct <- tan(slope.deg * pi / 180) * 100
slope.small <- terra::aggregate(slope.geo, fact = 8)
leaflet(munic) %>% addTiles() %>% setView(-75.9, 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>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Pendiente del terreno en Bolívar (%)")## 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'##Visualización alterna de la pendiente del terreno
slope.geo <- terra::terrain(dem, v = "slope", unit = "degrees")
slope.small <- terra::aggregate(slope.geo, fact = 4)## |---------|---------|---------|---------|========================================= leaflet(munic) %>%
addTiles() %>%
setView(-75.9, 8.7, 9) %>%
addRasterImage(
slope.small,
colors = palredgreen,
opacity = 0.8
) %>%
addPolygons(
color = "gray",
weight = 1,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste(
"Municipio:", munic$mpio_cnmbr, "<br>",
"Slope class:", munic$class, "<br>"
),
label = paste0(
munic$mpio_cnmbr,
" (",
round(munic$mean_slope, 2),
"°)"
)
) %>%
addLegend(
pal = palredgreen,
values = values(slope.small),
title = "Pendiente del terreno (°)"
)## 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'##Visualización de la orientación del terreno
aspect.geo <- terra::terrain(dem,
v = "aspect",
unit = "degrees")
pal.aspect <- colorNumeric(
palette = c(
"#d73027", # Norte
"#fc8d59", # Noreste
"#fee08b", # Este
"#d9ef8b", # Sureste
"#91cf60", # Sur
"#66c2a5", # Suroeste
"#3288bd", # Oeste
"#5e4fa2", # Noroeste
"#d73027" # Norte nuevamente para cerrar el círculo
),
domain = c(0, 360),
na.color = "transparent"
)aspect.geo <- terra::terrain(dem, v = "aspect", unit = "degrees")
aspect.small <- terra::aggregate(aspect.geo, fact = 4)## |---------|---------|---------|---------|========================================= # Paleta de colores
pal.aspect <- colorNumeric(
palette = terrain.colors(10),
domain = terra::values(aspect.small),
na.color = "transparent"
)
# Mapa
leaflet() %>%
addTiles() %>%
addRasterImage(
aspect.small,
colors = pal.aspect,
opacity = 0.8
) %>%
addLegend(
pal = pal.aspect,
values = terra::values(aspect.small),
title = "Orientación del terreno (°)"
)