#Introducción En este cuaderno sabremos cómo calcular atributos geomorfométricos del terreno utilizando Modelos Digitales de Elevación (MDE). Estos atributos son esenciales para comprender la topografía y su influencia en diversos procesos agropecuarios.
Empezamos limpiando la memoria
rm(list=ls())
A continuacion instalamos las librerias necesarias para las lecturas de los codigos, segun las caracteristicas de este proceso de modelos digitales de elevación.
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 TRUE
library(leaflet)
## Warning: package 'leaflet' was built under R version 4.3.3
library(terra)
## Warning: package 'terra' was built under R version 4.3.3
## terra 1.8.15
library(MultiscaleDTM)
## Warning: package 'MultiscaleDTM' was built under R version 4.3.3
library(exactextractr)
## Warning: package 'exactextractr' was built under R version 4.3.3
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.
(dem = suppressWarnings(terra::rast("Datos/elev_boyaca.tif")))
## class : SpatRaster
## dimensions : 287, 325, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -74.65833, -71.95, 4.658333, 7.05 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : elev_boyaca.tif
## name : elev_boyaca
Este comando permite leer el archivo ráster y obtener información básica del DEM, como su extensión y resolución. Para reducir la resolución del DEM y evitar problemas de memoria, estás utilizando el siguiente código:
dem2 = terra::aggregate(dem,2, "mean")
Este comando agrupa las celdas del DEM original en bloques de 2x2, promediando sus valores, lo que reduce la resolución y, por lo tanto, la carga de memoria.
(munic <- sf::st_read("Datos/Depto.gpkg"))
## Reading layer `Depto' from data source
## `C:\Users\santi\OneDrive\Documentos\Parcial_final_july\Datos\Depto.gpkg'
## using driver `GPKG'
## Simple feature collection with 123 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: MAGNA-SIRGAS
Al recortar el DEM ( dem2) a los límites del departamento usando el siguiente código:
print(st_crs(munic))
## Coordinate Reference System:
## User input: MAGNA-SIRGAS
## wkt:
## GEOGCRS["MAGNA-SIRGAS",
## DATUM["Marco Geocentrico Nacional de Referencia",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
## BBOX[-4.23,-84.77,15.51,-66.87]],
## ID["EPSG",4686]]
print(crs(dem2))
## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
munic <- st_transform(munic, crs(dem2))
print(st_crs(munic))
## Coordinate Reference System:
## User input: GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
print(crs(dem2))
## [1] "GEOGCRS[\"WGS 84\",\n ENSEMBLE[\"World Geodetic System 1984 ensemble\",\n MEMBER[\"World Geodetic System 1984 (Transit)\"],\n MEMBER[\"World Geodetic System 1984 (G730)\"],\n MEMBER[\"World Geodetic System 1984 (G873)\"],\n MEMBER[\"World Geodetic System 1984 (G1150)\"],\n MEMBER[\"World Geodetic System 1984 (G1674)\"],\n MEMBER[\"World Geodetic System 1984 (G1762)\"],\n MEMBER[\"World Geodetic System 1984 (G2139)\"],\n ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1]],\n ENSEMBLEACCURACY[2.0]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n CS[ellipsoidal,2],\n AXIS[\"geodetic latitude (Lat)\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n AXIS[\"geodetic longitude (Lon)\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433]],\n USAGE[\n SCOPE[\"Horizontal component of 3D system.\"],\n AREA[\"World.\"],\n BBOX[-90,-180,90,180]],\n ID[\"EPSG\",4326]]"
dem3 <- terra::crop(dem2, munic, mask = TRUE)
(munic <- sf::st_read("Datos/Depto.gpkg"))
## Reading layer `Depto' from data source
## `C:\Users\santi\OneDrive\Documentos\Parcial_final_july\Datos\Depto.gpkg'
## using driver `GPKG'
## Simple feature collection with 123 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: MAGNA-SIRGAS
El siguiente código transforma el DEM a este sistema de coordenadas:
(dem_plane = project(dem3, "EPSG:9377"))
## class : SpatRaster
## dimensions : 144, 163, 1 (nrow, ncol, nlyr)
## resolution : 1843.218, 1843.218 (x, y)
## extent : 4816120, 5116565, 2071916, 2337339 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : elev_boyaca
## min value : 1
## max value : 20
Para transformar el objeto vectorial, se utiliza:
(munic_plane = sf::st_transform(munic, "EPSG:9377"))
Para calcular la pendiente y el aspecto, se utiliza la función SlpAsp del paquete MultiscaleDTM .
Se define el parámetro w como el tamaño de la ventana utilizada para calcular los atributos. En este caso, w = c(3, 3)corresponde a una ventana de celdas de 3x3.
El parámetro methoddefine el tipo de vecindad utilizado para los cálculos, donde “queen” es un tipo de vecindad que incluye los vecinos directos e intermedios.
El siguiente código calcula la pendiente y el aspecto:
(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 : 144, 163, 2 (nrow, ncol, nlyr)
## resolution : 1843.218, 1843.218 (x, y)
## extent : 4816120, 5116565, 2071916, 2337339 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.0000000, 0.01479092
## max values : 0.2742479, 359.97650023
Luego, se separan los dos resultados en dos capas: pendiente y aspecto.En este caso para la pendiente:
(slope = subset(slp_asp, 1))
## class : SpatRaster
## dimensions : 144, 163, 1 (nrow, ncol, nlyr)
## resolution : 1843.218, 1843.218 (x, y)
## extent : 4816120, 5116565, 2071916, 2337339 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000000
## max value : 0.2742479
Por último, se observa la distribución de los valores de la pendiente.
terra::hist(slope,
main = "Pendiente Boyaca",
xlab = "Pendiente (en grados)",
ylab = "Frecuencia",
col = "#B4EEB4",
border = "gray",
cex.main = 1.8, # tamaño del título
cex.lab = 1.2, # tamaño subtítulos
cex.axis = 1, # tamaño números
sub = "Distribucion de la pendiente en el departamento de Boyaca",
col.sub = "darkred",
font.sub = 4)
Con esto, se obtiene la distribución de las pendientes para el
departamento de Boyacá
(aspect =subset(slp_asp, 2))
## class : SpatRaster
## dimensions : 144, 163, 1 (nrow, ncol, nlyr)
## resolution : 1843.218, 1843.218 (x, y)
## extent : 4816120, 5116565, 2071916, 2337339 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.01479092
## max value : 359.97650023
terra::hist(aspect,
main = "Orientacion de Boyaca",
xlab = "Orientación (en grados)",
ylab = "Frecuencia",
col = "powderblue",
border = "white",
cex.main = 1.8, # tamaño del título
cex.lab = 1.2, # tamaño subtítulos
cex.axis = 1, # tamaño números
sub = "Distribución de la orientación en el departamento de Boyaca",
col.sub = "darkgreen",
font.sub = 4)
(slope_perc = tan(slope*(pi/180))*100)
## class : SpatRaster
## dimensions : 144, 163, 1 (nrow, ncol, nlyr)
## resolution : 1843.218, 1843.218 (x, y)
## extent : 4816120, 5116565, 2071916, 2337339 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000000
## max value : 0.4786566
terra::hist(slope_perc,
main = "Pendiente de Boyaca",
xlab = "Pendiente (en porcentaje)",
ylab = "Frecuencia",
col = "#FFB5C5",
border = "white",
cex.main = 1.8, # tamaño del título
cex.lab = 1.2, # tamaño subtítulos
cex.axis = 1, # tamaño números
sub = "Distribución de la pendiente en porcentaje en el departamento de Boyaca",
col.sub = "darkblue",
font.sub = 4)
7.Reclasificación del raster de pendiente en Boyaca
Para analizar las pendientes en los municipios de Boyacá, se emplea la clasificación del IGAC diseñada con fines agrológicos. Esta clasificación organiza los valores de pendiente en distintos rangos, los cuales reflejan el nivel de dificultad para el manejo del terreno. Se utiliza una matriz para establecer estos rangos y asignarles una clase numérica. Posteriormente, esta matriz se aplica al ráster de pendiente (slope_perc) para generar un nuevo ráster reclasificado.
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)
Esta reclasificación facilita la identificación de áreas con pendientes similares, permitiendo categorizarlas de acuerdo con su nivel de inclinación.
7.1 Para comprender cómo varían las pendientes en cada municipio de Boyacá, se calculan dos estadísticas clave:
Pendiente media por municipio: Se determina el valor promedio de la pendiente (en porcentaje) dentro de los límites de cada municipio. Esto permite comparar la inclinación general entre las distintas áreas de Boyacá.
Clase de pendiente predominante: Se identifica la clase de pendiente más frecuente en cada municipio, basándose en la clasificación del IGAC. Esto ayuda a entender qué tipo de pendiente es dominante en cada zona del departamento.
(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%
|
|======== | 11%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
## [1] 0.03824796 0.10228322 0.12188619 0.17257851 0.09248400 0.16334170
## [7] 0.04408044 0.08056268 0.01714049 0.14620163 0.07584477 0.04778641
## [13] 0.01998147 0.12053613 0.06858964 0.12052519 0.04400976 0.08060946
## [19] 0.09334166 0.17445140 0.03495606 0.07320797 0.05173586 0.17509906
## [25] 0.04729670 0.04331277 0.14154030 0.05228501 0.05650099 0.08881623
## [31] 0.18582866 0.10216649 0.03833278 0.06514294 0.02868063 0.06656063
## [37] 0.06456957 0.07168100 0.14991330 0.13040848 0.05356707 0.20959651
## [43] 0.06338517 0.03033858 0.03072923 0.01495125 0.16282517 0.11404359
## [49] 0.17378142 0.05188702 0.11122964 0.12968807 0.12505963 0.11071897
## [55] 0.12435024 0.05218230 0.18595566 0.02713498 0.16824986 0.06311592
## [61] 0.04350354 0.03693572 0.18502790 0.12473734 0.15123615 0.10707642
## [67] 0.20377566 0.05237078 0.14232056 0.14135656 0.06100240 0.09674646
## [73] 0.15474094 0.05277750 0.15133597 0.11645130 0.03003514 0.13521867
## [79] 0.04729862 0.02783820 0.05754092 0.16355187 0.15168579 0.08792827
## [85] 0.17757371 0.06089607 0.14813828 0.09615395 0.13507497 0.06216850
## [91] 0.05217327 0.04530768 0.03414519 0.07493364 0.08399057 0.10441225
## [97] 0.03655179 0.04477508 0.13539398 0.04924262 0.14722648 0.03367342
## [103] 0.07117687 0.04674789 0.06600705 0.06087670 0.10294863 0.06395403
## [109] 0.03239425 0.05041589 0.06653506 0.05131621 0.15195832 0.04249551
## [115] 0.09125347 0.24388927 0.02403909 0.02738918 0.07652663 0.07915299
## [121] 0.04765162 0.08762407 0.14786936
rc <- rast("Datos/elev_boyaca.tif")
if (exists("rc")) {
print(st_crs(rc))
} else {
print("El objeto 'rc' no está definido.")}
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
print(st_crs(munic))
## Coordinate Reference System:
## User input: MAGNA-SIRGAS
## wkt:
## GEOGCRS["MAGNA-SIRGAS",
## DATUM["Marco Geocentrico Nacional de Referencia",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
## BBOX[-4.23,-84.77,15.51,-66.87]],
## ID["EPSG",4686]]
print(st_crs(munic))
## Coordinate Reference System:
## User input: MAGNA-SIRGAS
## wkt:
## GEOGCRS["MAGNA-SIRGAS",
## DATUM["Marco Geocentrico Nacional de Referencia",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
## BBOX[-4.23,-84.77,15.51,-66.87]],
## ID["EPSG",4686]]
print(st_crs(rc))
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
munic$class <- exactextractr::exact_extract(rc, munic, 'mode')
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:4326)
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
ocultar la barra
progress = FALSE
7.2 Visualización: Histograma de la pendiente media en Boyacá
Finalmente, para visualizar cómo se distribuyen los valores de pendientes promedio entre los municipios de Boyacá, se crea un histograma. Este gráfico permite identificar patrones generales, como si la mayoría de los municipios presentan pendientes bajas, moderadas o pronunciadas.
hist(munic$mean_slope,
main = "Pendiente Media en Boyacá",
xlab = "Pendiente (en porcentaje)",
ylab = "Frecuencia de municipios",
col = "#FFEC8B",
border = "white",
cex.main = 1.8, # Tamaño del título
cex.lab = 1.2, # Tamaño de las etiquetas de los ejes
cex.axis = 1, # Tamaño de los números en los ejes
sub = "Distribución de la pendiente media en porcentaje en el departamento de Boyacá",
col.sub = "darkgreen",
font.sub = 4)
El histograma permite visualizar de manera clara cuántos municipios presentan pendientes promedio dentro de rangos específicos, lo que simplifica la interpretación visual de la información geográfica.
7.3 Análisis de la pendiente en los municipios de Nariño Estadísticas zonales para los municipios de Boyacá:
A partir del dataset de los municipios de Boyacá, se determina la clase de pendiente predominante en cada municipio utilizando un ráster de pendiente reclasificado. Este cálculo se realiza mediante el método de la moda (mode), que identifica la categoría de pendiente más frecuente dentro de los límites de cada municipio.
munic$class <- exactextractr::exact_extract(rc, munic, 'mode')
## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:4326)
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|=== | 4%
|
|=== | 5%
|
|==== | 6%
|
|===== | 7%
|
|====== | 8%
|
|====== | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 17%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 23%
|
|================= | 24%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 30%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 36%
|
|========================== | 37%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================== | 42%
|
|============================== | 43%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 50%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 63%
|
|============================================= | 64%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 70%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 77%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 83%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
Esta información es fundamental para planificar actividades agrícolas en las zonas montañosas de Boyacá y para identificar áreas donde la pendiente podría representar un riesgo de erosión o dificultar el aprovechamiento del suelo.
Histograma de clases de pendiente predominante: Se crea un histograma para visualizar la distribución de las clases de pendiente predominante en los municipios de Boyacá. Esto permite identificar patrones generales en el relieve del departamento.
hist(munic$class,
main = "Pendiente reclasificada en Boyacá",
xlab = "Clase de pendiente",
ylab = "Frecuencia de municipios",
col = "#DDA0DD",
border = "white",
cex.main = 1.8, # Tamaño del título
cex.lab = 1.2, # Tamaño de las etiquetas de los ejes
cex.axis = 1, # Tamaño de los números en los ejes
sub = "Distribución de la clase de pendiente reclasificada en el departamento de Boyacá",
col.sub = "darkblue",
font.sub = 4)
(rc.geo = project(rc, "EPSG:4326"))
## class : SpatRaster
## dimensions : 287, 325, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -74.65833, -71.95, 4.658333, 7.05 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : elev_boyaca
## min value : 1
## max value : 20
La pendiente en Boyacá muestra una gran variabilidad debido a su ubicación en la cordillera de los Andes, por lo que emplear un sistema de referencia estándar permite integrar de manera más eficiente esta información con datos de otras áreas regionales.
(slope.geo = project(slope_perc, "EPSG:4326"))
## class : SpatRaster
## dimensions : 144, 163, 1 (nrow, ncol, nlyr)
## resolution : 0.01666666, 0.01666666 (x, y)
## extent : -74.66542, -71.94876, 4.652944, 7.052943 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000000
## max value : 0.4786566
En este análisis, se crea un mapa interactivo utilizando la librería leaflet en R para visualizar la pendiente del terreno en los municipios de Boyacá. Este tipo de visualización es útil para comprender las características topográficas del departamento y apoyar la toma de decisiones relacionadas con el uso del suelo, la gestión de riesgos naturales y la planificación territorial.
8.1. Definición de la paleta de colores Se define una paleta de colores que representa los distintos niveles de pendiente, desde las más bajas (verde) hasta las más altas (rojo). Esto permite identificar visualmente las áreas con mayor y menor inclinación del terreno en Boyacá.
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"),
values(slope.geo),
na.color = "transparent")
8.2. Configuración y creación del mapa interactivo El mapa se centra en los municipios de Boyacá, permitiendo visualizar:
Límites municipales: Cada municipio se resalta con un borde gris, y la opacidad de relleno se ajusta para destacar las capas subyacentes.
Pendiente del terreno: Representada por el ráster slope.geo, utilizando la paleta de colores definida previamente.
Información interactiva: Al hacer clic en un municipio, se desplegará un cuadro emergente con el nombre del municipio y la clase de pendiente predominante.
st_crs(munic)
## Coordinate Reference System:
## User input: MAGNA-SIRGAS
## wkt:
## GEOGCRS["MAGNA-SIRGAS",
## DATUM["Marco Geocentrico Nacional de Referencia",
## ELLIPSOID["GRS 1980",6378137,298.257222101,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["Colombia - onshore and offshore. Includes San Andres y Providencia, Malpelo Islands, Roncador Bank, Serrana Bank and Serranilla Bank."],
## BBOX[-4.23,-84.77,15.51,-66.87]],
## ID["EPSG",4686]]
munic <- st_transform(munic, 4326)
st_crs(munic) <- 4326
El código para generar este mapa sería:
class(munic)
## [1] "sf" "data.frame"
st_geometry(munic)
## Geometry set for 123 features
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.66496 ymin: 4.655196 xmax: -71.94885 ymax: 7.055557
## Geodetic CRS: WGS 84
## First 5 geometries:
## MULTIPOLYGON (((-73.34014 5.583085, -73.34011 5...
## MULTIPOLYGON (((-73.36793 5.01349, -73.36672 5....
## MULTIPOLYGON (((-72.76309 5.643222, -72.76309 5...
## MULTIPOLYGON (((-73.50487 5.843478, -73.50451 5...
## MULTIPOLYGON (((-72.91692 6.086122, -72.91712 6...
st_crs(munic)
## Coordinate Reference System:
## User input: EPSG:4326
## wkt:
## GEOGCRS["WGS 84",
## ENSEMBLE["World Geodetic System 1984 ensemble",
## MEMBER["World Geodetic System 1984 (Transit)"],
## MEMBER["World Geodetic System 1984 (G730)"],
## MEMBER["World Geodetic System 1984 (G873)"],
## MEMBER["World Geodetic System 1984 (G1150)"],
## MEMBER["World Geodetic System 1984 (G1674)"],
## MEMBER["World Geodetic System 1984 (G1762)"],
## MEMBER["World Geodetic System 1984 (G2139)"],
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ENSEMBLEACCURACY[2.0]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["geodetic latitude (Lat)",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["geodetic longitude (Lon)",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## USAGE[
## SCOPE["Horizontal component of 3D system."],
## AREA["World."],
## BBOX[-90,-180,90,180]],
## ID["EPSG",4326]]
st_crs(munic) <- 4326
sum(st_is_empty(munic))
## [1] 0
leaflet(munic) %>% addTiles() %>% setView(-72.6, 5.8, 8) %>%
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 Cordoba (%)")