Elevacion Raster Deparamento del Huila
Introducción
Requisitos
Las librerías y paquetes requeridos para la el desarrollo de este trabajo.
# INSTALL THIS PACKAGE FROM THE CONSOLE, NOT FROM THIS CHUNK
# paquetes = c("MultiscaleDTM", "exactextractr")
# install.packages(paquetes)Limpiaremos la memoria ocupada:
rm(list=ls())Ahora cargaremos las librerías requeridas:
library(elevatr)
library(sf)
library(leaflet)
library(terra)
library(MultiscaleDTM)
library(exactextractr)
library(paletteer)
library(gstat)Introduccion a DTM multiescala
Carga y ejecución de archivos
Miramos los datos contenidos en la carpeta para poder identificar cual es el archivo “.tif” que contiene los datos de elevación en ráster y queremos trabajar:
list.files("./datos")## [1] "Agua_areas_Huila.gpkg" "Depto_Hueila.gpkg" "Deptos_Colombia.gpkg"
## [4] "Elevacion_Huila.tif" "HUL_Cities.gpkg" "HUL_pop.tif"
## [7] "HUL_roads.gpkg" "Municipios_Huila.gpkg" "P1_Realpe.qgz"
## [10] "Rios_lines_Huila.gpkg" "water_areas_Huila.gpkg"Para nuestro caso el archivo deseado es “Elevacion_Huila.tif”. Tomaremos los datos-archivos geoespaciales y los procesaremos como datos de elevación y limites de los municipios:
(dem = terra::rast("datos/Elevacion_Huila.tif"))## class : SpatRaster
## dimensions : 3582, 3586, 1 (nrow, ncol, nlyr)
## resolution : 0.0006862561, 0.0006862561 (x, y)
## extent : -76.64062, -74.17971, 1.406085, 3.864255 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs
## source : Elevacion_Huila.tif
## name : Elevacion_Huila
## min value : -747
## max value : 5372Procesamiento de datos para la creación de un DEM
Después crearemos un modelo de elevación digital y reduciremos la resolución (2x2 celdas) del DEM para reducir el requerimiento de memoria:
dem2 = terra::aggregate(dem,2, "mean")## |---------|---------|---------|---------|========================================= Leeremos los datos de los municipios usando a lireria SF
(munic <- sf::st_read("datos/Municipios_Huila.gpkg"))## Reading layer `Municipios_Huila' from data source
## `C:\Users\usuario\Documents\GB2\EXAMEN_2\datos\Municipios_Huila.gpkg'
## using driver `GPKG'
## Simple feature collection with 37 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS: MAGNA-SIRGAS## Simple feature collection with 37 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## mpio_cdpmp mpio_cnmbr mpio_crslc mpio_narea
## 1 41001 NEIVA 1612 1269.81322
## 2 41006 ACEVEDO 1898 521.79729
## 3 41013 AGRADO 1837 273.73168
## 4 41016 AIPE 1912 795.63062
## 5 41020 ALGECIRAS Ordenanza 41 del 7 de Abril de 1924 589.68899
## 6 41026 ALTAMIRA 1855 180.97556
## 7 41078 BARAYA 1844 786.00202
## 8 41132 CAMPOALEGRE 1912 462.79934
## 9 41206 COLOMBIA 1886 1584.96992
## 10 41244 ELÍAS 1830 80.44174
## mpio_nano shape_Leng shape_Area geom
## 1 2023 3.0725978 0.103252186 MULTIPOLYGON (((-75.35274 3...
## 2 2023 1.1996640 0.042362075 MULTIPOLYGON (((-75.90919 1...
## 3 2023 0.9400852 0.022236902 MULTIPOLYGON (((-75.77785 2...
## 4 2023 2.1393957 0.064707602 MULTIPOLYGON (((-75.31817 3...
## 5 2023 1.2639679 0.047929546 MULTIPOLYGON (((-75.15409 2...
## 6 2023 0.6945213 0.014699156 MULTIPOLYGON (((-75.68774 2...
## 7 2023 1.4979021 0.063933424 MULTIPOLYGON (((-74.91257 3...
## 8 2023 1.0581873 0.037619594 MULTIPOLYGON (((-75.33338 2...
## 9 2023 2.3445241 0.128978175 MULTIPOLYGON (((-74.51853 3...
## 10 2023 0.5174644 0.006532031 MULTIPOLYGON (((-75.94098 2...Posteriormente tomaremos los datos de los limites de los municipios para poder recortar los datos de elevación del DEM y hacer que coincidan:
dem3 = terra::crop(dem2,munic, mask=TRUE)Transformacion de coordenadas
TOmaremos los datos de sistema de referencia espacial (CRS) del DEM y los municipios se transforman a un sistema de coordenadas proyectadas (EPSG:9377):
# revisar los parámetros de esta función de *terra* aquí:
# https://rdrr.io/github/rspatial/terra/man/project.html
(dem_plane = project(dem3, "EPSG:9377"))## class : SpatRaster
## dimensions : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : Elevacion_Huila
## min value : 315.8539
## max value : 5280.1484# revisar los parámetros de esta función de *sf* aquí:
# https://github.com/rstudio/cheatsheets/blob/main/sf.pdf
(munic_plane = sf::st_transform(munic, "EPSG:9377"))## Simple feature collection with 37 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 4596794 ymin: 1729794 xmax: 4843142 ymax: 1982829
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional
## First 10 features:
## mpio_cdpmp mpio_cnmbr mpio_crslc mpio_narea
## 1 41001 NEIVA 1612 1269.81322
## 2 41006 ACEVEDO 1898 521.79729
## 3 41013 AGRADO 1837 273.73168
## 4 41016 AIPE 1912 795.63062
## 5 41020 ALGECIRAS Ordenanza 41 del 7 de Abril de 1924 589.68899
## 6 41026 ALTAMIRA 1855 180.97556
## 7 41078 BARAYA 1844 786.00202
## 8 41132 CAMPOALEGRE 1912 462.79934
## 9 41206 COLOMBIA 1886 1584.96992
## 10 41244 ELÍAS 1830 80.44174
## mpio_nano shape_Leng shape_Area geom
## 1 2023 3.0725978 0.103252186 MULTIPOLYGON (((4738657 192...
## 2 2023 1.1996640 0.042362075 MULTIPOLYGON (((4676438 176...
## 3 2023 0.9400852 0.022236902 MULTIPOLYGON (((4691158 181...
## 4 2023 2.1393957 0.064707602 MULTIPOLYGON (((4742534 193...
## 5 2023 1.2639679 0.047929546 MULTIPOLYGON (((4760608 185...
## 6 2023 0.6945213 0.014699156 MULTIPOLYGON (((4701144 179...
## 7 2023 1.4979021 0.063933424 MULTIPOLYGON (((4787573 192...
## 8 2023 1.0581873 0.037619594 MULTIPOLYGON (((4740686 186...
## 9 2023 2.3445241 0.128978175 MULTIPOLYGON (((4831450 198...
## 10 2023 0.5174644 0.006532031 MULTIPOLYGON (((4672937 178...Cálculo de Pendiente (Slope) y Aspecto (Aspect):
# Explain what is w
# Explain what is method
# Change if needed
(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 : 1669, 1619, 2 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## names : slope, aspect
## min values : 0.00000, 1.145656e-04
## max values : 78.29077, 3.599998e+02Layer 1:
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 78.29077Cuales son los valotres de distribucion
terra::hist(slope,
main = "Huila's slope",
xlab = "Slope (in degrees)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
Layer 2:
(aspect =subset(slp_asp, 2))## class : SpatRaster
## dimensions : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : aspect
## min value : 1.145656e-04
## max value : 3.599998e+02terra::hist(aspect,
main = "Huila's aspect",
xlab = "Aspect (in degrees)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
Complementando, convertiremos la inclinacion en gradosa porcentaje de inclinacion
(slope_perc = tan(slope*(pi/180))*100)## class : SpatRaster
## dimensions : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 482.4902terra::hist(slope_perc,
main = "Huila's slope",
xlab = "Slope (in percentage)")## Warning: [hist] a sample of 37% of the cells was used (of which 71% was NA)
##Calculando la Estadistica zonal
# Reclassify the slope raster
#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)
m <- c(0, 3, 1,
3, 7, 2,
7, 12, 3,
12, 25, 4,
25, 50, 5,
50, 75, 6,
74, 490, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope_perc, m, right=TRUE)Verificaremos los datos
ahora generaremosla estadistica zonal usando exactextractr
(munic$mean_slope <- exactextractr::exact_extract(slope_perc, munic, 'mean'))## | | | 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] 23.069536 22.971500 13.887896 18.557335 32.051628 14.012997 23.509218
## [8] 18.159473 27.879963 26.955832 23.962423 23.030651 26.672739 18.523106
## [15] 31.574409 15.578115 22.617243 23.450455 27.235809 21.257706 20.056879
## [22] 18.243835 19.343857 20.165983 19.859833 22.496557 21.036348 25.655258
## [29] 32.813728 24.513773 21.133305 17.930473 19.953043 30.260111 21.859444
## [36] 7.383247 12.466011munic## Simple feature collection with 37 features and 8 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -76.62466 ymin: 1.552125 xmax: -74.41303 ymax: 3.843208
## Geodetic CRS: MAGNA-SIRGAS
## First 10 features:
## mpio_cdpmp mpio_cnmbr mpio_crslc mpio_narea
## 1 41001 NEIVA 1612 1269.81322
## 2 41006 ACEVEDO 1898 521.79729
## 3 41013 AGRADO 1837 273.73168
## 4 41016 AIPE 1912 795.63062
## 5 41020 ALGECIRAS Ordenanza 41 del 7 de Abril de 1924 589.68899
## 6 41026 ALTAMIRA 1855 180.97556
## 7 41078 BARAYA 1844 786.00202
## 8 41132 CAMPOALEGRE 1912 462.79934
## 9 41206 COLOMBIA 1886 1584.96992
## 10 41244 ELÍAS 1830 80.44174
## mpio_nano shape_Leng shape_Area geom mean_slope
## 1 2023 3.0725978 0.103252186 MULTIPOLYGON (((-75.35274 3... 23.06954
## 2 2023 1.1996640 0.042362075 MULTIPOLYGON (((-75.90919 1... 22.97150
## 3 2023 0.9400852 0.022236902 MULTIPOLYGON (((-75.77785 2... 13.88790
## 4 2023 2.1393957 0.064707602 MULTIPOLYGON (((-75.31817 3... 18.55733
## 5 2023 1.2639679 0.047929546 MULTIPOLYGON (((-75.15409 2... 32.05163
## 6 2023 0.6945213 0.014699156 MULTIPOLYGON (((-75.68774 2... 14.01300
## 7 2023 1.4979021 0.063933424 MULTIPOLYGON (((-74.91257 3... 23.50922
## 8 2023 1.0581873 0.037619594 MULTIPOLYGON (((-75.33338 2... 18.15947
## 9 2023 2.3445241 0.128978175 MULTIPOLYGON (((-74.51853 3... 27.87996
## 10 2023 0.5174644 0.006532031 MULTIPOLYGON (((-75.94098 2... 26.95583hist(munic$mean_slope,
main = "Huila's mean slope",
xlab = "Slope (in percentage)")
(munic$class <- exactextractr::exact_extract(rc, munic, 'mode'))## | | | 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] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4hist(munic$class,
main = "Cordoba's reclassified slope",
xlab = "Slope (as a category)")
convertiremosen coordenadas geograficas
# This is the reclassified slope
(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7# This is the percentage slope
(slope.geo = project(slope_perc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.0000
## max value : 452.9977Graficaremos los datos de elevacion
# expand the number of colors to improve your visualization
# https://r-graph-gallery.com/42-colors-names.html
palette<- colorNumeric(c("#C1766FFF", "#A65461FF", "#541F3FFF"), values(slope.geo),
na.color = "transparent")leaflet(munic) %>% addTiles() %>% setView(-75.3, 2.7, 7) %>%
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 = palette, opacity = 0.8) %>%
addLegend(pal = palette, values = values(slope.geo),
title = "Terrain slope in Huila (%)")## 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'Visualizcion de la grafica pero en grados
(slope = subset(slp_asp, 1))## class : SpatRaster
## dimensions : 1669, 1619, 1 (nrow, ncol, nlyr)
## resolution : 152.2093, 152.2093 (x, y)
## extent : 4596760, 4843187, 1729543, 1983581 (xmin, xmax, ymin, ymax)
## coord. ref. : MAGNA-SIRGAS 2018 / Origen-Nacional (EPSG:9377)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 78.29077# Reclassify the slope raster
#rc <- classify(slope_perc, c(0, 3, 7, 12, 25,50, 75), include.lowest=TRUE, brackets=TRUE)
m <- c(0, 1.7, 1,
1.7, 4, 2,
4, 6.9, 3,
6.9, 14.3, 4,
14.3, 28.6, 5,
28.6, 43, 6,
43, 79, 7)
m <- matrix(m, ncol=3, byrow = TRUE)
rc <- classify(slope, m, right=TRUE)(munic$mean_slope2 <- exactextractr::exact_extract(slope, munic, '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] 12.644789 12.797388 7.769434 10.266978 17.480240 7.859426 12.981675
## [8] 9.856176 15.294997 14.747098 13.205628 12.670143 14.697580 10.204793
## [15] 17.106588 8.650787 12.526482 12.855496 15.001105 11.832603 11.193542
## [22] 10.081303 10.848270 11.223944 11.026142 11.898920 11.653651 14.111075
## [29] 17.829418 13.594316 11.653077 9.940928 10.950447 16.399023 12.131503
## [36] 4.158818 7.014504(munic$class2 <- exactextractr::exact_extract(rc, munic, '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] 5 4 4 4 5 4 5 1 5 5 5 4 5 5 5 4 4 5 4 4 4 4 4 4 5 2 4 5 5 5 4 4 5 5 4 1 4(rc.geo = project(rc, "EPSG:4326"))## class : SpatRaster
## dimensions : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0
## max value : 7(slope.geo2 = project(slope, "EPSG:4326"))## class : SpatRaster
## dimensions : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : slope
## min value : 0.00000
## max value : 77.37261palredgreen <- colorNumeric(c("darkseagreen3","yellow2", "orange", "brown2", "darkred"), values(slope.geo2),
na.color = "transparent")leaflet(munic) %>% addTiles() %>% setView(-75.3, 2.7, 7) %>%
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$class2, "<br>")) %>%
addRasterImage(slope.geo2, colors = palredgreen, opacity = 0.8) %>%
addLegend(pal = palredgreen, values = values(slope.geo2),
title = "Terrain slope in Huila (grados)")## 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'Visualizacion de la rientacion de los terrenos
mapa 1
# Visualizar la orientación de los terrenos
# Usaremos la capa 'aspect' ya obtenida previamente
# Convertimos la capa de orientación a coordenadas geográficas
(aspect.geo = project(aspect, "EPSG:4326"))## class : SpatRaster
## dimensions : 1677, 1618, 1 (nrow, ncol, nlyr)
## resolution : 0.001372493, 0.001372493 (x, y)
## extent : -76.63096, -74.41026, 1.548552, 3.850222 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source(s) : memory
## name : aspect
## min value : 0.4234968
## max value : 359.7618713# Crear una paleta de colores para representar los 360 grados de orientación
pal_aspect <- colorNumeric(palette = "Spectral", domain = values(aspect.geo), na.color = "transparent")# Crear el mapa interactivo con Leaflet
leaflet(munic) %>%
addTiles() %>%
setView(-75.3, 2.7, 7) %>% # Centrar el mapa en Huila
addPolygons(
color = "gray",
weight = 1.0,
smoothFactor = 0.5,
opacity = 0.4,
fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr)
) %>%
addRasterImage(aspect.geo, colors = pal_aspect, opacity = 0.8) %>%
addLegend(pal = pal_aspect, values = values(aspect.geo),
title = "Orientación de los terrenos (grados)")## 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'mapa 2
# Reclasificación del aspecto en 8 direcciones cardinales
m_aspect <- c(0, 22.5, 1, # Norte
22.5, 67.5, 2, # Noreste
67.5, 112.5, 3, # Este
112.5, 157.5, 4, # Sureste
157.5, 202.5, 5, # Sur
202.5, 247.5, 6, # Suroeste
247.5, 292.5, 7, # Oeste
292.5, 337.5, 8, # Noroeste
337.5, 360, 1) # Norte (cerrando el ciclo)
# Crear el mapa reclasificado
m_aspect <- matrix(m_aspect, ncol=3, byrow=TRUE)
aspect_reclass <- classify(aspect, m_aspect, right=TRUE)# Colores para las diferentes direcciones cardinales
palette_aspect <- colorFactor(c("blue", "green", "yellow", "orange", "red", "purple", "brown", "pink"),
values(aspect_reclass), na.color = "transparent")# Visualización en leaflet
leaflet(munic) %>% addTiles() %>% setView(-75.3, 2.7, 7) %>%
addPolygons(color = "gray", weight = 1.0, smoothFactor = 0.5,
opacity = 0.4, fillOpacity = 0.10,
popup = paste("Municipio: ", munic$mpio_cnmbr, "<br>",
"Aspecto: ", munic$class, "<br>")) %>%
addRasterImage(aspect_reclass, colors = palette_aspect, opacity = 0.8) %>%
addLegend(pal = palette_aspect, values = values(aspect_reclass),
title = "Orientación de los terrenos")## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA## 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'