Atributos Geomorfometricos
Daniel Diaz Vargas
2025-11-04
1 INTRODUCCION
Este cuaderno tiene como propocito ilustrar el como calcular atributos geomorfometricos en el terreno a partir de modelos digitales de elevacion en formato de malla. como objetivo se establece entender los conceptos y herramientas basicas de geoinformatica, con la utilizacion de paquetes MultiscaleDTM de R.
2 Configuracion
rm(list=ls())LLAMAR LIBRERIAS:
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.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUElibrary(leaflet)
library(terra)## terra 1.8.70library(MultiscaleDTM)
library(exactextractr)Leemos el DEM con TERRA
`
list.files("datos01")## [1] "DEM _Practica.gpkg" "Depa Huila.tif.xml" "ELE HUILA.tif"
## [4] "ELE HUILA.tif.xml" "ele_huila.tif.aux.xml" "munic_huila.gpkg"
## [7] "tabla.png"4 Lectura de datos de entrada
(dem = terra::rast("datos01/ELE HUILA.tif"))## class : SpatRaster
## size : 2416, 1795, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -81.83333, -66.875, -4.225, 15.90833 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source : ELE HUILA.tif
## name : ELE HUILAReduzcamos la resolución del DEM para evitar problemas de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= dem2## class : SpatRaster
## size : 1208, 898, 1 (nrow, ncol, nlyr)
## resolution : 0.01666667, 0.01666667 (x, y)
## extent : -81.83333, -66.86667, -4.225, 15.90833 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## name : ELE HUILA
## min value : -2.75
## max value : 5305.75Ahora, leeremos sobre nuestros municipios:
(munic <- vect("datos01/munic_huila.gpkg"))## class : SpatVector
## geometry : polygons
## dimensions : 37, 11 (geometries, attributes)
## extent : -76.6402, -74.5054, 1.493471, 3.807901 (xmin, xmax, ymin, ymax)
## source : munic_huila.gpkg (col_adm2)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## names : ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2
## type : <num> <chr> <chr> <num> <chr> <num> <chr> <chr>
## values : 53 COL Colombia 17 Huila 616 Acevedo Municipio
## 53 COL Colombia 17 Huila 617 Agrado Municipio
## 53 COL Colombia 17 Huila 618 Aipe Municipio
## ENGTYPE_2 NL_NAME_2 VARNAME_2
## <chr> <chr> <chr>
## Municipality NA NA
## Municipality NA NA
## Municipality NA NAAhora, recortemos el objeto dem2 a los límites de nuestro departamento:
munic_huila <- st_read("datos01/munic_huila.gpkg")## Reading layer `col_adm2' from data source
## `D:\GB2\Proyecto 4\datos01\munic_huila.gpkg' using driver `GPKG'
## Simple feature collection with 37 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.6402 ymin: 1.493471 xmax: -74.5054 ymax: 3.807901
## Geodetic CRS: WGS 84huila_crs <- st_read("datos01/munic_huila.gpkg")## Reading layer `col_adm2' from data source
## `D:\GB2\Proyecto 4\datos01\munic_huila.gpkg' using driver `GPKG'
## Simple feature collection with 37 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.6402 ymin: 1.493471 xmax: -74.5054 ymax: 3.807901
## Geodetic CRS: WGS 84(huila_crs = crs((munic_huila)))## [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 MEMBER[\"World Geodetic System 1984 (G2296)\"],\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]]"atribuir crs de huila al dem2
crs(dem2)## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"crs(dem2 <- huila_crs)## [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 MEMBER[\"World Geodetic System 1984 (G2296)\"],\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]]"print("CRS actual del DEM:")## [1] "CRS actual del DEM:"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 MEMBER[\"World Geodetic System 1984 (G2296)\"],\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]]"print("CRS de referencia:")## [1] "CRS de referencia:"huila_crs## [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 MEMBER[\"World Geodetic System 1984 (G2296)\"],\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]]"dem2 = terra::rast("datos01/ELE HUILA.tif")crs(dem2)## [1] "GEOGCRS[\"WGS 84\",\n DATUM[\"unknown\",\n ELLIPSOID[\"WGS84\",6378137,298.257223563,\n LENGTHUNIT[\"metre\",1,\n ID[\"EPSG\",9001]]]],\n PRIMEM[\"Greenwich\",0,\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n CS[ellipsoidal,2],\n AXIS[\"latitude\",north,\n ORDER[1],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]],\n AXIS[\"longitude\",east,\n ORDER[2],\n ANGLEUNIT[\"degree\",0.0174532925199433,\n ID[\"EPSG\",9122]]]]"(dem3 = terra::crop(dem2,munic_huila, mask=TRUE))## Warning: [crop] CRS do not match## class : SpatRaster
## size : 278, 256, 1 (nrow, ncol, nlyr)
## resolution : 0.008333333, 0.008333333 (x, y)
## extent : -76.64167, -74.50833, 1.491667, 3.808333 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84
## source(s) : memory
## varname : ELE HUILA
## name : ELE HUILA
## min value : 348
## max value : 44545 transformacion de cordenadas
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## size : 278, 257, 1 (nrow, ncol, nlyr)
## resolution : 924.0865, 924.0865 (x, y)
## extent : 4594798, 4832288, 1722781, 1979677 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : ELE HUILA
## min value : 348.00
## max value : 4407.92También necesitamos realizar una transformación similar para el objeto vectorial:
(munic_plane = project(munic, "EPSG:9377"))## class : SpatVector
## geometry : polygons
## dimensions : 37, 11 (geometries, attributes)
## extent : 4595093, 4832889, 1723306, 1978934 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## names : ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2
## type : <num> <chr> <chr> <num> <chr> <num> <chr> <chr>
## values : 53 COL Colombia 17 Huila 616 Acevedo Municipio
## 53 COL Colombia 17 Huila 617 Agrado Municipio
## 53 COL Colombia 17 Huila 618 Aipe Municipio
## ENGTYPE_2 NL_NAME_2 VARNAME_2
## <chr> <chr> <chr>
## Municipality NA NA
## Municipality NA NA
## Municipality NA NA6 Cálculo de atributos del terreno
Veamos cómo calcular la pendiente y la orientación.
Abre la consola y escribe ?SlpAsp. En la ventana de
ayuda de la derecha, se mostrará la información correspondiente.
(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
## size : 278, 257, 2 (nrow, ncol, nlyr)
## resolution : 924.0865, 924.0865 (x, y)
## extent : 4594798, 4832288, 1722781, 1979677 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.01242603, 7.782563e-03
## max values : 23.84544089, 3.599989e+02Ahora, dividamos el objeto slp_asp en dos partes.
Capa 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## size : 278, 257, 1 (nrow, ncol, nlyr)
## resolution : 924.0865, 924.0865 (x, y)
## extent : 4594798, 4832288, 1722781, 1979677 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.01242603
## max value : 23.84544089¿Cuál es la distribución de los valores de la pendiente?
terra::hist(slope,
main = "Pendiente Huila",
xlab = "Pendiente (en grados)")
Capa 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## size : 278, 257, 1 (nrow, ncol, nlyr)
## resolution : 924.0865, 924.0865 (x, y)
## extent : 4594798, 4832288, 1722781, 1979677 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 7.782563e-03
## max value : 3.599989e+02terra::hist(aspect,
main = "Orientacion Huila",
xlab = "Orientacion (en grados)")
Además, convirtamos los grados de pendiente a porcentaje de pendiente:
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## size : 278, 257, 1 (nrow, ncol, nlyr)
## resolution : 924.0865, 924.0865 (x, y)
## extent : 4594798, 4832288, 1722781, 1979677 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.02168751
## max value : 44.20002550terra::hist(slope_perc,
main = "Pendiente Huila",
xlab = "Pendiente (en porcentaje)")
7 Cálculo de estadísticas zonales
Una operación de estadística zonal calcula estadísticas sobre los valores de celda de un ráster (un ráster de valores) dentro de las zonas definidas por otro conjunto de datos. En nuestro caso, nos interesa calcular el promedio de los valores de pendiente por municipio.
Primero, veamos la clasificación de pendientes del IGAC para fines agronómicos.
Reclasificar la raster de pendiente
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, calculemos las estadísticas zonales utilizando exactextractor:
munic2 <- st_as_sf(munic)(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic2, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 11.528454 8.846818 9.475933 14.747292 7.853603 10.562951 5.947561
## [8] 15.661534 11.152996 10.692496 11.499730 10.985253 7.327530 15.529167
## [15] 8.863760 9.377044 11.431463 16.009836 10.416643 10.388485 9.927352
## [22] 8.694324 11.138761 9.311407 9.204544 12.206645 8.967743 11.400238
## [29] 14.389353 10.551877 11.698293 8.156768 14.649766 10.680236 10.600506
## [36] 2.966304 6.114223hist(munic$mean_slope,
main = "Pendiente Promedio por Municipio",
xlab = "Pendiente (en Porcentaje)")
(munic$class <- exactextractr::exact_extract(rc, munic2, 'mode'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 3 2 4 4 3 4 1 4 4 4 4 4 2 4 4 3 4 4 4 3 4 2 3 4 4 4 3 4 4 3 4 1 4 4 4 1 1hist(munic$class,
main = "Pendiente Reclasificada Huila",
xlab = "Pendiente (Como Categoria)")
Transformemos la pendiente de nuevo a coordenadas geográficas:
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## size : 279, 257, 1 (nrow, ncol, nlyr)
## resolution : 0.008333243, 0.008333243 (x, y)
## extent : -76.64845, -74.50681, 1.489761, 3.814736 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 1
## max value : 5(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## size : 279, 257, 1 (nrow, ncol, nlyr)
## resolution : 0.008333243, 0.008333243 (x, y)
## extent : -76.64845, -74.50681, 1.489761, 3.814736 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.06320455
## max value : 44.200023658 Trazando la pendiente del terreno
En primer lugar, necesitaremos una paleta de colores para la pendiente:
palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo),
na.color = "transparent")munic2## Simple feature collection with 37 features and 11 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -76.6402 ymin: 1.493471 xmax: -74.5054 ymax: 3.807901
## Geodetic CRS: WGS 84
## First 10 features:
## ID_0 ISO NAME_0 ID_1 NAME_1 ID_2 NAME_2 TYPE_2 ENGTYPE_2
## 1 53 COL Colombia 17 Huila 616 Acevedo Municipio Municipality
## 2 53 COL Colombia 17 Huila 617 Agrado Municipio Municipality
## 3 53 COL Colombia 17 Huila 618 Aipe Municipio Municipality
## 4 53 COL Colombia 17 Huila 619 Algeciras Municipio Municipality
## 5 53 COL Colombia 17 Huila 620 Altamira Municipio Municipality
## 6 53 COL Colombia 17 Huila 621 Baraya Municipio Municipality
## 7 53 COL Colombia 17 Huila 622 Campoalegre Municipio Municipality
## 8 53 COL Colombia 17 Huila 623 Colombia Municipio Municipality
## 9 53 COL Colombia 17 Huila 624 Elías Municipio Municipality
## 10 53 COL Colombia 17 Huila 625 Garzón Municipio Municipality
## NL_NAME_2 VARNAME_2 geometry
## 1 <NA> <NA> POLYGON ((-75.83184 1.74618...
## 2 <NA> <NA> POLYGON ((-75.6673 2.283, -...
## 3 <NA> <NA> POLYGON ((-75.21477 3.36582...
## 4 <NA> <NA> POLYGON ((-75.30428 2.36857...
## 5 <NA> <NA> POLYGON ((-75.9259 1.9912, ...
## 6 <NA> <NA> POLYGON ((-74.9553 3.2666, ...
## 7 <NA> <NA> POLYGON ((-75.38387 2.53194...
## 8 <NA> <NA> POLYGON ((-74.9553 3.2666, ...
## 9 <NA> <NA> POLYGON ((-76.07361 1.95269...
## 10 <NA> <NA> POLYGON ((-75.76563 2.11705...Ahora, es el momento de graficar. Asegúrate de cambiar varias opciones en el código si es necesario:
leaflet(munic2) %>% addTiles() %>% setView(-76.6402, 1.493471,6.807901) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic2$NAME_2, "<br>",
"Slope class: ", munic2$NAME_2, "<br>")) %>%
addRasterImage(slope.geo, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo),
title = "Pendiente del terreno en Huila (%)")9 Visualización alternativa de la pendiente del terreno
Visualice la pendiente del terreno en grados y añada etiquetas con los nombres de los municipios junto a sus correspondientes valores de pendiente media. Mantenga los valores emergentes tal como se muestran en el gráfico anterior.
(slope_deg.geo = project(slope, "EPSG:4326"))## class : SpatRaster
## size : 279, 257, 1 (nrow, ncol, nlyr)
## resolution : 0.008333243, 0.008333243 (x, y)
## extent : -76.64845, -74.50681, 1.489761, 3.814736 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.03621353
## max value : 23.84544182(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic2, 'mean'))## Warning in .local(x, y, ...): Polygons transformed to raster CRS (EPSG:9377)## | | | 0% | |== | 3% | |==== | 5% | |====== | 8% | |======== | 11% | |========= | 14% | |=========== | 16% | |============= | 19% | |=============== | 22% | |================= | 24% | |=================== | 27% | |===================== | 30% | |======================= | 32% | |========================= | 35% | |========================== | 38% | |============================ | 41% | |============================== | 43% | |================================ | 46% | |================================== | 49% | |==================================== | 51% | |====================================== | 54% | |======================================== | 57% | |========================================== | 59% | |============================================ | 62% | |============================================= | 65% | |=============================================== | 68% | |================================================= | 70% | |=================================================== | 73% | |===================================================== | 76% | |======================================================= | 78% | |========================================================= | 81% | |=========================================================== | 84% | |============================================================= | 86% | |============================================================== | 89% | |================================================================ | 92% | |================================================================== | 95% | |==================================================================== | 97% | |======================================================================| 100%## [1] 11.528454 8.846818 9.475933 14.747292 7.853603 10.562951 5.947561
## [8] 15.661534 11.152996 10.692496 11.499730 10.985253 7.327530 15.529167
## [15] 8.863760 9.377044 11.431463 16.009836 10.416643 10.388485 9.927352
## [22] 8.694324 11.138761 9.311407 9.204544 12.206645 8.967743 11.400238
## [29] 14.389353 10.551877 11.698293 8.156768 14.649766 10.680236 10.600506
## [36] 2.966304 6.114223munic2 <- st_as_sf(munic)munic2$label_text <- sprintf(
"%s (%.2f%%)",
munic2$NAME_2,
munic2$mean_slope
)palredgreen_deg <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope_deg.geo),
na.color = "transparent")leaflet(munic2) %>%
addTiles() %>%
setView(-75.9, 2.3, 9) %>%
addRasterImage(slope_deg.geo, colors = palredgreen_deg, opacity = 0.8,
group = "Pendiente en Grados") %>%
addLegend(pal = palredgreen_deg, values = values(slope_deg.geo),
title = "Pendiente del terreno en Huila (°)") %>%
addPolygons(
color = "gray",
weight = 1.0,
smoothFactor = 0.5,
opacity = 0.8,
fillOpacity = 0.10,
popup = paste("Municipio: ", munic2$NAME_2, "<br>",
"Pendiente Media: ", round(munic2$mean_slope, 2), "%", "<br>",
"Clasificación Modal: ", munic2$class, "<br>"),label = ~label_text,
labelOptions = labelOptions(
noHide = TRUE, # Mostrar siempre la etiqueta
direction = 'center',# Centrar la etiqueta en el polígono
textOnly = TRUE, # Solo mostrar el texto (sin fondo o borde de globo)
style = list("font-weight" = "bold",
"padding" = "2px",
"font-size" = "10px")
)
)Describe el resultado y su significado.
El mapa interactivo generado combina el análisis ráster de la pendiente del terreno en grados con las estadísticas zonales por municipio, ofreciendo una visión integral de la geomorfología de Huila. La visualización de la pendiente, a través de una escala de color (desde verdes claros en zonas planas hasta rojos intensos en áreas escarpadas), permite identificar visualmente el gradiente topográfico. Las áreas rojas y oscuras, que dominan las zonas cordilleranas, señalan terrenos con alta inclinación y, por ende, con mayores limitaciones para el uso de la tierra, como riesgo de erosión y desafíos para la agricultura mecanizada o la construcción. Por otro lado, las extensiones de colores claros corresponden a los valles y planicies aluviales del río Magdalena, aptas para el desarrollo y la producción intensiva.
Los polígonos municipales complementan esta imagen con indicadores clave obtenidos mediante la estadística zonal. Las etiquetas permanentes muestran el nombre del municipio junto a su pendiente media en porcentaje, proporcionando un valor promedio de la inclinación a nivel administrativo. Más importante aún, los pop-ups revelan la Clasificación Modal de la Pendiente (según estándares como los del IGAC), indicando la categoría de aptitud de uso de la tierra más frecuente dentro de cada municipio. Esta integración de datos es fundamental para la toma de decisiones, ya que permite a los planificadores asociar de manera directa las características topográficas dominantes con los límites administrativos, facilitando la zonificación, la gestión de riesgos y la planificación del desarrollo rural y de infraestructura en el departamento.
10 Visualización del aspecto del terreno
El ráster ‘aspect’ está en grados y ya fue calculado en la sección 6
(aspect.geo = project(aspect, "EPSG:4326"))## class : SpatRaster
## size : 279, 257, 1 (nrow, ncol, nlyr)
## resolution : 0.008333243, 0.008333243 (x, y)
## extent : -76.64845, -74.50681, 1.489761, 3.814736 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : aspect
## min value : 1.26984
## max value : 357.52902Rango del aspecto: 0 a 360 grados.
La función colorNumeric de leaflet puede tomar una paleta de RColorBrewer, pero para el aspecto, definiremos una paleta que recorra el círculo cromático (e.g., rainbow o hsv).
Usaremos ‘hcl.colors’ de R para una mejor percepción.
pal_aspect <- leaflet::colorNumeric(
palette = hcl.colors(100, palette = "Spectral", rev = TRUE), # Una paleta cíclica adecuada
domain = c(0, 360), # El aspecto va de 0 a 360 grados
na.color = "transparent"
)1. Añadir el ráster de Aspecto del Terreno
leaflet(munic2) %>%
addTiles() %>%
setView(-75.9, 2.3, 9) %>%
# 1. Añadir el ráster de Aspecto del Terreno
addRasterImage(aspect.geo, colors = pal_aspect, opacity = 0.8,
group = "Aspecto en Grados") %>%
# 2. Añadir Leyenda del Aspecto
addLegend(pal = pal_aspect, values = values(aspect.geo),
title = "Orientación del Terreno (°)") %>% # <--- ¡Aquí estaba la función que seguía al pipe!
# 3. Añadir Polígonos de Municipio
addPolygons(
color = "gray",
weight = 1.0,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.0,
# Popup (emergente) con la información zonal
popup = paste("Municipio: ", munic2$NAME_2, "<br>",
"Pendiente Media: ", round(munic2$mean_slope, 2), "%", "<br>",
"Clasificación Modal: ", munic2$class, "<br>")
)Describe el resultado y su significado.
El mapa visualiza el Aspecto (Orientación) del Terreno en grados para el Huila, utilizando una paleta de colores circular que asigna un matiz distinto a cada dirección cardinal, desde el Norte (\(0^\circ\)/\(360^\circ\)) hasta el Oeste (\(270^\circ\)). En las zonas de llanura y valle (como la cuenca del Magdalena), el aspecto es uniforme o plano, mientras que las áreas cordilleranas muestran un patrón complejo de colores variados. Este patrón segmentado y frecuente revela la presencia de numerosas laderas orientadas en distintas direcciones (Norte, Sur, Este y Oeste), lo cual es característico de un relieve escarpado y muy diseccionado por la red de drenaje.La interpretación de este resultado es fundamentalmente climática y geomorfológica. El aspecto es un determinante clave del microclima local, pues controla la cantidad de radiación solar directa que recibe la ladera, influyendo en la temperatura y la humedad del suelo. Por ejemplo, las laderas orientadas al sur son típicamente más secas y cálidas, mientras que las laderas orientadas al norte son más frescas y húmedas, un factor crítico en la distribución de la vegetación y la aptitud de los cultivos (como el café). Además, la dirección del aspecto es esencial para modelar la dirección del flujo de agua superficial y predecir zonas con mayor susceptibilidad a procesos erosivos.