Procesamiento y análisis de datos de elevación en R

¿Por qué este cuaderno?

Este cuaderno ilustra varias funcionalidades para obtener, procesar y visualizar modelos digitales de elevación (DEM) en R. Tiene como objetivo ayudar a los estudiantes de Geomática Básica de la Universidad Nacional a familiarizarse con los DEM y las variables geomorfométricas.

Se usaron las siguientes librerias:

rgdal, raster, rgeos, terra, elevat, rasterVis, rgl, mapview Posteriormente se limpió la memoria y se cargaron las librerias necesarias

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

Introducción

La finalidad de los datos de elevación obtenidos es dar sentido a la aplicación de una amplia gama de aplicaciones, que incluyen, por ejemplo, visualización, hidrología y modelado ecológico. Entonces para poder hacer uso de la herramienta denominada DEM el cual a partir de datos de elevacion tipo raster se usa el paquete elevatr que toma los mosaicos de terreno de Amazon Web Services.

Obtención de datos de elevación ráster para el departamento Bolivar

Primero, se va a cargar un shapefile que represente los municipios de Bolivar. En este cuaderno utilizaré un shapefile descargado del geoportal del DANE.

Leamos el archivo de forma usando una función provista por el paquete sf:

munic <-  st_read("C:/Cuaderno3R/Bolivar/BOLIVAR_MPIO.shp")

## Reading layer `BOLIVAR_MPIO' from data source 
##   `C:\Cuaderno3R\BOLIVAR\BOLIVAR_MPIO.shp' using driver `ESRI Shapefile'
## Simple feature collection with 46 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.19063 ymin: 6.99916 xmax: -73.74578 ymax: 10.80147
## Geodetic CRS:  WGS 84

Ahora se verifica la geometría, el cuadro delimitador, el CRS y los atributos del objeto munic:

head(munic)

## Simple feature collection with 6 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.19063 ymin: 8.15458 xmax: -73.8718 ymax: 10.70319
## Geodetic CRS:  WGS 84
##   DPTO_CCDGO MPIO_CCDGO          MPIO_CNMBR
## 1         13      13001 CARTAGENA DE INDIAS
## 2         13      13006                ACHÍ
## 3         13      13030   ALTOS DEL ROSARIO
## 4         13      13042              ARENAL
## 5         13      13052              ARJONA
## 6         13      13062         ARROYOHONDO
##                             MPIO_CRSLC MPIO_NAREA MPIO_NANO DPTO_CNMBR
## 1                                 1599   593.2364      2017    BOLÍVAR
## 2     Ordenanza 15 de Abril 18 de 1934   951.7664      2017    BOLÍVAR
## 3 Ordenanza 30 de diciembre 13 de 1994   303.7462      2017    BOLÍVAR
## 4            ORD 18 DE MAYO 16 DE 1996   461.9323      2017    BOLÍVAR
## 5                                 1770   589.5791      2017    BOLÍVAR
## 6       ORD 41 DE DICIEMBRE 02 DE 1997   163.8084      2017    BOLÍVAR
##   Shape_Leng Shape_Area                       geometry
## 1  4.0370005 0.04912313 MULTIPOLYGON (((-75.29333 1...
## 2  1.8371512 0.07817905 MULTIPOLYGON (((-74.53501 8...
## 3  0.9849827 0.02495982 MULTIPOLYGON (((-74.11309 8...
## 4  1.5548421 0.03791813 MULTIPOLYGON (((-73.94174 8...
## 5  1.3582253 0.04861694 MULTIPOLYGON (((-75.30055 1...
## 6  0.7833122 0.01351407 MULTIPOLYGON (((-74.98125 1...

Ahora se da la creación de nuevo objeto con centroides de municipios: Posteriormente con la fialidad de encontrar un punto central para cada región y convertir las características simples en polígonos espaciales

sp.munic = as_Spatial(munic)

Entonces usamos la función gCentroid del paquete rgeos para cumplir con el primer objetivo

centers <- data.frame(gCentroid(sp.munic, byid = TRUE))
centers$region <- row.names(sp.munic)
centers$munic <- sp.munic$MPIO_CNMBR

centers

##            x         y region                 munic
## 1  -75.46828 10.452059      1   CARTAGENA DE INDIAS
## 2  -74.47806  8.622943      2                  ACHÍ
## 3  -74.19516  8.757469      3     ALTOS DEL ROSARIO
## 4  -74.10696  8.341018      4                ARENAL
## 5  -75.37122 10.164221      5                ARJONA
## 6  -75.02690 10.241718      6           ARROYOHONDO
## 7  -74.16271  8.836185      7      BARRANCO DE LOBA
## 8  -75.00119 10.218828      8               CALAMAR
## 9  -74.10014  7.239229      9            CANTAGALLO
## 10 -74.68841  9.235269     10                CICUCO
## 11 -74.89446  9.540048     11               CÓRDOBA
## 12 -75.31447 10.556637     12             CLEMENCIA
## 13 -75.15794  9.699820     13  EL CARMEN DE BOLÍVAR
## 14 -74.92811 10.034837     14              EL GUAMO
## 15 -73.91230  8.864819     15              EL PEÑON
## 16 -74.14707  8.983829     16       HATILLO DE LOBA
## 17 -74.74730  9.169513     17              MAGANGUÉ
## 18 -75.17978 10.178214     18               MAHATES
## 19 -74.22584  9.054590     19             MARGARITA
## 20 -75.34380  9.976017     20         MARIA LA BAJA
## 21 -74.46343  8.017100     21           MONTECRISTO
## 22 -74.54520  9.142511     22                MOMPÓS
## 23 -73.97135  8.253700     23               MORALES
## 24 -74.11589  8.506135     24                NOROSI
## 25 -74.38588  8.871588     25              PINILLOS
## 26 -73.85118  8.722018     26               REGIDOR
## 27 -73.99703  8.493023     27              RIOVIEJO
## 28 -75.07334 10.369721     28         SAN CRISTÓBAL
## 29 -75.20523 10.364682     29        SAN ESTANISLAO
## 30 -74.35870  9.091261     30          SAN FERNANDO
## 31 -75.10609  9.835362     31           SAN JACINTO
## 32 -74.66840  8.237510     32 SAN JACINTO DEL CAUCA
## 33 -75.06572  9.967958     33   SAN JUAN NEPOMUCENO
## 34 -74.01033  8.787692     34    SAN MARTÍN DE LOBA
## 35 -74.13566  7.422366     35             SAN PABLO
## 36 -75.26548 10.649813     36        SANTA CATALINA
## 37 -75.36116 10.469211     37            SANTA ROSA
## 38 -74.26273  7.772519     38    SANTA ROSA DEL SUR
## 39 -73.98372  7.851316     39                SIMITÍ
## 40 -75.11699 10.332090     40           SOPLAVIENTO
## 41 -74.64010  9.313447     41        TALAIGUA NUEVO
## 42 -74.29676  8.492345     42               TIQUISO
## 43 -75.37972 10.353108     43               TURBACO
## 44 -75.47038 10.216047     44               TURBANA
## 45 -75.26736 10.444913     45            VILLANUEVA
## 46 -74.88067  9.745147     46              ZAMBRANO

Ahora, se descargan los datos de elevación para el departamento de Bolivar: A continuación se va a tomar z la cual nos indica el zoom de los datos y en cuanto mayor sea el zoom, mayor será la resolución espacial

elevation <- get_elev_raster(munic, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

Mosaico y Proyección

Indicación de las unidades de elevación que están en metros.

¿Qué es el objeto de elevación? Tomamos el raster generado anteriormente denominado elevation y lo corremos

#elevation 

if (require(rgdal)) {
 rf <- writeRaster(elevation, filename=file.path("C:/Cuaderno3R/bolivar_dem_z10.tif"), format="GTiff", overwrite=TRUE)
}

El tamaño de mi DEM cuando el zoom es 12 es 457.437KB

elevation = raster("C:/Cuaderno3R/bolivar_dem_z10.tif")

Visualización de los datos de elevación

A continuación se va a disminuir la resolución con la finalidad de acelerar la tarea

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

crs(elevation2) <- "+proj=longlat +datum=WGS84"

Es importante a la hora de elebaorar un mapa que, se use una paleta de colores que permita analizar y leer el mapa de una forma adecuada, además de que sea agradable a la hora de visualizar

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

Es momento de usar la libreria leaflet para visualizar los datos de elevación:

leaflet(munic) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", munic$MPIO_CNMBR, "<br>",
                          "Km2: ", munic$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 Bolivar(mts)")

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Reproyección de los datos de elevación

Es de fundamental importancia reconocer que a la hora de usar un DEM se deben manejar coordenadas planas debido a que, en coordenadas geográficas en la dimensión vertical y horizontal se usan dos unidades de medida distintas siendo metros y grados respectivamente

##library(sp)
### WGS 84 / UTM zone 18N ++++++++++++++
spatialref <- CRS(SRS_string="EPSG:32618")

Teniendo en cuenta ahora la nueva configuración que se le dio al DEM, lo vamos a reproyectar con la función projectRaster del paquete raster package

#Se crea el raster incialmente a partir del uso de projectExtent
pr3 <- projectExtent(elevation, crs=spatialref)

En el paso anterior se creó un objeto raster con projectExtent, ahora es necesario que se ajustar el tamaño de la celda

res(pr3) <- 80
rep_elev <- projectRaster(elevation, pr3)

A coninuación se va a leer lo que se configuró para nuestro raster

rep_elev

## class      : RasterLayer 
## dimensions : 5766, 3886, 22406676  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 357510.8, 668390.8, 736684.4, 1197964  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-30_115318_5884_84031.grd 
## names      : bolivar_dem_z10 
## values     : -2790, 5572.246  (min, max)

Ahora reproyectemos el SpatialPolygonDataFrame el cual es un objeto de datos básico para procesos geoestadísticos y se usaprincipalmente para preparar datos, este es aquel que representa a los municipios de Bolivar.

(rep_munic = spTransform(sp.munic,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 46 
## extent      : 369265.3, 638066.3, 773712.6, 1194040  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 9
## names       : DPTO_CCDGO, MPIO_CCDGO, MPIO_CNMBR,                           MPIO_CRSLC,    MPIO_NAREA, MPIO_NANO, DPTO_CNMBR,    Shape_Leng,       Shape_Area 
## min values  :         13,      13001,       ACHÍ,                                 1599,   42.26880612,      2017,    BOLÍVAR, 0.30172488396, 0.00348843338908 
## max values  :         13,      13894,   ZAMBRANO, Ordenanza 42 del 27 de Abril de 1923, 2382.78237094,      2017,    BOLÍVAR, 4.03700051267,   0.195324627752

writeRaster(rep_elev, filename="c:/Cuaderno3R/rep_boivar_dem.tif", datatype='INT4S', overwrite=TRUE)

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

A continuación vamos a ver algunas aplicaciones estadísticas que se le pueden dar a nuestro DEM Inicialmente se corre un código que limpia la elevacion que sea inferior a 0

rep_elev[rep_elev < 0.0] <- 0.0

Teniendo la instrucción anterior, se va a leer un código que permita la creación de un histograma, el cual permite visualizar la distribución de una variable

hist(rep_elev)

Ahora teniendo una herramienta estadística como el histograma, se le va a ordenar a R hacer un resumen de los datos estadísticos del DEM

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

#El siguiente nos permite generar un vector de estadísticas unidimensional
metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)

Teniendo lo anterior ahora se va a representar:

(df_estadisticas <- data.frame(metricas, valores))

##   metricas   valores
## 1     mean  322.3025
## 2      min    0.0000
## 3      max 5572.2461
## 4      std  597.6396

##Obtención de variables geomorfométricas Inicialmente vamos a calcular la pendiente, el aspecto y el sombreado; donde la pendiente es la tasa de cambio de elevación de una celda a la siguiente en un (DEM)

slope = terrain(rep_elev,opt='slope', unit='degrees')
aspect = terrain(rep_elev,opt='aspect',unit='degrees')
hill = hillShade(slope,aspect,40,315)

Ahora se puede ver representandolo así primero:

slope

## class      : RasterLayer 
## dimensions : 5766, 3886, 22406676  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 357510.8, 668390.8, 736684.4, 1197964  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 78.27454  (min, max)

Segunda representación

aspect

## class      : RasterLayer 
## dimensions : 5766, 3886, 22406676  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 357510.8, 668390.8, 736684.4, 1197964  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : aspect 
## values     : 0, 360  (min, max)

A continuación vamos a reducir el tamaño de RasterLayer que se generó a partir de un objeto sp o ráster y nos arrojó nuestro RasterLayer de los mosaicos que se superponen al cuadro delimitador de la entrada

slope2 <- aggregate(slope, fact=4, fun=mean)

Ahora se va a proyectar el conjunto de datos de elevación

slope3 <- projectRasterForLeaflet(slope2, "ngb")

y como lo vamos a representar, le vamos a asignar o crear una paleta de colores que permita visualizar y entender la representación

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

Y a continuación si vamos a ordenar la representación:

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

## Warning in colors(.): Some values were outside the color scale and will be
## treated as NA

Ahora se traerá de nuevo el objeto slope2 creado anteriormente:

slope2

## class      : RasterLayer 
## dimensions : 1442, 972, 1401624  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 357510.8, 668550.8, 736524.4, 1197964  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-30_115432_5884_11031.grd 
## names      : slope 
## values     : 0, 70.83002  (min, max)

Como vamos a manejar un raster se usará la librería terra pero, ¿por qué no usar la libreria raster? porque es mas lenta que la terra aunque esta última sea mas reciente, es mas reocmendada usarla, teniendo en cuenta lo anterior convertiremos un RasterLayer en un SpatRaster que es un objeto terra) así:

(nslope2  <- as(slope2, "SpatRaster"))

## class       : SpatRaster 
## dimensions  : 1442, 972, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 357510.8, 668550.8, 736524.4, 1197964  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source      : r_tmp_2023-05-30_115432_5884_11031.grd 
## name        :    slope 
## min value   :  0.00000 
## max value   : 70.83002

A continuación vamos a hacer una clasificación de pendientes la cual permite entender la variabilidad del relieve, esta debe contar con intervalos de pendiente, nombres de clase y descripción. Pero antes crearemos una matriz de reclasificación.

m <- c(0, 3, 1,  3, 7, 2, 7,12,3,  12,
       25,  4, 25, 50, 5, 50, 75, 6)
(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

y posteriormente se realiza la reclasificación:

slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)

slope_rec

## class       : SpatRaster 
## dimensions  : 1442, 972, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 357510.8, 668550.8, 736524.4, 1197964  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        : slope 
## min value   :     0 
## max value   :     6

##Visualización Para esta sección se va a usar la libreria tmap para crear y guardar un mapa estático que indica la representación de la distribución en el territorio, en el cual vamos a ver la clasificación de pendientes.

slope.map <- tm_shape(slope_rec) +
  tm_raster(alpha = 0.6, style= "cat", 
            title="Slope") +  
  # Borders only
  tm_shape(munic) + 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)  

Y ahora si la representación de la pendiente:

slope.map

## stars object downsampled to 821 by 1218 cells. See tm_shape manual (argument raster.downsample)

Ahora, vamos a generar una copia de seguridad del mapa en el directorio donde se tiene guardada toda la información para el desarrollo del cuaderno y se guarda con determinadas medidas porque se va a usar una hoja tamaño carta.

tmap_save(slope.map, "C:/Cuaderno3R/bolivar_slope.png", width = 8.5, height =11, dpi=600)

## stars object downsampled to 821 by 1218 cells. See tm_shape manual (argument raster.downsample)

## Map saved to C:\Cuaderno3R\bolivar_slope.png

## Resolution: 5100 by 6600 pixels

## Size: 8.5 by 11 inches (600 dpi)

##Visualización interactiva Para esta sección se va a generar un mapa de pendientes interactivo y para esto se convierte el objeto sf incial en un objeto SpatVector:

(nmunic  <- as(rep_munic, "SpatVector"))

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 46, 9  (geometries, attributes)
##  extent      : 369265.3, 638066.3, 773712.6, 1194040  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       : DPTO_CCDGO MPIO_CCDGO        MPIO_CNMBR      MPIO_CRSLC
##  type        :      <chr>      <chr>             <chr>           <chr>
##  values      :         13      13001   CARTAGENA DE I~            1599
##                        13      13006              ACHÍ Ordenanza 15 d~
##                        13      13030 ALTOS DEL ROSARIO Ordenanza 30 d~
##  MPIO_NAREA MPIO_NANO DPTO_CNMBR Shape_Leng Shape_Area
##       <num>     <int>      <chr>      <num>      <num>
##       593.2      2017    BOLÍVAR      4.037    0.04912
##       951.8      2017    BOLÍVAR      1.837    0.07818
##       303.7      2017    BOLÍVAR      0.985    0.02496

Ahora, como se hizo anteriormente vamos a convertir el objeto raster en un objeto terra:

(nslope  <- as(slope, "SpatRaster"))

## class       : SpatRaster 
## dimensions  : 5766, 3886, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 357510.8, 668390.8, 736684.4, 1197964  (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   : 78.27454

Ahora, vamos a calcular las categorías de pendientes agregadas por municipio del departamento de Bolivar:

munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)

Nuestro objeto terra será el siguiente:

munic.slope

## class       : SpatRaster 
## dimensions  : 5766, 3886, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 357510.8, 668390.8, 736684.4, 1197964  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :      slope 
## min value   :  0.5844179 
## max value   : 14.0517475

A continuación el mapa interactivo de las pendientes:

tmap_mode("view")

## tmap mode set to interactive viewing

  tm_shape(munic) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "MPIO_CNMBR", size = 0.5, alpha=0.7) +
  tm_shape(munic.slope) +
  tm_raster(style="cont", n=7, alpha=0.6, title="Mean slope")+
tm_layout(legend.outside = TRUE)

## stars object downsampled to 821 by 1218 cells. See tm_shape manual (argument raster.downsample)

##Bibliografia Lizarazo, I., 2023. Procesamiento y análisis de datos de elevación en R. Disponible en https://rpubs.com/ials2un/dem_analysis

sessionInfo()

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