Examen 2
Oscar Duvan Diaz Amaya
2025-01-23
1. Introduccion
Este cuaderno ilustra cĆ³mo calcular atributos geomorfomĆ©tricos del terreno a partir de modelos de elevaciĆ³n digitales cuadriculados para el departamento de Cundinamarca.
2. Setup
Instalamos y llamamos librerias para el desarrollo de nuestro trabajo.
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)# Limpiamos la memoria
rm(list=ls())3. Introduccion a MultiscaleDTM
Con esta funcion podemos calcular distintos atributos geomorfomĆ©tricos del terreno del departamento de Cundinamarca por municipios, con mĆŗltiples escalas a partir de modelos digitales cuadriculados regularmente (DTM; es decir, rĆ”steres de elevaciĆ³n).
4. Leer datos de entrada
Tomamos los datos de altitud (metros) de los municipios de Cundinamarca.
(dem= terra::rast("cundinamarca/cundinamarca Elevacion_Cundinamarca.tif"))## class : SpatRaster
## dimensions : 3578, 3590, 1 (nrow, ncol, nlyr)
## resolution : 0.0006854461, 0.0006854461 (x, y)
## extent : -75.23438, -72.77362, 3.513228, 5.965754 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +no_defs
## source : cundinamarca Elevacion_Cundinamarca.tif
## name : cundinamarca Elevacion_Cundinamarca
## min value : -643
## max value : 4267dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Este codigo nos sirve para reducir la resoluciĆ³n del raster (hace que las celdas sean mĆ”s grandes y menos detalladas). Ademas suaviza los valores al calcular la media en bloques de 2x2 celdas.Y nos sirve para mejorar el rendimiento en cĆ”lculos espaciales, ya que el nuevo raster tiene menos celdas .
(munic <- sf::st_read("cundinamarca/cundinamarca mun_Cun.gpkg"))## Reading layer `mgn_adm_mpio_grafico' from data source
## `C:\Users\diazo\OneDrive\Documentos\GEOMATICA\PROYECTO_4\cundinamarca\cundinamarca mun_Cun.gpkg'
## using driver `GPKG'
## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: MAGNA-SIRGAS## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -74.89063 ymin: 3.730129 xmax: -73.05256 ymax: 5.837258
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 25 001 25001 CUNDINAMARCA AGUA DE DIOS
## 2 25 019 25019 CUNDINAMARCA ALBĆN
## 3 25 035 25035 CUNDINAMARCA ANAPOIMA
## 4 25 040 25040 CUNDINAMARCA ANOLAIMA
## 5 25 053 25053 CUNDINAMARCA ARBELĆEZ
## 6 25 086 25086 CUNDINAMARCA BELTRĆN
## 7 25 095 25095 CUNDINAMARCA BITUIMA
## 8 25 099 25099 CUNDINAMARCA BOJACĆ
## 9 25 120 25120 CUNDINAMARCA CABRERA
## 10 25 123 25123 CUNDINAMARCA CACHIPAY
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 78 Noviembre 29 de 1963 MUNICIPIO 85.95871 2023
## 2 Ordenanza 19 del 22 de Mayo de 1903 MUNICIPIO 50.83510 2023
## 3 DispocisiĆ³n de Septiembre 1 de 1864 MUNICIPIO 123.95549 2023
## 4 1882 MUNICIPIO 120.95146 2023
## 5 Decreto 32 de Enero 16 de 1886 MUNICIPIO 142.45716 2023
## 6 Ordenanza 5 del 12 de Noviembre de 1853 MUNICIPIO 176.85320 2023
## 7 1772 MUNICIPIO 61.04708 2023
## 8 1700 MUNICIPIO 102.31988 2023
## 9 Ordenanza 79 Noviembre 29 de 1963 MUNICIPIO 421.83959 2023
## 10 Ordenanza 6 Noviembre 26 de 1982 MUNICIPIO 53.45262 2023
## shape_Leng shape_Area geom
## 1 0.4963557 0.007002449 MULTIPOLYGON (((-74.67275 4...
## 2 0.5162142 0.004144411 MULTIPOLYGON (((-74.47852 4...
## 3 0.6705565 0.010100758 MULTIPOLYGON (((-74.55201 4...
## 4 0.5280079 0.009859268 MULTIPOLYGON (((-74.43043 4...
## 5 0.8065843 0.011603784 MULTIPOLYGON (((-74.42848 4...
## 6 0.9093670 0.014413273 MULTIPOLYGON (((-74.76139 4...
## 7 0.4907957 0.004976541 MULTIPOLYGON (((-74.51062 4...
## 8 0.4979169 0.008339597 MULTIPOLYGON (((-74.32819 4...
## 9 1.0188051 0.034347363 MULTIPOLYGON (((-74.46463 4...
## 10 0.4211104 0.004356744 MULTIPOLYGON (((-74.40513 4...En este caso tomamos los datos de los municipios de cundinamarca en un archivo shp o gpgk como en este caso para que nos sirva.
dem3 = terra::crop(dem2,munic, mask=TRUE)Aca lo que se hizo fue extraer solo una regiĆ³n de un rĆ”ster grande. Ademas hacemos mĆ”s eficiente el anĆ”lisis, eliminando datos innecesarios.
5. Transforming coordinates
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 1537, 1345, 1 (nrow, ncol, nlyr)
## resolution : 151.708, 151.708 (x, y)
## extent : 4790148, 4994195, 1970141, 2203316 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : cundinamarca Elevacion_Cundinamarca
## min value : 37.13568
## max value : 3987.95483(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 116 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4790252 ymin: 1970314 xmax: 4994174 ymax: 2203181
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## dpto_ccdgo mpio_ccdgo mpio_cdpmp dpto_cnmbr mpio_cnmbr
## 1 25 001 25001 CUNDINAMARCA AGUA DE DIOS
## 2 25 019 25019 CUNDINAMARCA ALBĆN
## 3 25 035 25035 CUNDINAMARCA ANAPOIMA
## 4 25 040 25040 CUNDINAMARCA ANOLAIMA
## 5 25 053 25053 CUNDINAMARCA ARBELĆEZ
## 6 25 086 25086 CUNDINAMARCA BELTRĆN
## 7 25 095 25095 CUNDINAMARCA BITUIMA
## 8 25 099 25099 CUNDINAMARCA BOJACĆ
## 9 25 120 25120 CUNDINAMARCA CABRERA
## 10 25 123 25123 CUNDINAMARCA CACHIPAY
## mpio_crslc mpio_tipo mpio_narea mpio_nano
## 1 Ordenanza 78 Noviembre 29 de 1963 MUNICIPIO 85.95871 2023
## 2 Ordenanza 19 del 22 de Mayo de 1903 MUNICIPIO 50.83510 2023
## 3 DispocisiĆ³n de Septiembre 1 de 1864 MUNICIPIO 123.95549 2023
## 4 1882 MUNICIPIO 120.95146 2023
## 5 Decreto 32 de Enero 16 de 1886 MUNICIPIO 142.45716 2023
## 6 Ordenanza 5 del 12 de Noviembre de 1853 MUNICIPIO 176.85320 2023
## 7 1772 MUNICIPIO 61.04708 2023
## 8 1700 MUNICIPIO 102.31988 2023
## 9 Ordenanza 79 Noviembre 29 de 1963 MUNICIPIO 421.83959 2023
## 10 Ordenanza 6 Noviembre 26 de 1982 MUNICIPIO 53.45262 2023
## shape_Leng shape_Area geom
## 1 0.4963557 0.007002449 MULTIPOLYGON (((4814462 204...
## 2 0.5162142 0.004144411 MULTIPOLYGON (((4836135 210...
## 3 0.6705565 0.010100758 MULTIPOLYGON (((4827912 207...
## 4 0.5280079 0.009859268 MULTIPOLYGON (((4841442 209...
## 5 0.8065843 0.011603784 MULTIPOLYGON (((4841536 203...
## 6 0.9093670 0.014413273 MULTIPOLYGON (((4804743 209...
## 7 0.4907957 0.004976541 MULTIPOLYGON (((4832575 210...
## 8 0.4979169 0.008339597 MULTIPOLYGON (((4852759 208...
## 9 1.0188051 0.034347363 MULTIPOLYGON (((4837472 200...
## 10 0.4211104 0.004356744 MULTIPOLYGON (((4844230 208...En este caso lo que se hizo fue transformar las coordenadas a ESPG:ā9377ā.
6. Computar atributos del terreno
slp_asp = SlpAsp(dem_plane,
w = c(3, 3),
unit = "degrees",
method = "queen",
metrics = c("slope", "aspect", "eastness", "northness"),
na.rm = FALSE,
include_scale = FALSE,
mask_aspect = TRUE,
filename = NULL,
overwrite = FALSE,
wopt = list()
)Con este codigo podemos calcular diversas mĆ©tricas de relieve como pendiente y aspecto basadas en un Modelo Digital de ElevaciĆ³n (DEM)
Capa 1
(slope = terra::subset(slp_asp, 1))## class : SpatRaster
## dimensions : 1537, 1345, 1 (nrow, ncol, nlyr)
## resolution : 151.708, 151.708 (x, y)
## extent : 4790148, 4994195, 1970141, 2203316 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 75.47234terra::hist(slope,
main = "Cundinamarca`s slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 48% of the cells was used (of which 53% was NA)
La grafica anterior permite visualizar cĆ³mo se distribuyen las pendientes en la regiĆ³n de Cundinamarca, el terreno es mas que todo plano debido a que la matoria de los datos son bajos.
Capa 2
(aspect= subset(slp_asp, 2))## class : SpatRaster
## dimensions : 1537, 1345, 1 (nrow, ncol, nlyr)
## resolution : 151.708, 151.708 (x, y)
## extent : 4790148, 4994195, 1970141, 2203316 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 0.0000
## max value : 359.9989terra::hist(aspect,
main = "Cundinamarca's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 48% of the cells was used (of which 53% was NA)
La grafica anterior evidencia la orientaciĆ³n (aspecto) del municipio de Cundinamarca, y nos evidencia que esta centrado en el pais debido a la variabilidad de datos. # Convertimos los grados de pendiente en porcentaje
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## dimensions : 1537, 1345, 1 (nrow, ncol, nlyr)
## resolution : 151.708, 151.708 (x, y)
## extent : 4790148, 4994195, 1970141, 2203316 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 385.9027terra::hist(slope_perc,
main = "Cundinamarca's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 48% of the cells was used (of which 53% was NA)
La grafica anterior nos muestra los valores de pendiente en porcentaje que se calculan con la fĆ³rmula: Pendiente en porcentaje = tan (pendiente en grados)Ć100 Esto significa que la pendiente se expresa como una relaciĆ³n entre el cambio de altura y la distancia horizontal, multiplicada por 100, como se vio en clase.
7. Calcular estadisticas zonales
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)Aca clasificamos por clases a los valores de los intervalos de pendiente correspondientes.
(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% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 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% | |======================================================= | 78% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 9.2789555 24.5619545 17.8890781 23.0709496 20.3946609 16.0540905
## [7] 31.9070873 16.0421371 20.4668007 20.3224583 8.3369131 30.4048481
## [13] 31.3308182 19.8978462 26.2327480 9.2998056 23.0421772 26.2894878
## [19] 14.5278978 15.0135031 10.2309313 16.4060421 21.2737465 29.1958122
## [25] 8.4426174 10.6956606 24.3487587 36.5306015 0.7304142 8.9482422
## [31] 19.5193539 35.0954018 11.7416220 25.9359512 32.0816460 12.4349642
## [37] 16.4772072 12.5185528 21.0935307 16.0828953 21.9922333 16.5229702
## [43] 25.2329445 48.5213966 43.4887161 22.0672913 24.3876438 20.9459171
## [49] 18.8080273 31.9185276 29.7840252 26.6676540 16.0642433 24.5273533
## [55] 4.1809587 28.7582703 25.6633453 4.6606512 21.3899345 9.3765745
## [61] 17.5166206 32.5096741 26.0836811 21.3995514 26.1463566 32.6371422
## [67] 25.5408764 8.4538078 22.0645447 7.9754252 24.0922394 29.9205189
## [73] 40.7420959 27.1739693 23.0858212 9.8030128 28.0060501 21.5570679
## [79] 27.3691025 25.3083839 24.2922192 24.6100788 14.7060375 15.6528111
## [85] 20.3410358 11.2807369 14.0337925 12.1420918 15.5307016 13.3509054
## [91] 25.9112644 15.4044952 19.4499454 15.0191574 16.6925240 24.1082783
## [97] 6.8990779 26.5358372 25.3646297 14.8658962 8.7521000 25.5627384
## [103] 29.4129009 24.8152695 10.9253435 23.1209049 27.4529915 28.5834007
## [109] 31.1060753 28.5162868 16.1921844 24.4757900 21.9290314 29.0733757
## [115] 22.7351818 15.3128786hist(munic$mean_slope,
main = "Cundinamarca's mean slope",
xlab = "Slope (in porcentaje)")
La grafica anterior nos muestra la clasificacion para la pendiente (en porcentaje) segun las clases en las cuales se dividieron los intervales de datos.
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |===== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |============== | 21% | |=============== | 22% | |================ | 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% | |======================================================= | 78% | |======================================================== | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================= | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%## [1] 2 4 4 4 4 4 5 4 4 4 1 5 5 4 5 1 4 4 4 4 1 4 4 5 1 1 4 5 1 1 4 5 1 4 5 1 4
## [38] 4 4 4 5 4 4 5 5 4 4 4 4 5 5 5 4 4 1 5 4 1 5 1 4 5 5 4 4 5 4 1 4 1 4 5 5 5
## [75] 4 1 4 4 4 4 4 4 4 4 4 1 4 1 4 4 4 4 4 4 4 4 1 4 4 4 1 5 5 4 4 4 5 5 5 5 4
## [112] 4 4 5 4 4hist(munic$class,
main = "Cundinamarca's reclassified slope",
xlab = "Slope (as a category)")
La grafica anetior representa cuantos municipios estan dentro de una clase de pendiente determinada segun las clases anteriores.
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1541, 1345, 1 (nrow, ncol, nlyr)
## resolution : 0.001370884, 0.001370884 (x, y)
## extent : -74.89607, -73.05223, 3.727532, 5.840064 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 363.2341(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1541, 1345, 1 (nrow, ncol, nlyr)
## resolution : 0.001370884, 0.001370884 (x, y)
## extent : -74.89607, -73.05223, 3.727532, 5.840064 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 363.23418. Trazado de pendiente del terreno
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")leaflet(munic) %>% addTiles() %>% setView(-75.9, 8.7, 9) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Slope class: ", munic$class, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Terrain slope in Cundinamarca (%)")## 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'