Geomorphometric terrain attributes in R
Ivan Lizarazo
22.01.2025
- 1. Introduction
- 2. Setup
- 3. Introduction to MultiscaleDTM
- 4. Reading input data
- 5. Transforming coordinates (10%)
- 6. Computing terrain attributes (20%)
- 7. Computing zonal statistics (20%)
- 8. Plotting terrain slope (20%)
- 9. Alternative visualization of terrain slope (15%)
- 10. Visualizing terrain aspect (15%)
- Bibliography
1. Introduction
This notebook illustrates how to calculate geomorphometric terrain attributes from gridded digital elevation models. It aims to teach basic geoinformatics concepts & tools to undergraduate agronomy students at Universidad Nacional de Colombia.
For calculating terrain attributes, we will use the MultiscaleDTM R package.
When creating your notebook, make sure to add as many text as needed in order to: (i) explain each code chunk; and (ii) describe the corresponding output. These two criteria will be considered to assessing the quality of your notebook.
2. Setup
The required packages should be installed in advance from the Console.
## SETUP
# INSTALL THIS PACKAGE FROM THE CONSOLE, NOT FROM THIS CHUNK
# paquetes = c("MultiscaleDTM", "exactextractr")
# install.packages(paquetes)Let’s start by cleaning the memory:
rm(list=ls())Now, let’s load the required libraries:
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)When running a chunk, pay attention to any message indicating error or warning eventual conflicts.
Conflicts may be caused by functions, from differents packages, which have the same names (e.g. terra::extract and tidyr::extract). In order to avoid such problems, we need to call several functions using the long way, i.e. package_name::function.
3. Introduction to MultiscaleDTM
This R package calculates multi-scale geomorphometric terrain attributes from regularly gridded digital terrain models (DTM; i.e. elevation or bathymetry rasters) via a specified window.
In case you need to refresh how to compute geomorphometric attributes, review the course slides.
4. Reading input data
First, we will load a raster DEM for our department. In this notebook, I will use a DEM previously obtained using this notebook.
Let’s read the DEM using terra:
# change the path & DEM filename as needed
(dem = terra::rast("cordoba/elev_cordoba_z10.tif"))## class : SpatRaster
## dimensions : 3568, 3090, 1 (nrow, ncol, nlyr)
## resolution : 0.000682546, 0.000682546 (x, y)
## extent : -76.64062, -74.53156, 7.013738, 9.449062 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : elev_cordoba_z10.tif
## name : elev_cordoba_z10
## min value : -1241
## max value : 3724Let’s reduce DEM resolution in order avoid memory issues:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Now, we will read our municipalities using the sf library:
(munic <- sf::st_read("cordoba/cordoba.gpkg"))## Reading layer `mgn_adm_mpio_grafico' from data source
## `/Users/ials/Documents/unal/GB2024/data/cordoba/cordoba.gpkg'
## using driver `GPKG'
## Simple feature collection with 30 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.51453 ymin: 7.347087 xmax: -74.78095 ymax: 9.447748
## Geodetic CRS: MAGNA-SIRGASNote that both datasets are in geographic coordinates.
Now, let’s clip the dem2 object to our department’s boundaries:
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Transforming coordinates (10%)
In order to compute terrain attributes, the DEM needs to be in plane coordinates. As you may remember, IGAC changed in 2020 to a sistema de coordenadas planas con origen unico nacional. What is its EPSG code?
Modify the following two code chunks to transform the input data to the Colombia’s national coordinate system.
# 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
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_cordoba_z10
## min value : -58.53565
## max value : 2688.35400We also need to conduct a similar transformation for the vector object:
# check the parameters of this sf function here
# https://github.com/rstudio/cheatsheets/blob/main/sf.pdf
(munic_plane = sf::st_transform(munic, "EPSG:9377"))6. Computing terrain attributes (20%)
Let’s find how to compute slope and aspect.
Go to the Console and type ?SlpAsp. Then, in the Help window to the right, the corresponding info will display:
Now, let’s call the SlpAsp function:
# Explain what is w
# Explain what is method
# Change if needed
(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
## dimensions : 1550, 1277, 2 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0000, 1.234529e-04
## max values : 57.1558, 3.599994e+02Now, let’s split the slp_asp object in two parts.
Layer 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 57.1558What is the slope values distribution?
terra::hist(slope,
main = "Cordoba's slope",
xlab = "Slope (in degrees)")
Layer 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.234529e-04
## max value : 3.599994e+02terra::hist(aspect,
main = "Cordoba's aspect",
xlab = "Aspect (in degrees)")
In addition, let’s convert slope degrees to slope percentage:
#
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## dimensions : 1550, 1277, 1 (nrow, ncol, nlyr)
## resolution : 150.7599, 150.7599 (x, y)
## extent : 4611972, 4804493, 2370248, 2603926 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.000
## max value : 154.907terra::hist(slope_perc,
main = "Cordoba's slope",
xlab = "Slope (in percentage)")
#terra::hist(slope_perc)7. Computing zonal statistics (20%)
A zonal statistics operation is one that calculates statistics on cell values of a raster (a value raster) within the zones defined by another dataset. In our case, we are interested in computing average of slope values per municipality.
First, let’s know the IGAC’s slope classification for agrology purposes:
In case you want to understand rationale behind this classification check this link
Now, let’s put our hands to work:
# 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)Now, let’s compute zonal statistics using exactextractr:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning 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.6827564hist(munic$mean_slope,
main = "Cordoba's mean 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% | |== | 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 1hist(munic$class,
main = "Cordoba's reclassified slope",
xlab = "Slope (as a category)")
#terra::hist(slope_perc)Let’s transform the slope back to geographic coordinates:
# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1556, 1291, 1 (nrow, ncol, nlyr)
## resolution : 0.001365132, 0.001365132 (x, y)
## extent : -76.53401, -74.77163, 7.33663, 9.460775 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7# This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1556, 1291, 1 (nrow, ncol, nlyr)
## resolution : 0.001365132, 0.001365132 (x, y)
## extent : -76.53401, -74.77163, 7.33663, 9.460775 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.000
## max value : 133.0098. Plotting terrain slope (20%)
First, we will need a color palette for slope:
# 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")Now, it is plotting time. Make sure you change several options in the code if you need it:
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 = "Terrain slope in Cordoba (%)")Click on different sites for getting municipalities’ names and most frequent percentage slope values at each.
9. Alternative visualization of terrain slope (15%)
Visualize slope terrain in degrees and add labels showing municipalities names along their corresponding mean slope values. Keep the popup values as depicted in the previous plot.
# WRITE YOUR CODE IN THIS CODE CHUNK
#
#
#Describe the outcome and its meaning.
10. Visualizing terrain aspect (15%)
Visualize terrain aspect using a palette similar to this one.
# WRITE YOUR CODE IN THIS CODE CHUNK
#
#
#Describe the outcome and its meaning.
Once your notebook is complete, publish it at rpubs, review its contents, fix any errors and/or make any changes needed to demonstrate that you are proficient in computing & visualizing terrain attributes.
Send by e-mail the link to your published notebook before 1st February 2025 at 11:59 pm
Bibliography
If you reuse this notebook’s code please cite this work as:
Lizarazo, I., 2025. Geomorphometric terrain attributes in R. Available at https://rpubs.com/ials2un/geomorphometric
Reproducibility
sessionInfo()## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Sonoma 14.7.1
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/Bogota
## tzcode source: internal
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] exactextractr_0.10.0 MultiscaleDTM_0.8.3 terra_1.7-78
## [4] leaflet_2.2.1 sf_1.0-17 elevatr_0.99.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.4 sass_0.4.9 generics_0.1.3 class_7.3-22
## [5] KernSmooth_2.23-22 lattice_0.22-6 digest_0.6.37 magrittr_2.0.3
## [9] rgl_1.1.3 evaluate_1.0.0 grid_4.3.2 fastmap_1.2.0
## [13] jsonlite_1.8.9 e1071_1.7-16 DBI_1.2.3 promises_1.3.0
## [17] purrr_1.0.2 fansi_1.0.6 scales_1.3.0 crosstalk_1.2.1
## [21] codetools_0.2-19 jquerylib_0.1.4 cli_3.6.3 shiny_1.9.1
## [25] rlang_1.1.4 units_0.8-5 munsell_0.5.1 base64enc_0.1-3
## [29] cachem_1.1.0 yaml_2.3.10 raster_3.6-26 tools_4.3.2
## [33] dplyr_1.1.4 colorspace_2.1-1 httpuv_1.6.15 png_0.1-8
## [37] vctrs_0.6.5 R6_2.5.1 mime_0.12 proxy_0.4-27
## [41] lifecycle_1.0.4 classInt_0.4-10 htmlwidgets_1.6.4 pkgconfig_2.0.3
## [45] progressr_0.14.0 bslib_0.8.0 pillar_1.9.0 later_1.3.2
## [49] glue_1.8.0 Rcpp_1.0.13 highr_0.11 tidyselect_1.2.1
## [53] tibble_3.2.1 xfun_0.47 rstudioapi_0.16.0 knitr_1.48
## [57] farver_2.1.2 xtable_1.8-4 htmltools_0.5.8.1 rmarkdown_2.28
## [61] compiler_4.3.2 sp_2.1-4