ELEVACION
Maria Alejandra
true
Intruducción
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.
- Setup
rm(list=ls())library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(exactextractr)
library(MultiscaleDTM)Leer el modelo de elevacion digital
(dem = terra::rast("ELEVACION/ANTIO_dem.tif"))## class : SpatRaster
## dimensions : 413, 389, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -77.125, -73.88333, 5.425, 8.866667 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : ANTIO_dem.tif
## name : ANTIO_demdem2 = terra::aggregate(dem,2, "mean")(munic <- sf::st_read("ELEVACION/MUN_ANTI.gpkg"))## Reading layer `mgn_adm_mpio_grafico' from data source
## `C:\Users\Alejandra\Documents\pf\PROYECTO4\DATOS\ELEVACION\MUN_ANTI.gpkg'
## using driver `GPKG'
## Simple feature collection with 125 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
## Geodetic CRS: MAGNA-SIRGAS## Simple feature collection with 125 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -77.12783 ymin: 5.418558 xmax: -73.88128 ymax: 8.873974
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp mpio_cnmbr Area_km2
## 1 05 001 05001 MEDELLÍN 374.56650
## 2 05 002 05002 ABEJORRAL 506.85538
## 3 05 004 05004 ABRIAQUÍ 296.53170
## 4 05 021 05021 ALEJANDRÍA 128.89200
## 5 05 030 05030 AMAGÁ 84.06465
## 6 05 031 05031 AMALFI 1208.84471
## 7 05 034 05034 ANDES 402.04874
## 8 05 036 05036 ANGELÓPOLIS 81.80967
## 9 05 038 05038 ANGOSTURA 338.33405
## 10 05 040 05040 ANORÍ 1413.31828
## geom
## 1 MULTIPOLYGON (((-75.66974 6...
## 2 MULTIPOLYGON (((-75.46938 5...
## 3 MULTIPOLYGON (((-76.08351 6...
## 4 MULTIPOLYGON (((-75.0332 6....
## 5 MULTIPOLYGON (((-75.67587 6...
## 6 MULTIPOLYGON (((-74.91836 7...
## 7 MULTIPOLYGON (((-75.86822 5...
## 8 MULTIPOLYGON (((-75.69149 6...
## 9 MULTIPOLYGON (((-75.27173 6...
## 10 MULTIPOLYGON (((-74.90935 7...(dem3 = terra::crop(dem2,munic, mask=TRUE))## class : SpatRaster
## dimensions : 207, 195, 1 (nrow, ncol, nlyr)
## resolution : 0.01666667, 0.01666667 (x, y)
## extent : -77.125, -73.875, 5.416667, 8.866667 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : ANTIO_dem
## min value : 0
## max value : 3761- 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?
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 208, 196, 1 (nrow, ncol, nlyr)
## resolution : 1842.536, 1842.536 (x, y)
## extent : 4542819, 4903956, 2157042, 2540290 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : ANTIO_dem
## min value : 0.6166098
## max value : 3656.6940918We also need to conduct a similar transformation for the vector object:
(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 125 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4544725 ymin: 2157214 xmax: 4902701 ymax: 2540303
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp mpio_cnmbr Area_km2
## 1 05 001 05001 MEDELLÍN 374.56650
## 2 05 002 05002 ABEJORRAL 506.85538
## 3 05 004 05004 ABRIAQUÍ 296.53170
## 4 05 021 05021 ALEJANDRÍA 128.89200
## 5 05 030 05030 AMAGÁ 84.06465
## 6 05 031 05031 AMALFI 1208.84471
## 7 05 034 05034 ANDES 402.04874
## 8 05 036 05036 ANGELÓPOLIS 81.80967
## 9 05 038 05038 ANGOSTURA 338.33405
## 10 05 040 05040 ANORÍ 1413.31828
## geom
## 1 MULTIPOLYGON (((4704762 226...
## 2 MULTIPOLYGON (((4726714 221...
## 3 MULTIPOLYGON (((4659221 230...
## 4 MULTIPOLYGON (((4775207 226...
## 5 MULTIPOLYGON (((4703922 223...
## 6 MULTIPOLYGON (((4788293 236...
## 7 MULTIPOLYGON (((4682435 219...
## 8 MULTIPOLYGON (((4702253 224...
## 9 MULTIPOLYGON (((4749105 232...
## 10 MULTIPOLYGON (((4789365 238...#atributos=SlpAsp(MultiscaleDTM)
(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 : 208, 196, 2 (nrow, ncol, nlyr)
## resolution : 1842.536, 1842.536 (x, y)
## extent : 4542819, 4903956, 2157042, 2540290 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 9.006649e-04, 5.220056e-03
## max values : 1.792480e+01, 3.599943e+02(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 208, 196, 1 (nrow, ncol, nlyr)
## resolution : 1842.536, 1842.536 (x, y)
## extent : 4542819, 4903956, 2157042, 2540290 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 9.006649e-04
## max value : 1.792480e+01terra::hist(slope,
main = "Antioquia's slope",
xlab = "Slope (in degrees)")
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 208, 196, 1 (nrow, ncol, nlyr)
## resolution : 1842.536, 1842.536 (x, y)
## extent : 4542819, 4903956, 2157042, 2540290 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 5.220056e-03
## max value : 3.599943e+02terra::hist(aspect,
main = "Antioquia's aspect",
xlab = "Aspect (in degrees)")
#
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## dimensions : 208, 196, 1 (nrow, ncol, nlyr)
## resolution : 1842.536, 1842.536 (x, y)
## extent : 4542819, 4903956, 2157042, 2540290 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.001571957
## max value : 32.346916759terra::hist(slope_perc,
main = "Antioquia's slope",
xlab = "Slope (in percentage)")
#terra::hist(slope_perc)# 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)## | | | 0% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 10.7041664 10.6345587 8.9466019 4.0995612 9.4517746 6.2421207
## [7] 11.9016571 11.1935673 5.1079688 6.1736937 10.8155918 12.6221418
## [13] 3.5149839 1.0484123 12.0423212 8.1422052 14.4538040 5.7812872
## [19] 11.2842922 11.1984520 10.5730038 10.0931597 14.3411341 15.0252104
## [25] 1.6770254 11.3654308 9.3151522 9.3139639 10.1389380 2.8473873
## [31] 15.6334209 1.2593139 8.2670517 5.1766248 0.3621982 2.6534741
## [37] 7.8004942 8.8749819 4.7108736 12.2691402 14.6201429 8.8845243
## [43] 6.3628283 10.0800400 2.9797909 3.5400000 9.5500212 10.4934978
## [49] 7.7424369 10.2127571 14.3926783 7.6983910 7.6615405 9.4228868
## [55] 3.9440317 3.2128665 12.4204750 8.1694794 5.7525082 11.9966402
## [61] 10.3795080 9.8752403 6.4867201 10.7288980 6.7359338 2.2791073
## [67] 16.9922199 2.9432967 1.8123683 10.7881680 3.6968403 7.1056490
## [73] 14.0619335 0.8601292 1.1075258 18.1658783 3.8431520 15.5554171
## [79] 8.6265011 2.2281306 2.3938746 1.1341115 2.4941216 5.1679792
## [85] 2.2748320 16.4687309 9.6522436 9.4759169 13.1577377 5.7235403
## [91] 6.3824172 15.6645098 4.9916821 0.9038266 6.2976727 3.2816935
## [97] 0.9955718 6.0484014 4.0532765 2.9090943 8.5258884 4.0131092
## [103] 7.3747005 2.9765301 2.7065382 7.7175527 12.1313829 14.6914682
## [109] 5.4275265 10.3890467 11.6484213 14.2961197 1.2175245 9.9546728
## [115] 7.7310972 9.5903368 9.4232435 2.8154550 9.5543423 3.3225145
## [121] 2.5806246 7.5960279 4.0361981 0.3794824 1.4480009## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)hist(munic$mean_slope,
main = "Antioquia'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% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |=== | 4% | |=== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |=========== | 15% | |=========== | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================= | 24% | |================= | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |========================= | 35% | |========================= | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |=============================== | 44% | |=============================== | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |======================================= | 55% | |======================================= | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================= | 64% | |============================================= | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 4 4 3 2 3 2 4 3 2 2 4 4 1 1 4 3 4 2 4 4 3 3 4 4 1 4 4 3 3 1 4 1 3 2 1 1 3
## [38] 3 2 4 4 4 2 3 1 1 4 4 2 3 4 3 2 4 2 1 4 2 2 4 4 2 2 4 3 1 4 1 1 3 1 1 4 1
## [75] 1 4 1 4 3 1 1 1 1 2 1 4 2 3 4 2 2 4 2 1 2 1 1 2 2 1 3 1 2 1 1 2 4 4 2 3 4
## [112] 4 1 4 2 3 4 1 3 1 1 3 1 1 1## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)hist(munic$class,
main = "Antioquia's reclassified slope",
xlab = "Slope (as a category)")
#terra::hist(slope_perc)# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 209, 197, 1 (nrow, ncol, nlyr)
## resolution : 0.01666664, 0.01666664 (x, y)
## extent : -77.15609, -73.87276, 5.405126, 8.888454 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 5# This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 209, 197, 1 (nrow, ncol, nlyr)
## resolution : 0.01666664, 0.01666664 (x, y)
## extent : -77.15609, -73.87276, 5.405126, 8.888454 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.003827545
## max value : 30.566703796# 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(-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 Antioquia (%)")## 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'