Atributos geomorfométricos del terreno del Caquetá
- 1. Introducción
- 2. Instalación de programas
- 3. Introducción al DTM
- 4. Lectura de datos inicial
- 5. Transformar coordenadas
- 6. Cálculo de atributos del terreno
- 7. Cálculo de estadísticas zonales
- 8. Trazado de la pendiente del terreno
- 9. Visualización alternativa de la pendiente del terreno
- 10. Visualización de la orientación del terreno
- 11. Bibliografía
1. Introducción
Este Notebook se realiza como forma de presentación del segundo examen del curso de Geomática básica de la Universidad Nacional. A continuación se expondran los respectivos códigos para obtener visualizaciones de atributos geomorfométricos, tales como la elevación y la orientación, a partir de modelos dígitales de elevación, específicamente en el departamento del Caquetá.
Se realizara una descripción detallada de cada proceso para lograr un buen entendimiento de los códigos a utilizar, junto con la teoría necesaria para interpretar adecuadamente las distintas visualizaciones que arroje el programa.
2. Instalación de programas
Para la obtención de los atributos que deseamos visualizar vamos a utlizar los siguientes paquetes: - Elevatr: Este paquete nos permite descargar datos de elevación dígital (DEM), muy útil para estudios topográficos. - Sf: Funciona para cargar, manipular y visualizar shapefiles, convertir coordenadas, y hacer análisis geoespacial. Trabaja con datos vectoriales (puntos, lineas, polígonos). - Leaflet: Permite crear mapas interactivos con marcadores y datos espaciales. - Terra: Este paquete procesa datos raster y vectores de forma eficiente, Ayuda mucho a la hora de realizar, por ejemplo, análisis de relieve. - MultiscaleDTM: Extrae métricas avanzadas de un DEM para realizar un análisis más detallado de un terreno. - Exactextractr: Extrae valores de un raster de forma eficiente y precisa.
Para descargar los paquetes descritos anteriormente escribiremos el siguiente código directamente en la consola: paquetes = c(“MultiscaleDTM”, “exactextractr”, “elevatr”, “sf”, “leaflet”, “terra”, “MultiscaleDTM”, “exactextractr”) install.packages(paquetes)
A continuación procederemos a leer todos los paquetes y a limpiar la memoria, ahora sí desde un chunk
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.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.8.10library(MultiscaleDTM)
library(exactextractr)rm(list=ls())Es posible que al momento de la descarga de los programas o la lectura aparezcan mensajes de Aviso referentes a las versiones de los distintos programas
3. Introducción al DTM
Vamos a utilizar la función SlpAsp() proporcionada por el paquete MultiscaleDTM, esta función se utiliza para calcular dos métricas geomorfológicas importantes a partir de un Modelo Digital de Elevación (DEM): La elevación y la orientación
- Pendiente (Slope): Calcula la pendiente en grados o radianes a partir de un DEM.
- Orientación (Aspect): Calcula la dirección en la que una pendiente se orienta en grados (0-360)
El cálculo tanto de la pendiente como de la orientación se realiza usando una ventana de píxeles de entrada, que tradicionalmente es de 3x3, junto con unas funciones que permite trabajar con diferentes métodos y escalas para obtener una descripción más precisa de las características del terreno. Estas funciones son el método de torre (rook), el método reina (queen) o el método de límite (boundary). El que trabajeremos para obtener los datos es el método reina que trabaja con cuatro celdas (La de arriba, la de abajo, izquierda y derecha de la celda focal) y adicionalmente otras cuatro celdas de las esquinas adicionales a la celda focal.
4. Lectura de datos inicial
Los primeros datos que cargaremos será un DEM raster previamente realizado en otra sesión. Este DEM se encuentra en un archivo .tif y se cargará con ayuda de la biblioteca terra
(dem = terra::rast("C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/elevación_caquetá.tif"))## class : SpatRaster
## dimensions : 6142, 8193, 1 (nrow, ncol, nlyr)
## resolution : 0.0006865479, 0.0006865479 (x, y)
## extent : -76.64062, -71.01574, -1.054321, 3.162456 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : elevación_caquetá.tif
## name : elevación_caquetá
## min value : -1010
## max value : 5372Usaremos a continuación un código que nos permitira reducir la resolución del DEM con el fin de ahorrar memoria
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Ahora subiremos los datos del departamento, en este caso los del Caquetá, con ayuda de la biblioteca sf. Este archivo es un geopack que previamente debe estar descargado y localizado en la misma carpeta de archivos desde donde se creo este Notebook.
(munic <- sf::st_read("C:/Users/Karol Ruiz/Documents/GEOMATICA 2024-2/PROYECTO4/DATOS/mun_caquetá.gpkg"))## Multiple layers are present in data source C:\Users\Karol Ruiz\Documents\GEOMATICA 2024-2\PROYECTO4\DATOS\mun_caquetá.gpkg, reading layer `municipios_colombia'.
## Use `st_layers' to list all layer names and their type in a data source.
## Set the `layer' argument in `st_read' to read a particular layer.## Warning in CPL_read_ogr(dsn, layer, query, as.character(options), quiet, :
## automatically selected the first layer in a data source containing more than
## one.## Reading layer `municipios_colombia' from data source
## `C:\Users\Karol Ruiz\Documents\GEOMATICA 2024-2\PROYECTO4\DATOS\mun_caquetá.gpkg'
## using driver `GPKG'
## Simple feature collection with 16 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.938334
## Geodetic CRS: MAGNA-SIRGAS## Simple feature collection with 16 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.30622 ymin: -0.70584 xmax: -71.25385 ymax: 2.938334
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 18 001 18001 CAQUETÁ FLORENCIA
## 2 18 029 18029 CAQUETÁ ALBANIA
## 3 18 094 18094 CAQUETÁ BELÉN DE LOS ANDAQUÍES
## 4 18 150 18150 CAQUETÁ CARTAGENA DEL CHAIRÁ
## 5 18 205 18205 CAQUETÁ CURILLO
## 6 18 247 18247 CAQUETÁ EL DONCELLO
## 7 18 256 18256 CAQUETÁ EL PAUJÍL
## 8 18 410 18410 CAQUETÁ LA MONTAÑITA
## 9 18 460 18460 CAQUETÁ MILÁN
## 10 18 479 18479 CAQUETÁ MORELIA
## geom
## 1 MULTIPOLYGON (((-75.47547 1...
## 2 MULTIPOLYGON (((-75.89531 1...
## 3 MULTIPOLYGON (((-75.78705 1...
## 4 MULTIPOLYGON (((-74.78722 1...
## 5 MULTIPOLYGON (((-76.09066 1...
## 6 MULTIPOLYGON (((-75.36167 2...
## 7 MULTIPOLYGON (((-75.36638 2...
## 8 MULTIPOLYGON (((-75.47712 1...
## 9 MULTIPOLYGON (((-75.32574 1...
## 10 MULTIPOLYGON (((-75.77185 1...Aquí se puede visualizar la tabla de atributos que quedo guardado con el nombre munic. El mensaje de Aviso indica que se selecciono la primera capa, puesto que el archivo a utilizar tenía varias capas.
Ahora recortaremos el DEM a los límites de nuestro departamento mediante el siguiente código:
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Transformar coordenadas
Hay que modificar las coordenadas del DEM para que estén coordenadas planas y así poder trabajar los datos. Con los siguientes dos códigos cambiaremos las coordenadas al sistema de coordenadas Nacional de Colombia. Este primer código es para modificar las coordenadas del DEM
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 2646, 3688, 1 (nrow, ncol, nlyr)
## resolution : 152.4091, 152.4091 (x, y)
## extent : 4632112, 5194197, 1480000, 1883274 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elevación_caquetá
## min value : -135.9608
## max value : 3607.0352El siguiente código es para modificar las coordenadas del objeto Munic que es un vector.
(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 16 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4632126 ymin: 1480062 xmax: 5194256 ymax: 1882858
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 18 001 18001 CAQUETÁ FLORENCIA
## 2 18 029 18029 CAQUETÁ ALBANIA
## 3 18 094 18094 CAQUETÁ BELÉN DE LOS ANDAQUÍES
## 4 18 150 18150 CAQUETÁ CARTAGENA DEL CHAIRÁ
## 5 18 205 18205 CAQUETÁ CURILLO
## 6 18 247 18247 CAQUETÁ EL DONCELLO
## 7 18 256 18256 CAQUETÁ EL PAUJÍL
## 8 18 410 18410 CAQUETÁ LA MONTAÑITA
## 9 18 460 18460 CAQUETÁ MILÁN
## 10 18 479 18479 CAQUETÁ MORELIA
## geom
## 1 MULTIPOLYGON (((4724662 172...
## 2 MULTIPOLYGON (((4677907 170...
## 3 MULTIPOLYGON (((4690016 175...
## 4 MULTIPOLYGON (((4801248 173...
## 5 MULTIPOLYGON (((4656116 167...
## 6 MULTIPOLYGON (((4737450 181...
## 7 MULTIPOLYGON (((4736906 180...
## 8 MULTIPOLYGON (((4724479 172...
## 9 MULTIPOLYGON (((4741294 169...
## 10 MULTIPOLYGON (((4691681 173...6. Cálculo de atributos del terreno
Toda la información necesaria sobre la función SlpAsp la puedes encontrar escribiendo directamente en la consola ?SlpAsp. Se abrirá una ventana en la derecha con toda la información.
Tras esto, usaremos la función con el método queen y unidad en degrees (Grados). Vamos a hallar tanto la elevación como la orientación (metrics) y se eliminaran los los valores NA para evitar problemas de cálculo (na.rm=TRUE),no sé incluira la escala (include_scale=FALSE), pero si se incluirá una máscara a los valores de orientación para evitar errores en zonas planas (mask_aspect=TRUE)
(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 : 2646, 3688, 2 (nrow, ncol, nlyr)
## resolution : 152.4091, 152.4091 (x, y)
## extent : 4632112, 5194197, 1480000, 1883274 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 1.362937e-04, 1.454152e-04
## max values : 5.267276e+01, 3.599999e+02Ahora dividiremos el objeto en 2 capas.
La primera capa será exclusivamente de la elevación. La obtendremos como se muestra a continuación:
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 2646, 3688, 1 (nrow, ncol, nlyr)
## resolution : 152.4091, 152.4091 (x, y)
## extent : 4632112, 5194197, 1480000, 1883274 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 1.362937e-04
## max value : 5.267276e+01Ahora visualizaremos la elevación del departamento mediante un histograma, con ayuda de la función terra:
terra::hist(slope,
main = "Elevación Caquetá",
xlab = "Elevación (grados)")## Warning: [hist] a sample of 10% of the cells was used (of which 59% was NA)
En este histograma se puede evidenciar que el departamento del Caquetá
tiene zonas con bastante elevación que probablemente corresponde a la
zona del departamento donde pasa la cordillera, pero en general el resto
del departamento tiene muy poca elevación.
A continuación dividiremos la segunda capa con la orientación, mediante el siguiente código:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 2646, 3688, 1 (nrow, ncol, nlyr)
## resolution : 152.4091, 152.4091 (x, y)
## extent : 4632112, 5194197, 1480000, 1883274 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.454152e-04
## max value : 3.599999e+02Al igual que la elevación aquí también visualizaremos los datos mediante un histograma
terra::hist(aspect,
main = "Orientación Caquetá",
xlab = "Orientación (grados)")## Warning: [hist] a sample of 10% of the cells was used (of which 59% was NA)
De este histograma podemos interpretar que el departamento del Caquetá
tiene una orientación con una tendencia equitativa o parecida en todos
los grados, es decir, que podemos encontrar pendientes orientadas hacía
todos los lados o direcciones.
Adicionalmente, convertiremos los grados de la elevación a porcentaje, mediante la siguiente ecuación. Recordemos que la relación que nos proporciona la elevación se da en términos de tangente del ángulo de elevación, por lo cual debemos hallar la tangente del ángulo y multiplicarlo por 100 para determinar el porcentaje.
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## dimensions : 2646, 3688, 1 (nrow, ncol, nlyr)
## resolution : 152.4091, 152.4091 (x, y)
## extent : 4632112, 5194197, 1480000, 1883274 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 2.378774e-04
## max value : 1.311394e+02Ahora visualizaremos otro histograma con los datos de la elevación en porcentaje
terra::hist(slope_perc,
main = "Elevación Caquetá",
xlab = "Elevación (porcentaje)")## Warning: [hist] a sample of 10% of the cells was used (of which 59% was NA)
Como podemos ver este histograma es muy parecido al anterior en el que
visualizabamos la elevación en grados, sin embargo, los intervalos en el
eje X sobre el que se mide la variable varian por el cambio de
unidad.
7. Cálculo de estadísticas zonales
Una operación de estadísticas zonales es aquella que calcula estadísticas sobre los valores de celda de un ráster dentro de las zonas definidas por otro dataset. En nuestro caso, nos interesa calcular el promedio de los valores de pendiente por municipio.
En este caso nos guiaremos con la clasificación del IGAC en efectos agrológicos que se pueden consultar en el siguiente link: (“https://www.igac.gov.co/sites/default/files/listadomaestro/in-agr-pc02-03_elaboracion_de_cartografia_geomorfologica_0.pdf”)
A continuación vamos a reclasificar los valores para hallar el valor de la pendeinte por municipio
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 vamos a calcular la estadística zonal con la función exactextractr:
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |==== | 6% | |========= | 12% | |============= | 19% | |================== | 25% | |====================== | 31% | |========================== | 38% | |=============================== | 44% | |=================================== | 50% | |======================================= | 56% | |============================================ | 62% | |================================================ | 69% | |==================================================== | 75% | |========================================================= | 81% | |============================================================= | 88% | |================================================================== | 94% | |======================================================================| 100%## [1] 21.042652 1.809873 20.804121 2.110796 1.994576 13.908103 10.098153
## [8] 5.944215 2.320801 4.048290 16.770296 21.260427 6.011272 2.827471
## [15] 2.246294 2.270885Ahora con estos valores crearemos otro histograma de la elevación del departamento donde los intervalos serán de 5 en 5 desde 0 hasta 25
hist(munic$mean_slope,
main = "Pendiente media del Caquetá",
xlab = "Media de elevación (porcentaje)")
Ahora reclasificamos otra vez y hacemos otro histograma pero esta vez
con los intervalos de 1 en 1 desde 1 hasta 5
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |==== | 6% | |========= | 12% | |============= | 19% | |================== | 25% | |====================== | 31% | |========================== | 38% | |=============================== | 44% | |=================================== | 50% | |======================================= | 56% | |============================================ | 62% | |================================================ | 69% | |==================================================== | 75% | |========================================================= | 81% | |============================================================= | 88% | |================================================================== | 94% | |======================================================================| 100%## [1] 5 1 5 1 1 1 1 1 1 1 1 5 1 1 1 1hist(munic$class,
main = "Caqueta's reclassified slope",
xlab = "Slope (as a category)")
Ahora transformamos la pendiente de nuevo en coordenadas geográficas
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 2658, 3684, 1 (nrow, ncol, nlyr)
## resolution : 0.001373092, 0.001373092 (x, y)
## extent : -76.30992, -71.25145, -0.7061153, 2.943562 (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
## dimensions : 2658, 3684, 1 (nrow, ncol, nlyr)
## resolution : 0.001373092, 0.001373092 (x, y)
## extent : -76.30992, -71.25145, -0.7061153, 2.943562 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.01579605
## max value : 125.635887158. Trazado de la pendiente del terreno
Ahora realizaremos un mapa interactivo para visualizar la elevación, pero en primer lugar elegiremos una paleta de colores adecuada. Podemos agregar o quitar colores para mejorar la visualización de la elevación
palredgreen <- colorNumeric(c("darkseagreen3", "green", "yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")Ahora si vamos a desarrollar el mapa. Inicialmente añadiremos un línea de código que nos va a ayudar a reducir la resolución del mapa para que no sea tan pesado, después iremos poniendo todas las específicaciones de valores, la paleta de colores, los datos que deseamos visualizar, los marcadores y el título.
slope.geo <- aggregate(slope.geo, fact = 2, fun = mean) ## |---------|---------|---------|---------|========================================= 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 = "Elevación del terreno Caquetá (%)")## 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'