Geomorphometric terrain attributes
Jair Leandro Diaz Aldana;Ingenieria Agronómica
11/06/2025
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.13.1, GDAL 3.11.0, PROJ 9.6.0; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.8.54library(MultiscaleDTM)
library(exactextractr)list.files("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE")## [1] "CASANARE.gpkg" "elev_CASANARE_z10.tif"(dem = terra::rast("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/elev_CASANARE_z10.tif"))## class : SpatRaster
## size : 3573, 5128, 1 (nrow, ncol, nlyr)
## resolution : 0.0006856127, 0.0006856127 (x, y)
## extent : -73.125, -69.60918, 4.214913, 6.664608 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_CASANARE_z10.tif
## name : elev_CASANARE_z10
## min value : -789
## max value : 5332dem2 = terra::aggregate(dem,2,"mean")## |---------|---------|---------|---------|========================================= |———|———|———|———|=========================================
(munic <- sf::st_read("C:/Users/Janus/Desktop/GB2/proyecto3/datos/departamentos/CASANARE/CASANARE.gpkg"))## Reading layer `municipios' from data source
## `C:\Users\Janus\Desktop\GB2\proyecto3\datos\departamentos\CASANARE\CASANARE.gpkg'
## using driver `GPKG'
## Simple feature collection with 19 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -73.07777 ymin: 4.287476 xmax: -69.83591 ymax: 6.346111
## Geodetic CRS: MAGNA-SIRGAS(dem3 = terra::crop(dem2,munic, mask=TRUE))## Warning: [crop] CRS do not match## class : SpatRaster
## size : 1502, 2365, 1 (nrow, ncol, nlyr)
## resolution : 0.001371225, 0.001371225 (x, y)
## extent : -73.07838, -69.83543, 4.286903, 6.346483 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source(s) : memory
## name : elev_CASANARE_z10
## min value : -89.5
## max value : 4423.0(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_CASANARE_z10
## min value : -8.611589
## max value : 4416.087891(munic_plane = sf::st_transform(munic, "EPSG:9377"))(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 : 1506, 2370, 2 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 2.159213e-05, 4.141929e-04
## max values : 5.750944e+01, 3.599996e+02(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 2.159213e-05
## max value : 5.750944e+01terra::hist(slope,
main = "Casanare's slope",
xlab = "Slope (in degrees")## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 4.141929e-04
## max value : 3.599996e+02terra::hist(aspect,
main = "Casanare's aspect",
xlab = " Aspect (in degrees)")## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1506, 2370, 1 (nrow, ncol, nlyr)
## resolution : 151.8373, 151.8373 (x, y)
## extent : 4991306, 5351161, 2031677, 2260344 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 3.768537e-05
## max value : 1.570256e+02terra::hist(slope_perc,
main = "Casanare's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 28% of the cells was used (of which 46% was NA)
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)## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%## [1] 5.9245358 9.2382898 30.9212017 1.4186283 40.3001289 0.5959594
## [7] 9.1764860 4.7089849 0.5160555 0.7966918 2.8408487 28.1054153
## [13] 15.5886774 28.3318501 0.5730320 25.7137718 4.7331200 0.6306346
## [19] 1.4323320Warning 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] 2.5160389 1.2292731 1.4036171 5.2911134 0.9474660 0.5868814 [7] 1.5853550 3.0762432 0.3953781 1.7273880 2.0587800 4.6855459 [13] 1.4440339 10.1471510 3.8976312 3.0785139 2.0715215 3.6382058 [19] 16.4352589 2.8243833 2.1665022 2.3823540 3.1995339 1.3972237 [25] 3.4256206 11.3830862 1.5984520 14.8620186 2.7408032 5.6827564
hist(munic$mean_slope,
main = "Casanare´s slope",
xlab = "Slope (in percentage")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |==== | 5% | |======= | 11% | |=========== | 16% | |=============== | 21% | |================== | 26% | |====================== | 32% | |========================== | 37% | |============================= | 42% | |================================= | 47% | |===================================== | 53% | |========================================= | 58% | |============================================ | 63% | |================================================ | 68% | |==================================================== | 74% | |======================================================= | 79% | |=========================================================== | 84% | |=============================================================== | 89% | |================================================================== | 95% | |======================================================================| 100%## [1] 1 1 5 1 5 1 1 1 1 1 1 5 1 5 1 5 1 1 1Warning 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] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 4 1 1
hist(munic$class,
main = "Casanare´s reclassified slope",
xlab = "Slope (as a category)")
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 7(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1514, 2373, 1 (nrow, ncol, nlyr)
## resolution : 0.001371239, 0.001371239 (x, y)
## extent : -73.07864, -69.82469, 4.280104, 6.35616 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 2.620884e-03
## max value : 1.570256e+02palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")leaflet(munic) %>% addTiles() %>% setView(-72.0, 5.3, 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 (%)")## 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'This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.