Elevasion en R, aplicado al departamento del Tolima
Myriam Yusely Fonseca Salguero
2025-02-02
Instalar y cargar las librerias
# install.packages("MultiscaleDTM")
# install.packages("MultiscaleDTM")
library(MultiscaleDTM)## Cargando paquete requerido: terra## terra 1.8.15library(exactextractr)
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)
library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:terra':
##
## intersect, union## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncargar archivo de polígonos
munic = st_read("Division_muncipios_completos.shp")## Reading layer `Division_muncipios_completos' from data source
## `C:\Users\yusel\Downloads\tar1 (2)\tar1\Division_muncipios_completos.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1122 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -81.73562 ymin: -4.225938 xmax: -66.85689 ymax: 13.39473
## Geodetic CRS: WGS 84TOLIMA_SHP = munic %>%
filter(Depto == "Tolima")3.Introduction to MultiscaleDTM
Cargar el archivo .tiff del TOLIMA siguiendo la sesión anterior usando terra
dem = terra::rast("TOLIMA_ELEVATION.tif")
dem## class : SpatRaster
## dimensions : 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. : lon/lat WGS 84 (EPSG:4326)
## source : TOLIMA_ELEVATION.tif
## name : TOLIMA_ELEVATION
## min value : -427
## max value : 5372agregar por la media de los pixeles los modelos de elevación reduciendo la resolución.
dem2 = terra::aggregate(dem, 2, "mean")## |---------|---------|---------|---------|========================================= Cargar mi capa de polígonos de los municipios del TOLIMA usando la libreria sf
(munic <- sf::st_read("TOLIMA_MAGNA_SIRGAS.shp"))## Reading layer `TOLIMA_MAGNA_SIRGAS' from data source
## `C:\Users\yusel\Downloads\tar1 (2)\tar1\TOLIMA_MAGNA_SIRGAS.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 47 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS: WGS 84## Simple feature collection with 47 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.10574 ymin: 2.871081 xmax: -74.47482 ymax: 5.319342
## Geodetic CRS: WGS 84
## First 10 features:
## MpCodigo MpNombre Depto Cod_Mun geometry
## 1 73001 Ibagué Tolima 73001 POLYGON ((-75.38684 4.69905...
## 2 73024 Alpujarra Tolima 73024 POLYGON ((-74.99537 3.51211...
## 3 73026 Alvarado Tolima 73026 POLYGON ((-74.98068 4.70338...
## 4 73030 Ambalema Tolima 73030 POLYGON ((-74.74744 4.94928...
## 5 73043 Anzoátegui Tolima 73043 POLYGON ((-75.20635 4.71674...
## 6 73055 Armero Tolima 73055 POLYGON ((-74.88728 5.14702...
## 7 73067 Ataco Tolima 73067 POLYGON ((-75.32703 3.66528...
## 8 73124 Cajamarca Tolima 73124 POLYGON ((-75.48451 4.58181...
## 9 73148 Carmen De Apicalá Tolima 73148 POLYGON ((-74.74449 4.24728...
## 10 73152 Casabianca Tolima 73152 POLYGON ((-75.07656 5.12370...Cortar mi modelo de elevación a mi capa de polígonos en el TOLIMA para que me quede solo la información que esta dentro de esa geometría. Y lo que esté por fuera lo deja como NA.
dem3 = terra::crop(dem2,munic, mask=TRUE)5. Transforming coordinates
El EPSG es un identificador único para sistemas de coordenadas y proyecciones geográficas. Cada código representa un Sistema de Referencia de Coordenadas (CRS) con el correspondiente datum, proyección y unidad de medida. Es muy importante para que las operaciones entre las capas se de correctamente. Si están en diferente sistema de coordenadas no se pueden operar.
Entonces, reproyectamos al sistema de coordenadas planas con origen único nacional
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 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 : TOLIMA_ELEVATION
## min value : -1.783039
## max value : 5277.944824(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 47 features and 4 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 4654882 ymin: 1875660 xmax: 4836359 ymax: 2146078
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## MpCodigo MpNombre Depto Cod_Mun geometry
## 1 73001 Ibagué Tolima 73001 POLYGON ((4735321 2077691, ...
## 2 73024 Alpujarra Tolima 73024 POLYGON ((4778423 1946330, ...
## 3 73026 Alvarado Tolima 73026 POLYGON ((4780381 2078030, ...
## 4 73030 Ambalema Tolima 73030 POLYGON ((4806322 2105144, ...
## 5 73043 Anzoátegui Tolima 73043 POLYGON ((4755352 2079582, ...
## 6 73055 Armero Tolima 73055 POLYGON ((4790881 2127047, ...
## 7 73067 Ataco Tolima 73067 POLYGON ((4741618 1963353, ...
## 8 73124 Cajamarca Tolima 73124 POLYGON ((4724438 2064764, ...
## 9 73148 Carmen De Apicalá Tolima 73148 POLYGON ((4806459 2027541, ...
## 10 73152 Casabianca Tolima 73152 POLYGON ((4769890 2124534, ...6. Computing terrain atributes
- Multiscale Slope and Aspect
¿Qué es w? w es el tamaño de la ventana de análisis, en este caso es en filas y columnas (3,3) es decir que la pendiente y el aspecto se calculan usando una matriz de vecinos de 3x3 pixeles
¿Qué es method? method es la forma en como se eligen los pixeles vecinos para el cálculo, en el caso de la opción “queen” indica que tomará los 8 vecinos cercanos al rededor de la celda. Puede ser cambiado a “rook” que serian los vecinos en cruz
Calculamos la pendiente y el aspecto a partir del DEM
(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 : 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+02Tomar solo la capa de pendiente
#capa1
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 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.1466terra::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)
El histograma muestra una gran proporcion que se acerca a 0 es decir son terrenos planos aunque va aumentando a 40 esto nos indica que estos terrenos tendran inclinaciones pronunciadas , esto indica que TOLIMA es mayormente plano.
Tomar solo la capa aspecto
#capa2
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 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's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 47% of the cells was used (of which 51% was NA)
Pendiente en porcentajes
slope_perc = tan(slope*(pi/180))*100
slope_perc## class : SpatRaster
## dimensions : 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. Calcular estadísticas zonales
Creamos la clasificación del IGAC para las pendientes. y clasificamos la pendiente en porcentaje con el valor superior incluido en cada categoria
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 calculamos la media de cada grupo clasificado
(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.327534 22.369406 15.429620 3.670041 33.040310 9.219160 27.107374
## [8] 36.053688 15.664418 35.969097 28.613638 18.602770 7.297594 18.538206
## [15] 20.129709 1.252070 18.143984 1.509809 24.971958 2.159533 34.834579
## [22] 16.666601 22.115627 11.497708 32.145004 11.752780 16.980551 27.755978
## [29] 16.803055 18.586514 20.754772 8.394489 34.529865 15.837936 10.464499
## [36] 41.639069 27.062588 30.451162 1.194180 36.180962 11.381084 31.771980
## [43] 13.298852 17.442602 11.510077 35.072704 18.109520Histograma para ver la distribución de esas medias
hist(munic$mean_slope,
main = "TOLIMA's mean slope",
xlab = "Slope (in percentage)", breaks = 18)
calcular la moda por municipio
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%munic$class## [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 4distribución de las modas
hist(munic$class,
main = "TOLIMA's reclassified slope",
xlab = "Slope (as a category)", breaks = 12)
pendiente en coordenadas geográficas
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 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(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 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+028. Plotting terrain slope
Ajustar la paleta de colores con los valores de la pendiente reproyectada y que los NA sean trasnparentes
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2","darkred"), values(slope.geo), na.color = "transparent", reverse = T)reproyectar capa de municipios para que se puedan visualizar en el mismo sistema de coordenadas
munic <- st_transform(munic, crs = 4326)graficar los polígonos y la pendiente con una letenda de los valores
leaflet(munic) %>%
addTiles() %>%
setView(-75.23421, 5.615986, 8) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$MpNombre, "<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Terrain slope in TOLIMA (%)")los popups tienen el nombre del municipio y la clasificiación de la pendiente asociada al IGAC, se observa la mayor parte del departamento en la clase 1 ligeramente plano, pendiente de 0 a 3%
9.Alternative visualization of terrain slope
munic$slope <- exactextractr::exact_extract(slope, 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%munic$slope## [1] 15.2772141 12.4048328 8.5452795 2.0956366 17.8026676 5.1712618
## [7] 14.8903093 19.4962254 8.6316547 19.3072281 15.4923220 10.3422184
## [13] 4.1039290 10.3418045 11.2390137 0.7152724 10.1067667 0.8588074
## [19] 13.7026653 1.2235999 18.7577572 9.3145256 12.1820030 6.4175906
## [25] 17.3899002 6.5920024 9.3969412 15.1148882 9.2905664 10.1960773
## [31] 11.4818830 4.7357507 18.6768017 8.8097067 5.8058267 22.0470467
## [37] 14.8280821 16.5663929 0.6840876 19.4581661 6.3707232 17.1571484
## [43] 7.3425555 9.7436514 6.4165058 18.9415569 10.1223640hist(munic$slope,
main = "TOLIMA's reclassified slope",
xlab = "Slope (as a category)", breaks = 12)
aspect.geo = project(aspect, "EPSG:4326")
aspect.geo## class : SpatRaster
## dimensions : 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 : aspect
## min value : 0.2454857
## max value : 359.9559631palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent", reverse = T)
slope_reduced <- aggregate(slope.geo, fact = 2, fun = mean)
leaflet(munic) %>%
addTiles() %>%
setView(-75.23421, 5.615986, 8) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$MpNombre, "<br>",
"Slope class: ", munic$slope, "<br>")) %>%
addRasterImage(slope_reduced, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope_reduced),
title = "Terrain slope mean in TOLIMA (degrees)")Es similar a la salida anterior indicando que la mayor parte del departamento tiene unas pendientes bajas, la mayor parte es terreno plano o ligeramente inclinado, aunque se observa también como hacia Bogotá (cordillera oriental) pendientes altas.
10. Visualizing terrain aspect
munic$aspect <- exactextractr::exact_extract(aspect, 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%munic$aspect## [1] 163.4453 198.4066 137.9192 149.9774 143.1216 160.0189 185.3571 156.2115
## [9] 197.2844 184.3899 160.1416 158.8511 168.2382 210.9852 190.8018 125.2017
## [17] 124.9517 147.8497 166.7902 138.9479 158.2712 156.1744 160.5640 131.8467
## [25] 154.7618 137.9451 230.0149 137.2533 177.9112 153.5352 159.8838 159.4818
## [33] 177.9445 214.1127 198.3944 162.8340 152.2108 157.9723 160.3389 178.9442
## [41] 167.6194 152.4690 260.1845 161.6029 128.1215 150.0227 248.7309hist(munic$aspect,
main = "TOLIMA's reclassified aspect",
xlab = "aspect degrees", breaks = 12)
aspect.geo = project(aspect, "EPSG:4326")
aspect.geo## class : SpatRaster
## dimensions : 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 : aspect
## min value : 0.2454857
## max value : 359.9559631aspect_reduced <- aggregate(aspect.geo, fact = 2, fun = mean)
pal_aspect <- colorNumeric(
palette = c("blue", "lightblue", "green", "yellow", "orange", "red"),
domain = c(0,360),
na.color = "transparent"
)
leaflet(munic) %>%
addTiles() %>%
setView(-75.23421, 5.615986, 8) %>%
addPolygons(
color = "gray",
weight = 1.0,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste("Municipio: ", munic$MpNombre, "<br>",
"Aspect: ", munic$aspect, "<br>")
) %>%
addRasterImage(aspect_reduced, colors = pal_aspect, opacity = 0.8) %>%
addLegend(pal = pal_aspect, values = values(aspect_reduced),
title = "Terrain Aspect in TOLIMA (degrees)")