Procesamiento de datos de elevación para el departamento del Guaviare

Introducción

Este cuaderno tiene como objeto ilustrar diferentes funciones para obtener, procesar y visualizar modelos de elevación digital (DEM) en R.

Configuración

En esta sección escribimos código para cargar todas las librerias que necesitamos. Puesto que era una cantidad considerable de paquetes que se debia instalar, se guardan en un vector llamado “paquetes” para que de este modo al utilizar la función “install.packages”, sea posible descargarlos todos al mismo tiempo.

# Algunos paquetes a instalar son paquetes = c("rgdal","raster", "rgeos", terra","elevatr","rasterVis", "rgl", "mapview")

paquetes<-c("rgdal","raster", "rgeos", "terra","elevatr","rasterVis", "rgl", "mapview")

library(tidyverse)
library(rgdal)
library(raster)
library(rgeos)
library(terra)
library(elevatr)
library(rasterVis)
library(rgl)
library(mapview)
library(sf)
library(leaflet)
library(tmap)

Antes de continuar, y con el fin de facilitarle el trabajo al computador (y también a nosotros), se recomienda limpiar la memoria de trabajo…

rm(list = ls())

Datos de elevación tipo raster requerida para el departamento

Primero, se debe descargar un shapefile que represente los municipios de nuestro departamento. Cabe resaltar que el documento a leer tiene información geográfica de todos los departamentos

Para ello podemos utilizar una función del paquete sf

muni<-st_read("C:\\Users\\María Isabel\\Desktop\\GBI\\GB2\\cuaderno_dem\\ADMINISTRATIVO\\MGN_MPIO_POLITICO.shp")
## Reading layer `MGN_MPIO_POLITICO' from data source 
##   `C:\Users\María Isabel\Desktop\GBI\GB2\cuaderno_dem\ADMINISTRATIVO\MGN_MPIO_POLITICO.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1121 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.73562 ymin: -4.229406 xmax: -66.84722 ymax: 13.39473
## Geodetic CRS:  MAGNA-SIRGAS

Puesto que únicamente estamos interesados en los municipios del guaviare, filtramos la información

muni_gua<-dplyr::filter(muni,DPTO_CNMBR == "GUAVIARE")

muni_gua
## Simple feature collection with 4 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.66391 ymin: 0.6554126 xmax: -69.99511 ymax: 2.924987
## Geodetic CRS:  MAGNA-SIRGAS
##   DPTO_CCDGO MPIO_CCDGO MPIO_CDPMP DPTO_CNMBR            MPIO_CNMBR
## 1         95        001      95001   GUAVIARE SAN JOSÉ DEL GUAVIARE
## 2         95        015      95015   GUAVIARE               CALAMAR
## 3         95        025      95025   GUAVIARE            EL RETORNO
## 4         95        200      95200   GUAVIARE            MIRAFLORES
##                                MPIO_CRSLC MPIO_NAREA MPIO_CSMBL MPIO_VGNC
## 1 Decreto Nal 1165 del 7 de Junio de 1976   16757.07          4      2021
## 2   Ordenanza 001 del 7 de Agosto de 1992   13069.59          4      2021
## 3   Ordenanza 001 del 7 de Agosto de 1992   12909.12          4      2021
## 4   Ordenanza 001 del 7 de Agosto de 1992   12808.08          4      2021
##   MPIO_TIPO Shape_Leng Shape_Area                       geometry
## 1 MUNICIPIO  14.085782   1.360391 MULTIPOLYGON (((-71.31266 2...
## 2 MUNICIPIO   7.529464   1.061812 MULTIPOLYGON (((-72.6716 2....
## 3 MUNICIPIO  12.315742   1.047240 MULTIPOLYGON (((-72.73135 2...
## 4 MUNICIPIO   7.385424   1.039462 MULTIPOLYGON (((-72.01357 1...
head(muni_gua)
## Simple feature collection with 4 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.66391 ymin: 0.6554126 xmax: -69.99511 ymax: 2.924987
## Geodetic CRS:  MAGNA-SIRGAS
##   DPTO_CCDGO MPIO_CCDGO MPIO_CDPMP DPTO_CNMBR            MPIO_CNMBR
## 1         95        001      95001   GUAVIARE SAN JOSÉ DEL GUAVIARE
## 2         95        015      95015   GUAVIARE               CALAMAR
## 3         95        025      95025   GUAVIARE            EL RETORNO
## 4         95        200      95200   GUAVIARE            MIRAFLORES
##                                MPIO_CRSLC MPIO_NAREA MPIO_CSMBL MPIO_VGNC
## 1 Decreto Nal 1165 del 7 de Junio de 1976   16757.07          4      2021
## 2   Ordenanza 001 del 7 de Agosto de 1992   13069.59          4      2021
## 3   Ordenanza 001 del 7 de Agosto de 1992   12909.12          4      2021
## 4   Ordenanza 001 del 7 de Agosto de 1992   12808.08          4      2021
##   MPIO_TIPO Shape_Leng Shape_Area                       geometry
## 1 MUNICIPIO  14.085782   1.360391 MULTIPOLYGON (((-71.31266 2...
## 2 MUNICIPIO   7.529464   1.061812 MULTIPOLYGON (((-72.6716 2....
## 3 MUNICIPIO  12.315742   1.047240 MULTIPOLYGON (((-72.73135 2...
## 4 MUNICIPIO   7.385424   1.039462 MULTIPOLYGON (((-72.01357 1...

Debido a que el CRS que maneja el archivo no corresponde con el que normalmente se trabaja que es el WGS 84 (y que es imprescindible para utilizar ciertas funciones) se procede a convertir. Para ello primero se crea el objeto con el CRS deseado

new_crs<-4326

Luego se le asigna dicho CRS a los municipios

muni_guav_new<-st_transform(muni_gua,new_crs)

muni_guav_new
## Simple feature collection with 4 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -73.66391 ymin: 0.6554126 xmax: -69.99511 ymax: 2.924987
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO MPIO_CDPMP DPTO_CNMBR            MPIO_CNMBR
## 1         95        001      95001   GUAVIARE SAN JOSÉ DEL GUAVIARE
## 2         95        015      95015   GUAVIARE               CALAMAR
## 3         95        025      95025   GUAVIARE            EL RETORNO
## 4         95        200      95200   GUAVIARE            MIRAFLORES
##                                MPIO_CRSLC MPIO_NAREA MPIO_CSMBL MPIO_VGNC
## 1 Decreto Nal 1165 del 7 de Junio de 1976   16757.07          4      2021
## 2   Ordenanza 001 del 7 de Agosto de 1992   13069.59          4      2021
## 3   Ordenanza 001 del 7 de Agosto de 1992   12909.12          4      2021
## 4   Ordenanza 001 del 7 de Agosto de 1992   12808.08          4      2021
##   MPIO_TIPO Shape_Leng Shape_Area                       geometry
## 1 MUNICIPIO  14.085782   1.360391 MULTIPOLYGON (((-71.31266 2...
## 2 MUNICIPIO   7.529464   1.061812 MULTIPOLYGON (((-72.6716 2....
## 3 MUNICIPIO  12.315742   1.047240 MULTIPOLYGON (((-72.73135 2...
## 4 MUNICIPIO   7.385424   1.039462 MULTIPOLYGON (((-72.01357 1...

Ahora debemos crear un nuevo objeto con los centroides de los municipios. El objetivo es encontrar un punto central para cada región, para esto primero convertimos simple features en polígonos espaciales

sp_muni_guav<-as_Spatial(muni_guav_new)

Luego usamos la función gCentroid del paquete rgeos

centers<-data.frame(gCentroid(sp_muni_guav,byid = TRUE))

centers$region<-row.names(sp_muni_guav)
centers$muni_guav_new<-sp_muni_guav$MPIO_CNMBR

centers
##           x        y region         muni_guav_new
## 1 -71.91798 2.484494      1 SAN JOSÉ DEL GUAVIARE
## 2 -73.03837 1.596191      2               CALAMAR
## 3 -71.58570 2.086327      3            EL RETORNO
## 4 -72.01811 1.363975      4            MIRAFLORES

Ahora debemos descargar datos de elevación para el Guaviare

elevation<-get_elev_raster(muni_guav_new,z=10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.
elevation
## class      : RasterLayer 
## dimensions : 4094, 5633, 23061502  (nrow, ncol, ncell)
## resolution : 0.00068651, 0.00068651  (x, y)
## extent     : -73.82812, -69.96101, 0.3518836, 3.162456  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : file297418e35aa6.tif 
## names      : file297418e35aa6

Nota 1: Z es el zoom de los datos, es decir, entre mayor sea su valor, mayor resolución espacial tendrá Nota 2: Los datos de elevación se encuentran en metros

Cabe resaltar que dicho Raster, es solo temporal, en ese caso, si se desea guardar en el directorio local de trabajo, es posible hacerlo con la función writeRaster de la libreria rgdal

if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("C:\\Users\\María Isabel\\Desktop\\GBI\\GB2\\cuaderno_dem"), format="GTiff", overwrite=TRUE)
}

Visualización de los datos de elevación

Para propósitos de visualización, una baja resolución podría optimizar dicha tarea. Esto se hace únicamente si los datos son muy pesados

elevation2 <- aggregate(elevation, fact=4, fun=mean)

Es imprescindible fijar una paleta de colores adecuada

pal<-colorNumeric(c("forestgreen","yellow","tan","brown"),values(elevation2),na.color="transparent")

Ahora utilizaremos la libreria leaflet para visualizar los datos de elevación

leaflet(muni_guav_new) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", muni_guav_new$MPIO_CNMBR, "<br>",
                          "Km2: ", muni_guav_new$MPIO_NAREA, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(elevation2, colors = pal, opacity = 0.9)  %>%
  addLegend(pal = pal, values = values(elevation2),
    title = "Datos de elevación para el Guaviare (mts)")

Reproyectando los datos de elevación

Puesto que los datos de coordenadas geográficas se encuentran en grados pero la elevación en metros, es imprescindible llevar a cabo la conversión a coordenadas planas en metros para que de este modo se puedan utilizar las unidades adecuadas. Para esto, en primer lugar se define un CRS objetivo

# WGS 84 a UTM zona 18 N

spatialref<-CRS(SRS_string = "EPSG:32618")

Luego se reproyecta el DEM usando la función projectRaster del paquete raster

# Primero se crea un objeto raster usando la función projectExtent

pr3<-projectExtent(elevation,crs = spatialref)

# Ajustamos el tamaño de la celda

res(pr3)<-80

# Proyectamos

rep_elev<- projectRaster(elevation,pr3)

rep_elev
## class      : RasterLayer 
## dimensions : 3900, 5390, 21021000  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-29_171046_10612_03454.grd 
## names      : layer 
## values     : -486.0538, 883.0737  (min, max)

Ahora, reproyectamos la tabla de polígonos espaciales, en donde se representan los municipios de nuestro departamento

rep_munic<-spTransform(sp_muni_guav, spatialref)
rep_munic
## class       : SpatialPolygonsDataFrame 
## features    : 4 
## extent      : 648629.8, 1057215, 72539.32, 323974  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 12
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CDPMP, DPTO_CNMBR,            MPIO_CNMBR,                              MPIO_CRSLC,    MPIO_NAREA, MPIO_CSMBL, MPIO_VGNC, MPIO_TIPO,    Shape_Leng,    Shape_Area 
## min values  :         95,        001,      95001,   GUAVIARE,               CALAMAR, Decreto Nal 1165 del 7 de Junio de 1976, 12808.0770572,          4,      2021, MUNICIPIO, 7.38542439659, 1.03946228704 
## max values  :         95,        200,      95200,   GUAVIARE, SAN JOSÉ DEL GUAVIARE,   Ordenanza 001 del 7 de Agosto de 1992, 16757.0657353,          4,      2021, MUNICIPIO, 14.0857819173, 1.36039119785

Guardemos el DEM creado.

writeRaster(rep_elev,filename = "./rep_guav_DEM.tif",datatype = 'INT4S', overwrite = TRUE)

Estadísticas básicas de los datos de elevación

Con el fin de eliminar datos poco significativos, se utiliza el siguiente código para eliminar elevaciones por debajo de los 0 m

rep_elev[rep_elev < 0.0] <-0.0

Ahora podemos crear un histograma que nos permita visualizar la información de elevación

hist(rep_elev, xlab = "Elevación (metros)", main = "Elevación para el departamento del Guaviare", col = "darkgreen")

Se observa inicialmente que las elevaciones del departamento se concentran principalmente entre los 100 a 300 metros.

Por otro lado, se puede obtener un resumen de estadísticas para el DEM

promedio<-cellStats(rep_elev, 'mean')
minimo<-cellStats(rep_elev,'min')
maximo<-cellStats(rep_elev,'max')
desviacion<-cellStats(rep_elev, 'sd')

# Vector único de estadísticas

metricas<-c('mean','min','max','sd')

valores<-c(promedio,minimo,maximo,desviacion)

# Tabla con el resultado

(df_estadisticas<-data.frame(metricas,valores))
##   metricas   valores
## 1     mean 226.62019
## 2      min   0.00000
## 3      max 883.07373
## 4       sd  50.27808

Obtención de variables geomorfométricas

Particularmente para la pendiente y el aspecto, este último, también conocido como orientación

Para iniciar primero se obtienen los datos relacionados con la información que deseo visualizar en un mapa

(pendiente<-terrain(rep_elev,opt='slope', unit='degrees'))
## class      : RasterLayer 
## dimensions : 3900, 5390, 21021000  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 74.66911  (min, max)
(orientacion<-terrain(rep_elev,opt='aspect', unit='degrees'))
## class      : RasterLayer 
## dimensions : 3900, 5390, 21021000  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : aspect 
## values     : 0, 360  (min, max)

Pendiente

Para poder observar la pendiente, es una buena idea en primer lugar, reducir el tamaño de la capa raster

pendiente2<-aggregate(pendiente, fact = 4, fun = mean)

Ahora se proyecta el conjunto de datos de elevación

pendiente3<- projectRasterForLeaflet(pendiente2,"ngb")

Se crea una paleta de colores

pal2<-colorNumeric(c("#537188", "#CBB279", "#E1D4BB", "#EEEEEE"), values(pendiente3), na.color = "transparent")

Ahora si se puede observar el mapa de pendientes

leaflet(muni_guav_new) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", muni_guav_new$MPIO_CNMBR, "<br>",
                          "Km2: ", muni_guav_new$MPIO_NAREA, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(pendiente3, colors = pal2, opacity = 1.0)  %>%
  addLegend(pal = pal2, values = values(pendiente3),
    title = "Pendiente para el Guaviare (%)")

La clasificación anterior fue automática, no obstante, el objetivo es tener algún punto de referencia para podernos guiar mejor, particularmente trabajaremos con la clasificación dada por el Instituto Geográfico Colombiano Agustin Codazzi (IGAC), en donde se puede observar siete categorías (solo trabajaremos con 6). En ese sentido se hace necesario reclasificar. Para esto, se desea trabajar con la libreria terra puesto que esta es más reciente y rápida. Para poder utilizarla se debe convertir el objeto raster a un objeto terra

(nspendiente2<- as(pendiente2, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 975, 1348, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 630212.2, 1061572, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :      slope 
## min value   :  0.1428026 
## max value   : 61.4023728

Hecho lo anterior, podemos iniciar con la reclasificación de pendiente

# Creamos un vector con las pendientes

m<-c(0, 3, 1, 3, 7, 2, 7, 12, 3, 12,
     25, 4, 25, 50, 5, 50, 75, 6)

# Creamos la matriz de reclasificación 

(rclmat<-matrix(m, ncol = 3, byrow = TRUE))
##      [,1] [,2] [,3]
## [1,]    0    3    1
## [2,]    3    7    2
## [3,]    7   12    3
## [4,]   12   25    4
## [5,]   25   50    5
## [6,]   50   75    6
# LLevamos a cabo la reclasificación

pendiente_rec<-terra::classify(nspendiente2, rclmat, right = TRUE)

pendiente_rec
## class       : SpatRaster 
## dimensions  : 975, 1348, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 630212.2, 1061572, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        : slope 
## min value   :     1 
## max value   :     6

Visualización estática

Para que el mapa de pendiente quede estático, se utilizará las funcionalidades de tmap

mapa.pendiente <- tm_shape(pendiente_rec) +
  tm_raster(alpha = 0.6, style= "cat", 
            title="Pendiente") +  
  # Borders only
  tm_shape(muni_guav_new) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.3, alpha=0.7) +
  tm_scale_bar(text.size = 0.5) +
  tm_compass(position = c("right", "top"), size = 2) +
  tm_graticules(alpha = 0.2) +
  tm_layout(outer.margins = 0.01, legend.position = c("left", "top"), title= 'Slope classes', title.position = c('right', 'top')) +
  tm_credits("Data source: Mapzen & DANE", size=0.3)  

mapa.pendiente
## stars object downsampled to 1175 by 850 cells. See tm_shape manual (argument raster.downsample)

Guardemos el mapa para futura consulta

tmap_save(mapa.pendiente, "./mapa_pendiente_guaviare.png", width = 8.5, height = 11, dpi = 600)
## stars object downsampled to 1175 by 850 cells. See tm_shape manual (argument raster.downsample)
## Map saved to C:\Users\María Isabel\Desktop\GBI\GB2\cuaderno_dem\mapa_pendiente_guaviare.png
## Resolution: 5100 by 6600 pixels
## Size: 8.5 by 11 inches (600 dpi)

Visualización interactiva por departamento

Para ello utilizaremos la información espacial de los municipios. En primer lugar debemos convertir el objeto sf en un SpatVector

(nmunic<-as(rep_munic, "SpatVector"))
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 4, 12  (geometries, attributes)
##  extent      : 648629.8, 1057215, 72539.32, 323974  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       : DPTO_CCDGO MPIO_CCDGO MPIO_CDPMP DPTO_CNMBR      MPIO_CNMBR
##  type        :      <chr>      <chr>      <chr>      <chr>           <chr>
##  values      :         95        001      95001   GUAVIARE SAN JOSÉ DEL G~
##                        95        015      95015   GUAVIARE         CALAMAR
##                        95        025      95025   GUAVIARE      EL RETORNO
##       MPIO_CRSLC MPIO_NAREA MPIO_CSMBL MPIO_VGNC MPIO_TIPO Shape_Leng
##            <chr>      <num>      <chr>     <int>     <chr>      <num>
##  Decreto Nal 11~  1.676e+04          4      2021 MUNICIPIO      14.09
##  Ordenanza 001 ~  1.307e+04          4      2021 MUNICIPIO      7.529
##  Ordenanza 001 ~  1.291e+04          4      2021 MUNICIPIO      12.32
##  Shape_Area
##       <num>
##        1.36
##       1.062
##       1.047

Ahora convertimos la pendiente que está en formato raster a terra

(nspendiente<-as(pendiente, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 3900, 5390, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 74.66911

Luego calculamos las categorías acumuladas para cada municipio (promedio)

pendiente.municipio<- zonal(nspendiente, nmunic, fun = mean, as.raster = TRUE)

pendiente.municipio
## class       : SpatRaster 
## dimensions  : 3900, 5390, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :    slope 
## min value   : 1.497147 
## max value   : 1.992928

Es hora de visualizar el mapa interactivo

tmap_mode("view")
## tmap mode set to interactive viewing
  tm_shape(muni_guav_new) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(pendiente.municipio) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Pendiente media")+
tm_layout(legend.outside = TRUE)
## stars object downsampled to 1176 by 851 cells. See tm_shape manual (argument raster.downsample)

Aspecto

Cabe recordar que el aspecto es la orientación que tienen las pendientes en un terreno dado. Su visualización en un mapa se obtiene al ilustrar las orientaciones de las pendientes desde 0 (siendo este el norte) hasta 360 grados (nuevamente llegando al norte), en sentido de las manecillas del reloj. Los diferentes valores de aspecto se asocian a colores o tonos determinados.

Para poder observar la orientación, es una buena idea reducir el tamaño de la capa raster

orientacion2<-aggregate(orientacion, fact = 4, fun = mean)

Ahora se proyecta el conjunto de datos de elevación

orientacion3<- projectRasterForLeaflet(orientacion2,"ngb")

Se crea una paleta de colores

pal3<-colorNumeric(c("#FFB84C", "#F266AB","#A459D1", "#2CD3E1"), values(orientacion3), na.color = "transparent")

Ahora si se puede observar el mapa de orientación

leaflet(muni_guav_new) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", muni_guav_new$MPIO_CNMBR, "<br>",
                          "Km2: ", muni_guav_new$MPIO_NAREA, "<br>")) %>%
  addLabelOnlyMarkers(data = centers,
                    lng = ~x, lat = ~y, label = ~region,
                    labelOptions = labelOptions(noHide = TRUE, direction = 'top', textOnly = TRUE)) %>%
  addRasterImage(orientacion3, colors = pal3, opacity = 1.0)  %>%
  addLegend(pal = pal3, values = values(orientacion3),
    title = "Orientación para el Guaviare")

A nivel general, se observa una tendencia de valores entre los 150º a 200º, esto quiere decir que en su mayoría las pendientes están orientadas hacia el suroriente

Visualización interactiva por departamento

Ahora convertimos el aspecto que está en formato raster a terra

(nsorientacion<-as(orientacion, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 3900, 5390, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        : aspect 
## min value   :      0 
## max value   :    360

Luego calculamos las categorías acumuladas para cada municipio (promedio)

orientacion.municipio<- zonal(nsorientacion, nmunic, fun = mean, as.raster = TRUE)

orientacion.municipio
## class       : SpatRaster 
## dimensions  : 3900, 5390, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 630212.2, 1061412, 38912.52, 350912.5  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :   aspect 
## min value   : 176.1057 
## max value   : 178.1265

Es hora de visualizar el mapa interactivo

tmap_mode("view")
## tmap mode set to interactive viewing
  tm_shape(muni_guav_new) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(orientacion.municipio) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Orientacion media")+
tm_layout(legend.outside = TRUE)
## stars object downsampled to 1176 by 851 cells. See tm_shape manual (argument raster.downsample)

Teniendo en cuenta los valores medios de orientación por departamento y el rango en el cual se observan, entonces si se podría hablar de una tendencia de orientación de las pendientes hacia el suroriente.

Bibliografía

Lizarazo, I., 2023. Elevation data processing and analysis in R. Available at https://rpubs.com/ials2un/dem_analysis

sessionInfo()
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tmap_3.3-3       leaflet_2.1.2    sf_1.0-11        mapview_2.11.0  
##  [5] rgl_1.1.3        rasterVis_0.51.5 lattice_0.20-45  elevatr_0.4.2   
##  [9] terra_1.7-29     rgeos_0.6-3      raster_3.6-20    rgdal_1.6-6     
## [13] sp_1.6-0         lubridate_1.9.2  forcats_1.0.0    stringr_1.5.0   
## [17] dplyr_1.1.0      purrr_1.0.1      readr_2.1.4      tidyr_1.3.0     
## [21] tibble_3.1.8     ggplot2_3.4.1    tidyverse_2.0.0 
## 
## loaded via a namespace (and not attached):
##  [1] satellite_1.0.4         httr_1.4.5              progress_1.2.2         
##  [4] webshot_0.5.4           RColorBrewer_1.1-3      tools_4.2.2            
##  [7] bslib_0.4.2             utf8_1.2.3              R6_2.5.1               
## [10] KernSmooth_2.23-20      DBI_1.1.3               colorspace_2.1-0       
## [13] withr_2.5.0             prettyunits_1.1.1       tidyselect_1.2.0       
## [16] curl_5.0.0              compiler_4.2.2          progressr_0.13.0       
## [19] leafem_0.2.0            cli_3.6.0               sass_0.4.5             
## [22] scales_1.2.1            classInt_0.4-9          hexbin_1.28.3          
## [25] proxy_0.4-27            digest_0.6.31           rmarkdown_2.20         
## [28] dichromat_2.0-0.1       base64enc_0.1-3         jpeg_0.1-10            
## [31] pkgconfig_2.0.3         htmltools_0.5.4         highr_0.10             
## [34] fastmap_1.1.1           htmlwidgets_1.6.2       rlang_1.0.6            
## [37] rstudioapi_0.14         farver_2.1.1            jquerylib_0.1.4        
## [40] generics_0.1.3          zoo_1.8-12              jsonlite_1.8.4         
## [43] crosstalk_1.2.0         magrittr_2.0.3          s2_1.1.2               
## [46] interp_1.1-4            Rcpp_1.0.10             munsell_0.5.0          
## [49] fansi_1.0.4             abind_1.4-5             lifecycle_1.0.3        
## [52] stringi_1.7.12          leafsync_0.1.0          yaml_2.3.7             
## [55] tmaptools_3.1-1         grid_4.2.2              parallel_4.2.2         
## [58] crayon_1.5.2            deldir_1.0-9            stars_0.6-1            
## [61] hms_1.1.2               knitr_1.42              pillar_1.8.1           
## [64] markdown_1.5            codetools_0.2-18        stats4_4.2.2           
## [67] wk_0.7.2                XML_3.99-0.14           glue_1.6.2             
## [70] evaluate_0.20           leaflet.providers_1.9.0 latticeExtra_0.6-30    
## [73] png_0.1-8               vctrs_0.5.2             tzdb_0.3.0             
## [76] gtable_0.3.1            slippymath_0.3.1        cachem_1.0.7           
## [79] xfun_0.37               mime_0.12               lwgeom_0.2-11          
## [82] e1071_1.7-13            class_7.3-20            viridisLite_0.4.1      
## [85] units_0.8-1             timechange_0.2.0        ellipsis_0.3.2