MOVIMIENTOS DE LADERA: FACTORES NATURALES
Asignatura G171: Procesos Geomorfológicos
Carolina Garmendia y Domingo F. Rasilla
2025-04-11
INTRODUCCIÓN
En primer lugar, se inactiva la notación científica.
options(scipen = 999)Para la extracción de las variables naturales que serán evaluadas, es necesario contar con los siguientes ficheros:
library(sf)## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUEsf_region <- st_read("D:/G171_Procesos_Geomorfologicos_2022/2_Procesos_Ladera/Practica_Laderas/sf_region.geojson")## Reading layer `sf_region' from data source
## `D:\G171_Procesos_Geomorfologicos_2022\2_Procesos_Ladera\Practica_Laderas\sf_region.geojson'
## using driver `GeoJSON'
## Simple feature collection with 93 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -126163 ymin: 5366748 xmax: 938906 ymax: 6173275
## Projected CRS: WGS 84 / UTM zone 10Nbbox_new <- st_read("D:/G171_Procesos_Geomorfologicos_2022/2_Procesos_Ladera/Practica_Laderas/bbox_new.geojson")## Reading layer `bbox_new' from data source
## `D:\G171_Procesos_Geomorfologicos_2022\2_Procesos_Ladera\Practica_Laderas\bbox_new.geojson'
## using driver `GeoJSON'
## Simple feature collection with 1 feature and 0 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -864334 ymin: 5062152 xmax: 1616492 ymax: 6703487
## Projected CRS: WGS 84 / UTM zone 10Nplot(st_geometry(sf_region))
plot(st_geometry(bbox_new), add = TRUE)
EXTRACCIÓN DE VARIABLES TOPOGRÁFICAS.
Importación de un modelo digital del terreno.
Para analizar las condiciones topográficas favorables a la génesis de
los movimientos de ladera, se utilizará la información contenida en un
modelo digital del terreno. Para ello se utilizará el paquete
elevatr, a través del que se puede
descargar directamente un mdt con diferentes niveles de resolución.
Este paquete puede descargarse directamente desde desde el repositorio CRAN.
library(terra)## Warning: package 'terra' was built under R version 4.3.3## terra 1.8.21#install.packages("elevatr")
library(elevatr)## Warning: package 'elevatr' was built under R version 4.3.3## 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.La función elevatr::get_elev_raster()descarga un mdt en
formato ráster desde el sitio web “AWS Open Data Terrain Tiles”. Para
hacerlo necesitamos varios argumentos:
locations es el marco espacial del ámbito analizado.
z el nivel de resolución del mdt veáse `
clip objeto espacial para recortar
mdt <- elevatr::get_elev_raster(locations = bbox_new, # El área espacial que nos interesa
z = 7, # La resolución del mdt
clip = "locations") # El ámbito espacial que nos interesa## Mosaicing & Projecting## Clipping DEM to locations## Note: Elevation units are in meters.Como el modelo digital del terreno es suministrado en un formato
diferente al SpatRaster de terra(), se
transforma con el siguiente código.
mdt <- as(mdt, "SpatRaster") Si existen valores inferiores a 0 (topografía submarina) deben eliminarse.
mdt <- app(mdt,
function(x){x[x < 0] <- NA;return(x)})
names(mdt) <- "altitud"Reproyección del objeto espacial
A continuación, se procede a reproyectar el objeto raster
mdt, que originalmente fue descargado en forma de
coordenadas geográficas, a coordenadas proyectadas correspondientes a la
zona 10N UTM.
mdt <- project(mdt,
"EPSG:32610")
crs(mdt, describe = TRUE)## name authority code
## 1 WGS 84 / UTM zone 10N EPSG 32610
## area
## 1 Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)
## extent
## 1 -126, -120, 0, 84Cartografiar el mapa
Para elaborar un mapa que muestre los movimientos de ladera
superpuestos a un mdt, usando ggplot(), es necesario
convertir el objeto ráster mdt en un dataframe, mediante la
función as.data.frame(), especificando las dimensiones x e
y. En este dataframe, el nombre de la variable será
altitud.
mdt_df <- as.data.frame(mdt,
xy = TRUE)
colnames(mdt_df)[3] <- "altitud"Además, las filas con valores ausentes (“NA”) serán eliminadas,
usando la función complete.cases().
mdt_df <- mdt_df[complete.cases(mdt_df), ]
mdt_df <- subset(mdt_df,
altitud >= 0)A continuación ya se pueden representar gráficamente los movimientos de ladera con un mapa topográfico de fondo.
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3ggplot() +
geom_raster(data = mdt_df, aes(x = x, y = y, fill = altitud)) +
geom_sf(data = sf_region, color = "red", fill = NA) +
coord_sf() +
scale_fill_viridis_c(limits = c(0, 4000)) +
labs(title = "Mapa topográfico de Columbia Britanica", x = "Longitud", y = "Latitud", fill = "Altitud (metros)")
Tras la elaboración del mapa, el dataframe mdt_df puede
ser eliminado del Global Environment.
rm(mdt_df)Obtención de parámetros topográficos derivados de un DEM
Con el DEM listo, se pueden derivar otras variables relacionadas con
la altitud, como la pendiente, la rugosidad etc… Para ello se recurre a
la función terrain(), que facilita su cálculo. Los
argumentos que acompañan a la función son muy sencillos:
mdt nombre del objeto ráster sobre el que se aplica la función.
opt nombre del parámetro topográfico a calcular (slope -pendiente-, aspect -orientación-, TPI, TRI, roughness, flowdir).
unit unidad en la que se mide ese parámetro (sólo importante para el cálculo de la pendiente y de la orientación).
Las variables a obtener son:
- La pendiente
mdt_pendiente <- terrain(mdt, v="slope", unit = "degrees")
plot(mdt_pendiente,
col = rainbow(20), main= "Pendiente (º)")
- La orientación de las vertientes
mdt_orientacion <- terrain(mdt, v = "aspect", unit = "degrees")
plot(mdt_orientacion,
col = rainbow(20),
main= "Orientación (º)")
- El TRI (Índice de rugosidad del terreno) es la media de las diferencias absolutas entre el valor de una celda y el valor de las 8 celdas que la rodean (la cantidad de diferencia de elevación entre las celdas adyacentes de un DEM). Proporciona una medida relativa de los cambios de elevación entre una celda de cuadrícula específica y los vecinos.
mdt_tri <- terrain(mdt, v = "TRI")
plot(mdt_tri,
main= "Índice de rugosidad del terreno (TRI)",
col=rainbow(50))
- El TPI (Índice de posición topográfica) es la diferencia entre el valor de una celda y el valor medio de las 8 celdas que la rodean. Los valores positivos de TPI representan ubicaciones que son más altas que el promedio de sus alrededores, según lo definido por la vecindad (crestas). Los valores negativos de TPI representan ubicaciones que están más bajas que sus alrededores (valles); los valores de TPI cercanos a cero son áreas planas (donde la pendiente es cercana a cero) o áreas de pendiente constante.
mdt_tpi <- terrain(mdt, v = "TPI")
plot(mdt_tpi,
main= "Índice de Posición Topográfica (TPI)",
col=rainbow(50))
- La curvatura es el cambio en el grado de pendiente para una distancia dada. Cualquier valor positivo indica una pendiente convexa (abultada hacia arriba), mientras que los valores negativos corresponden a una curva convexa (hacia abajo). Cualquier píxel observado que muestre una curvatura positiva indica la posibilidad de un flujo se disperse lejos de un área central, mientras que los valores negativos indican acumulación. En R, la curvatura se puede calcular mediante el paquete spatialEco, que deberá cargarse a continuación
#install.package(spatialEco)
library(spatialEco) #Para el cálculo de la curvatura## Warning: package 'spatialEco' was built under R version 4.3.2##
## Attaching package: 'spatialEco'## The following objects are masked from 'package:terra':
##
## is.empty, shift, sieveLa curvatura se puede calcular siguiendo varios métodos; el elegido
en este caso es el método de McNab en el argumento
type = "mcnab" de la función
spatialEco::curvature.
mdt_crv <- curvature(mdt, type="mcnab") ## |---------|---------|---------|---------|========================================= plot(mdt_crv,
main= "Curvatura",
col=rainbow(50))
A continuación, se agruparán todos los ficheros ráster relacionados con la topografía en un único objeto ráster que contendrá 7 capas, una por cada variable.
topografia <- c(mdt, mdt_pendiente, mdt_orientacion, mdt_tpi, mdt_tri, mdt_crv)
names(topografia) <- c("altitud", "pendiente", "orientacion", "TPI", "TRI", "curvatura")Este objeto ráster deberá ser guardado
writeRaster(topografia, "D:/G171_Procesos_Geomorfologicos_2022/2_Procesos_Ladera/Practica_Laderas/topografia.tif", overwrite = TRUE)## |---------|---------|---------|---------|========================================= Como ya no son necesarios los objetos ráster originales, serán eliminados del Global environment.
rm(mdt, mdt_pendiente, mdt_orientacion, mdt_tpi, mdt_tri, mdt_crv)Generación de zonas de influencia (“buffer”) alrededor de cada movimiento de ladera
A continuaciónes se obtendrán los valores de esos parámetros en los puntos de muestreo (la localización de los movimientos de ladera), creándose un dataframe para análisis posteriores. Una forma de hacerlo es crear un búfer alrededor de los puntos de muestreo y resumir los valores en el búfer, usando la media, la moda, el máximo, etc.
buffer <- st_buffer(sf_region,
500) # Tamaño del buffer: 500m
class(buffer) # Tipo de dato devuelto## [1] "sf" "data.frame"Podemos verificar que se ha realizado correctamente superponiendo al MDT la localización de los movimientos y los buffer circundantes.
plot(topografia[[1]],col = terrain.colors(100))
plot(st_geometry(sf_region), add = TRUE)
plot(st_geometry(buffer), add = TRUE)
Por último, se creará un dataframe con los valores medios de las
variables topográficas dentro del área de buffer, mediante la función
terra::extract().
df_topografia <- extract(topografia,
buffer,
fun = mean)Como la orientación es una variable circular, el valor 0º es igual al de 360º. Por esta razón, es conveniente transformar los valores en grados en cuadrantes, que se convierten en categorías que resumen toda la gama de orientaciones en sólo 4 direcciones principales, siguiendo el siguiente cuadro:
| Orientación en grados | Dirección |
|---|---|
| 0 a 45 y 315 a 360 | Norte |
| 45 a 135 | Este |
| 135 a 225 | Sur |
| 225 a 315 | Oeste |
Para ello podemos utilizar la siguiente función:
cuadrantes_4 <- function(grados) {
categorias <- c("N", "E", "S", "W") # "N" se repite para incluir 360°
cortes <- seq(-45, 360, by = 90) # Definir los rangos para cada dirección
# Convertir los grados en categorías
factor_categorico <- cut(grados, breaks = cortes, labels = categorias, include.lowest = TRUE)
return(factor_categorico)
}
names(df_topografia)[4] <- "borrame"
df_topografia$orientacion <- cuadrantes_4(df_topografia$borrame)
df_topografia <- subset(df_topografia, select =-borrame) EXTRACCIÓN DE VARIABLES CLIMÁTICAS: LA PRECIPITACIÓN ANUAL.
La base de datos WorldClim proporciona datos de diferentes variables climáticas con varias resoluciones espaciales. Se puede extraer directamente desde internet con el paquete **geodata*. A continuación se extraerá la precipitación media anual.
library(geodata)## Warning: package 'geodata' was built under R version 4.3.3precipitacion <- worldclim_global(var="prec",
res = 5,
path=tempdir())
precipitacion <- project(precipitacion,
"EPSG:32610")## |---------|---------|---------|---------|========================================= crs(precipitacion, describe = TRUE)## name authority code
## 1 WGS 84 / UTM zone 10N EPSG 32610
## area
## 1 Between 126°W and 120°W, northern hemisphere between equator and 84°N, onshore and offshore. Canada - British Columbia (BC); Northwest Territories (NWT); Nunavut; Yukon. United States (USA) - Alaska (AK)
## extent
## 1 -126, -120, 0, 84Recorte y enmascaramiento:
precipitacion_recortado <- crop(precipitacion,
bbox_new)
precipitacion_enmascarado <- mask(precipitacion_recortado,
bbox_new)
precipitacion_mensual <- project(precipitacion_enmascarado,
topografia)## |---------|---------|---------|---------|========================================= El fichero original contiene varias capas, cada una de ellas
correspondiente a un mes del año. Para calcular el total anual de
precipitaciones debe aplicarse la función app() del paquete
terra, a través de la cual se puede obtener la
precipitación anual.
precipitacion_anual <- terra::app(precipitacion_mensual, sum)## |---------|---------|---------|---------|========================================= A continuación se representa el mapa de la precipitación anual.
plot(precipitacion_anual,
main= "Precipitación anual",
col=rainbow(50))
Este objeto ráster podrá ser guardado
writeRaster(precipitacion_anual, "D:/G171_Procesos_Geomorfologicos_2022/2_Procesos_Ladera/Practica_Laderas/precipitacion_anual.tif", overwrite = TRUE)Se utiliza terra::extract() para extraer los datos de
lluvia anual correspondientes a la zona en la que ocurrión cada
movimiento de ladera.
df_precipitacion <- extract(precipitacion_anual, sf_region)
names(df_precipitacion) <- c("ID", "Precipitacion")CREACIÓN DE UN DATAFRAME CON LOS VALORES TOPOGRÁFICOS Y DE PRECIPITACIÓN.
Por último, todos los datos extraidos de la topografía y de la
precipitación se unirán en un único dataframe que llamaremos
df_naturales.
library(plyr)## Warning: package 'plyr' was built under R version 4.3.2df_naturales <- join_all(list(df_topografia, df_precipitacion),
by='ID',
type='left')
View(df_naturales)A continuación se grabará este dataframe como fichero
.csv.
write.csv2(df_naturales,
"D:/G171_Procesos_Geomorfologicos_2022/2_Procesos_Ladera/Practica_Laderas/df_naturales.csv",
row.names = FALSE) Finalmente, se elimina todo el conjunto de ficheros innecesarios.
rm(list = ls())