Atributos geomorfométricos del terreno en R
Karenth Yuliana Menjura Coy
#Tabla de Contenido 1. Introducción 2. Configuración 3. Introducción a MultiscaleDTM 4. Lectura de datos de entrada 5. Transformación de coordenadas (10%) 6. Cálculo de atributos del terreno (20%) 7. Cálculo de estadísticas zonales (20%) 8. Trazado de la pendiente del terreno (20%) 9. Visualización alternativa de la pendiente del terreno (15%) 10. Visualización de la orientación del terreno (15%) Bibliografía Reproducibilidad
#1. Introducción Este cuaderno ilustra cómo calcular los atributos geomorfométricos del terreno a partir de modelos digitales de elevación cuadriculados. Su objetivo es enseñar conceptos y herramientas geoinformáticas básicas a estudiantes de pregrado de agronomía de la Universidad Nacional de Colombia.
Para calcular los atributos del terreno, utilizaremos el paquete MultiscaleDTM de R.
Al crear su cuaderno, asegúrese de agregar el texto necesario para: (i) explicar cada fragmento de código; y (ii) describir el resultado correspondiente. Estos dos criterios se tendrán en cuenta para evaluar la calidad de su cuaderno.
#2. Configuración Los paquetes necesarios deben instalarse previamente desde la consola. #SETUP #INSTALL THIS PACKAGE FROM THE CONSOLE, NOT FROM THIS CHUNK paquetes = c(“MultiscaleDTM”, “exactextractr”) install.packages(paquetes)
Empecemos por limpiar la memoria:
rm(list=ls())Ahora, carguemos las bibliotecas necesarias:
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.10.2, PROJ 9.5.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.8.54library(MultiscaleDTM)
library(exactextractr)Al ejecutar un fragmento, preste atención a cualquier mensaje que indique un error o advierta sobre posibles conflictos.
Los conflictos pueden deberse a funciones de diferentes paquetes con el mismo nombre (por ejemplo, terra::extract y tidyr::extract). Para evitar estos problemas, es necesario llamar a varias funciones utilizando la ruta larga, es decir, package_name::function.
#3. Introducción a MultiscaleDTM Este paquete R calcula atributos geomorfométricos del terreno multiescala a partir de modelos digitales del terreno (MDT; es decir, rásteres de elevación o batimetría) con cuadrícula regular a través de una ventana específica.
En caso de que necesites refrescar la forma de calcular atributos geomorfométricos, revisa las diapositivas del curso.
#4. Lectura de datos de entrada Primero, cargaremos un DEM ráster para nuestro departamento. En este cuaderno, usaré un DEM obtenido previamente con este cuaderno.
Leamos el DEM con terra:
(dem = terra::rast("C:/Users/menju/OneDrive/Documentos/GB2/Rstudio2/elev_ELEVACION_z10.tif"))## class : SpatRaster
## size : 4092, 3078, 1 (nrow, ncol, nlyr)
## resolution : 0.000685414, 0.000685414 (x, y)
## extent : -76.28906, -74.17936, 2.811272, 5.615986 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elev_ELEVACION_z10.tif
## name : elev_ELEVACION_z10
## min value : -427
## max value : 5372Reduzcamos la resolución de DEM para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Ahora, leeremos nuestros municipios utilizando la biblioteca sf:
(munic <- sf::st_read("C:/Users/menju/OneDrive/Documentos/GB2/P1/DATOS/MUNI_TOL.gpkg"))## Reading layer `MUNICIPIOS_TOL' from data source
## `C:\Users\menju\OneDrive\Documentos\GB2\P1\DATOS\MUNI_TOL.gpkg'
## using driver `GPKG'
## Simple feature collection with 47 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS: MAGNA-SIRGASTenga en cuenta que ambos conjuntos de datos están en coordenadas geográficas.
Ahora, recortemos el objeto dem2 a los límites de nuestro departamento:
dem3 = terra::crop(dem2,munic, mask=TRUE)## Warning: [crop] CRS do not match#5. Transformación de coordenadas (10%) Para calcular los atributos del terreno, el MDE debe estar en coordenadas planas. Como recordará, el IGAC cambió en 2020 a un sistema de coordenadas planas con origen único nacional. ¿Cuál es su código EPSG?
Modifique los siguientes dos fragmentos de código para transformar los datos de entrada al sistema de coordenadas nacional de Colombia.
# 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 : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_ELEVACION_z10
## min value : -1.783039
## max value : 5277.944824También necesitamos realizar una transformación similar para el objeto vectorial:
# 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. Cálculo de atributos del terreno (20%) Veamos cómo calcular la pendiente y la orientación.
Vaya a la consola y escriba ?SlpAsp. A continuación, en la ventana de ayuda a la derecha, se mostrará la información correspondiente:
(Ver Consola)
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
## size : 1787, 1197, 2 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0000, 7.917609e-04
## max values : 69.1466, 3.599997e+02Ahora, dividamos el objeto slp_asp en dos partes.
Capa 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 69.1466¿Cuál es la distribución de valores de pendiente?
terra::hist(slope,
main = "Tolima's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
Capa 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 7.917609e-04
## max value : 3.599997e+02terra::hist(aspect,
main = "Tolima aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
Además, convirtamos los grados de pendiente en porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 1787, 1197, 1 (nrow, ncol, nlyr)
## resolution : 151.7803, 151.7803 (x, y)
## extent : 4654850, 4836531, 1875472, 2146704 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 262.5146terra::hist(slope_perc,
main = "Tolima's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
#7. Cálculo de estadísticas zonales Una operación de estadísticas zonales calcula estadísticas sobre los valores de celda de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En nuestro caso, nos interesa calcular el promedio de los valores de pendiente por municipio.
Primero, conozcamos la clasificación de pendientes del IGAC para fines agrológicos:
(Observar Tabla)
Si quieres entender la lógica de esta clasificación, consulta este enlace.
Ahora, manos a la obra: (Reclacificar un Ráster)
# 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)Ahora, calculemos las estadísticas zonales usando exactextractr:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |= | 2% | |=== | 4% | |==== | 6% | |====== | 9% | |======= | 11% | |========= | 13% | |========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================== | 26% | |=================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |=============================== | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |======================================= | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |=================================================== | 72% | |==================================================== | 74% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 85% | |============================================================= | 87% | |=============================================================== | 89% | |================================================================ | 91% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%## [1] 28.327206 22.369404 15.429619 3.669801 33.040310 9.219158 27.107445
## [8] 36.053734 15.664359 35.969097 28.613630 18.607035 7.297594 18.536375
## [15] 20.129711 1.252090 18.143921 1.507208 24.973742 2.159666 34.832108
## [22] 16.668184 22.115511 11.497709 32.145004 11.752955 16.978456 27.755978
## [29] 16.803055 18.587324 20.755070 8.394407 34.529945 15.837931 10.465549
## [36] 41.639145 27.062590 30.451038 1.194180 36.182438 11.381176 31.772017
## [43] 13.304363 17.442604 11.510288 35.072704 18.109522## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)hist(munic$mean_slope,
main = "Tolima'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% | |= | 2% | |=== | 4% | |==== | 6% | |====== | 9% | |======= | 11% | |========= | 13% | |========== | 15% | |============ | 17% | |============= | 19% | |=============== | 21% | |================ | 23% | |================== | 26% | |=================== | 28% | |===================== | 30% | |====================== | 32% | |======================== | 34% | |========================= | 36% | |=========================== | 38% | |============================ | 40% | |============================== | 43% | |=============================== | 45% | |================================= | 47% | |================================== | 49% | |==================================== | 51% | |===================================== | 53% | |======================================= | 55% | |======================================== | 57% | |========================================== | 60% | |=========================================== | 62% | |============================================= | 64% | |============================================== | 66% | |================================================ | 68% | |================================================= | 70% | |=================================================== | 72% | |==================================================== | 74% | |====================================================== | 77% | |======================================================= | 79% | |========================================================= | 81% | |========================================================== | 83% | |============================================================ | 85% | |============================================================= | 87% | |=============================================================== | 89% | |================================================================ | 91% | |================================================================== | 94% | |=================================================================== | 96% | |===================================================================== | 98% | |======================================================================| 100%## [1] 5 4 1 1 5 1 5 5 1 5 5 4 1 4 4 1 4 1 5 1 5 4 4 1 5 4 4 5 1 5 4 1 5 4 1 5 4 5
## [39] 1 5 1 5 1 4 1 5 4## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)hist(munic$class,
main = "Tolima's reclassified slope",
xlab = "Slope (as a category)")
Transformemos nuevamente la pendiente a coordenadas geográficas:
# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 1793, 1199, 1 (nrow, ncol, nlyr)
## resolution : 0.001370805, 0.001370805 (x, y)
## extent : -76.11484, -74.47124, 2.868096, 5.325948 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.2429543
## max value : 259.9135742#This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 1793, 1199, 1 (nrow, ncol, nlyr)
## resolution : 0.001370805, 0.001370805 (x, y)
## extent : -76.11484, -74.47124, 2.868096, 5.325948 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1.254625e-04
## max value : 2.599136e+02#8. Trazado de la pendiente del terreno (20%) Primero, necesitaremos una paleta de colores para la pendiente:
# 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")Ahora es el momento de trazar el gráfico. Asegúrate de cambiar varias opciones en el código si es necesario:
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 Tolima (%)")## 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'