Limpiar la memoria de R.

rm(list=ls())

Cargar las librerias previamente instaladas en la consola.

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)

2. Introducción a elevatr

LOs datos de elevación son usados en gran cantidad de aplicaciones, incluyendo visualizacion, hidrologia y el modelo ecologico. El paquete “elevatr” funciona para estandarizar el acceso a los datos de elevación desde la web APIs. Durante el desarrollo del siguiente cuaderno, se obtendran tanto mapas estáticos como interactivos, los cuales reflejarán los municipios del departamento del meta con sus respectivas elevaciones y tipos de pendiente.

3. Obtener datos de elvación raster del departamento

Primero, se carga un shapefile o un geopackage que represente los municipios del departamento que se esta trabajando, en este caso el Meta. En este cuaderno, se usará un geopackage, el cual fue descargado como shapefile originalmente. Sin embargo, este archivo abarcaba todos los municipios de Colombia, por lo que se recortó por medio del programa Qgis para obtener solo los municipios del departamento del Meta. Mediante este archivo, se crea un objeto, el cual sera llamado “munic”.

munic <-  st_read("C:/GB2/CUADERNO_3/CUADERNO_3/MUNI_META.gpkg")
## Reading layer `MUNI_META' from data source 
##   `C:\GB2\CUADERNO_3\CUADERNO_3\MUNI_META.gpkg' using driver `GPKG'
## Simple feature collection with 29 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.89921 ymin: 1.604238 xmax: -71.07753 ymax: 4.899101
## Geodetic CRS:  MAGNA-SIRGAS

Se revisa la geometria, las coordenadas, CRS y los atributos del objeto “munic”.

head(munic)
## Simple feature collection with 6 features and 10 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.39891 ymin: 3.704314 xmax: -72.73891 ymax: 4.736325
## Geodetic CRS:  MAGNA-SIRGAS
##   DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR MPIO_CDPMP VERSION       AREA
## 1         50        001          VILLAVICENCIO      50001    2018 1285879354
## 2         50        006          ACACÃ\u008dAS      50006    2018 1123737929
## 3         50        110 BARRANCA DE UPÃ\u008dA      50110    2018  403818863
## 4         50        124               CABUYARO      50124    2018  913662509
## 5         50        150      CASTILLA LA NUEVA      50150    2018  514685967
## 6         50        223               CUBARRAL      50223    2018 1157865246
##    LATITUD  LONGITUD Shape_Leng Shape_Area                           geom
## 1 4.091669 -73.49292   2.039775 0.10471406 MULTIPOLYGON (((-73.69288 4...
## 2 4.009644 -73.72391   2.074207 0.09150745 MULTIPOLYGON (((-73.80268 4...
## 3 4.519076 -72.99549   1.312018 0.03289441 MULTIPOLYGON (((-73.01607 4...
## 4 4.315245 -72.95269   2.008848 0.07440326 MULTIPOLYGON (((-73.08865 4...
## 5 3.805154 -73.53887   1.349329 0.04189962 MULTIPOLYGON (((-73.7252 3....
## 6 3.828242 -74.07530   2.075091 0.09426999 MULTIPOLYGON (((-74.20143 4...

El siguiente paso es crear un nuevo objeto con los centroides de los municipios, el cual se llamará “centers”. El objetivo es encontrar un punto central para cada región.

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

Tabla de atributos del nuevo objeto creado “centers”.

centers
##            x        y region                  munic
## 1  -73.49292 4.091669      1          VILLAVICENCIO
## 2  -73.72391 4.009644      2          ACACÃ\u008dAS
## 3  -72.99549 4.519076      3 BARRANCA DE UPÃ\u008dA
## 4  -72.95269 4.315245      4               CABUYARO
## 5  -73.53887 3.805154      5      CASTILLA LA NUEVA
## 6  -74.07530 3.828242      6               CUBARRAL
## 7  -73.31478 4.232526      7                CUMARAL
## 8  -73.71442 4.353770      8            EL CALVARIO
## 9  -73.88801 3.615625      9            EL CASTILLO
## 10 -73.83162 3.706970     10              EL DORADO
## 11 -73.59625 3.382370     11          FUENTE DE ORO
## 12 -73.74696 3.496563     12                GRANADA
## 13 -73.95984 3.947776     13                 GUAMAL
## 14 -74.12455 3.154771     14                MESETAS
## 15 -74.43066 3.047727     15                  URIBE
## 16 -74.09628 3.614715     16         LEJANÃ\u008dAS
## 17 -72.72109 2.752260     17       PUERTO CONCORDIA
## 18 -72.64570 4.014299     18          PUERTO LÓPEZ
## 19 -73.23671 3.193093     19          PUERTO LLERAS
## 20 -73.13780 2.758084     20            PUERTO RICO
## 21 -73.47367 4.237467     21               RESTREPO
## 22 -73.27583 3.847617     22   SAN CARLOS DE GUAROA
## 23 -73.81635 3.289851     23      SAN JUAN DE ARAMA
## 24 -73.66041 4.471346     24            SAN JUANITO
## 25 -72.98537 3.535909     25       SAN MARTÃ\u008dN
## 26 -73.66652 2.811617     26           VISTAHERMOSA
## 27 -74.09488 2.161864     27            LA MACARENA
## 28 -71.93806 3.117523     28        MAPIRIPÃ\u0081N
## 29 -71.63157 4.005034     29    PUERTO GAITÃ\u0081N
# z denotes the zoom level of the data
# the higher the zoom the higher the spatial resolution
#elevation <- get_elev_raster(munic, z = 12)
elevation <- get_elev_raster(munic, z = 10)
## Mosaicing & Projecting
## Note: Elevation units are in meters.

elevation
## class      : RasterLayer 
## dimensions : 5115, 6148, 31447020  (nrow, ncol, ncell)
## resolution : 0.0006861734, 0.0006861734  (x, y)
## extent     : -75.23438, -71.01578, 1.406056, 4.915833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +no_defs 
## source     : file2fe469c25d87.tif 
## names      : file2fe469c25d87
# uncomment if needed 
# write to geotiff file (depends on rgdal)
# NOTE THAT z10 means Zoom Level 10
# CHANGE IT IF NEEDED
#if (require(rgdal)) {
rf <- writeRaster(elevation, filename=file.path("C:/GB2/CUADERNO_3/CUADERNO_3/META_dem_z10.tif"), format="GTiff", overwrite=TRUE)
#}

se guarda el archivo en el directorio de trabajo.

(elevation = raster("C:/GB2/CUADERNO_3/CUADERNO_3/META_dem_z10.tif"))
## class      : RasterLayer 
## dimensions : 5115, 6148, 31447020  (nrow, ncol, ncell)
## resolution : 0.0006861734, 0.0006861734  (x, y)
## extent     : -75.23438, -71.01578, 1.406056, 4.915833  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +ellps=GRS80 +no_defs 
## source     : META_dem_z10.tif 
## names      : META_dem_z10 
## values     : -838, 4359  (min, max)

4. Visualizar los datos de elevación

El siguente paso es visualizar los datos de elevacion y se van a convertir en un nuevo objeto, el cual se llamará “elevation2”.

#Use it only when zoom data was very high
elevation2 <- aggregate(elevation, fact=4, fun=mean)
crs(elevation2) <- "+proj=longlat +datum=WGS84"

Para una buena visualización, se debe tener una buena paleta de colores. En este caso se usara la siguiente:

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

Ahora, con la libreria “leaflet” se pueden ver los datos de elvació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 Elevacion para el Meta (mts)")

5. Reproyectar los datos de elevación

Cuando se trabaja con DEM, es mejor usar las coordenadas de mapa, ya que las coordenadas geográficas utilizan decimales en las dimensiones horizontales, y en las dimensiones verticales utilizan metros, por lo que esto resulta no ser muy conveniente para trabajar durante el desarrollo de este cuaderno.

Para reproyectar los datos de elevacion se tienen 2 pasos:

Primero, se define el objetico CRS:

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

Luego, se reproyecta el DEM usando la función projectRaster del paquete Raster:

# Proyectar raster
# Primero, se crea un objeto raster
# using projectExtent (no values are transferred)
pr3 <- projectExtent(elevation, crs=spatialref)
# Adjust the cell size to your cell size
res(pr3) <- 80
# now project
rep_elev <- projectRaster(elevation, pr3)

Que se obtiene?, lo siguiente:

rep_elev
## class      : RasterLayer 
## dimensions : 4866, 5871, 28568286  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 473927.7, 943607.7, 155399.3, 544679.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : r_tmp_2023-05-31_124906_12260_81421.grd 
## names      : layer 
## values     : -576.3135, 4248.04  (min, max)

Ahora, se van a reproyectar los polígonos espaciales que representan los municipios del departamento trabajado (Meta).

(rep_munic = spTransform(sp.munic,spatialref))
## class       : SpatialPolygonsDataFrame 
## features    : 29 
## extent      : 511200.9, 936266.2, 177345.3, 542784.7  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 10
## names       : DPTO_CCDGO, MPIO_CCDGO,    MPIO_CNMBR, MPIO_CDPMP, VERSION,          AREA,       LATITUD,       LONGITUD,    Shape_Leng,    Shape_Area 
## min values  :         50,        001, ACACÍAS,      50001,    2018, 117341083.432, 2.16186391692, -74.4306613394, 0.43640580003, 0.00955216819 
## max values  :         50,        711,  VISTAHERMOSA,      50711,    2018, 17250083166.3, 4.51907612664, -71.6315742854, 8.63318071192, 1.40216640362

Para evitar cualquier inconveniente. Se procede a guardar el DEM. Como recomendación, asegurarse de indicar la dirección donde desea guardarlo, es decir, el directorio de trabajo donde se encuentran todos los datos que se han utilizado hasta el momento para este cuaderno y el nombre del archivo.

writeRaster(rep_elev, filename="C:/GB2/CUADERNO_3/CUADERNO_3/rep_META_dem.tif", datatype='INT4S', overwrite=TRUE)

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

Una exploración de las estadísticas del DEM que se esta trabajando.

Este código “elimina” la elevacion menor a 0.0:

rep_elev[rep_elev < 0.0] <- 0.0

Ahora, se crea un histograma con los valores de elevacion:

# histograma
hist(rep_elev)

Se pueden calcular las estadísticas básicas del DEM:

# works for large files
promedio <- cellStats(rep_elev, 'mean')
minimo <- cellStats(rep_elev, 'min')
maximo <- cellStats(rep_elev, 'max')
desviacion  <- cellStats(rep_elev, 'sd')

El siguiente paso es crear un vertor de estadísticas unidimensional para poder realizar los cálculos.

metricas <- c('mean', 'min', 'max', 'std')
valores <- c(promedio, minimo, maximo, desviacion)

Como resultado se tiene que:

(df_estadisticas <- data.frame(metricas, valores))
##   metricas   valores
## 1     mean  572.4375
## 2      min    0.0000
## 3      max 4248.0396
## 4      std  769.4605

7. Obtención de variables geomorfométricas

Antes de inciar con este paso, es importante tener claros los conceptos de DEM, slope and aspect.

Primero, se calcula slope (pendiente), aspect (orientacion) y finalmente hillshade (sombreado).

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

Vease la primera variable (slope).

slope
## class      : RasterLayer 
## dimensions : 4866, 5871, 28568286  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 473927.7, 943607.7, 155399.3, 544679.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 87.38766  (min, max)

Ahora, la segunda (aspect).

aspect
## class      : RasterLayer 
## dimensions : 4866, 5871, 28568286  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 473927.7, 943607.7, 155399.3, 544679.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : aspect 
## values     : 0, 360  (min, max)

Para trazar la pendiente (slope) es conveniente reducir el tamaño de la capa Raster:

#This chunk is optional
#Use it only when zoom data was very high
slope2 <- aggregate(slope, fact=4, fun=mean)

Se puede proyectar el conjunto de datos de elevación:

slope3 <- projectRasterForLeaflet(slope2, "ngb")

Luego, se crea la paleta de colores.

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

Ahora, el momento de trazar ha llegado.

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 meta (%)")

Recuerde que es el objeto “slope2”:

slope2
## class      : RasterLayer 
## dimensions : 1217, 1468, 1786556  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 473927.7, 943687.7, 155239.3, 544679.3  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0.03053321, 85.52554  (min, max)

Para manejar datos tipo raster en R, se pueden utilizar las librerias (raster y terra)

Se va a convertir un objeto raster a un objeto terra. Este nuevo objeto se llamará “nslope2”.

(nslope2  <- as(slope2, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 1217, 1468, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 473927.7, 943687.7, 155239.3, 544679.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :       slope 
## min value   :  0.03053321 
## max value   : 85.52554275

Mediante la clasificación de la pendiente se puede entender la variablididad del relieve. Se necesita incluir una tabla con intervalos, clases y descripción según lo indica el instituto geográfico de Colombia (IGAC)

Se pude crear una clasificacion de pendientes con rangos comunes. Para esto, primero se crea 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

La matriz anterior esta compuesta por rangos comunes de pendiente, los cuales permitirán identificar el valor de pendiente y si se profundiza más también la clase y la descripción de la misma.

Como queda la reclasificación?,Algo así:

# for values >= 0 (instead of > 0), do
slope_rec <- terra::classify(nslope2, rclmat, right=TRUE)

Observe el resultado:

slope_rec
## class       : SpatRaster 
## dimensions  : 1217, 1468, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 473927.7, 943687.7, 155239.3, 544679.3  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :    slope 
## min value   :  1.00000 
## max value   : 85.52554

8. Visualización estática.

Para este paso, se usará la libreria “tmap” para crear y guardar un mapa estático de clasificación de pendiente.

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)

EL siguiente paso es ver como quedó el mapa.

slope.map

El departamento del Meta es principalmente conocido por sus llanuras, por lo que sus clases de pendientes sera bajas en casi toda su distribución geógrafica.

Se procede a guardar el mapa anterior.

tmap_save(slope.map, "C:/GB2/CUADERNO_3/meta_slope.png", width = 8.5, height =11, dpi=600)
## stars object downsampled to 1098 by 910 cells. See tm_shape manual (argument raster.downsample)
## Map saved to C:\GB2\CUADERNO_3\meta_slope.png
## Resolution: 5100 by 6600 pixels
## Size: 8.5 by 11 inches (600 dpi)

9. Visualización Interactiva.

Tambien se puede crear un mapa interactivo de la siguiente manera:

Se convierte un objeto sf en un objeto “vector espacial”:

(nmunic  <- as(rep_munic, "SpatVector"))
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 29, 10  (geometries, attributes)
##  extent      : 511200.9, 936266.2, 177345.3, 542784.7  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       : DPTO_CCDGO MPIO_CCDGO             MPIO_CNMBR MPIO_CDPMP VERSION
##  type        :      <chr>      <chr>                  <chr>      <chr>   <int>
##  values      :         50        001          VILLAVICENCIO      50001    2018
##                        50        006          ACACÃ\u008dAS      50006    2018
##                        50        110 BARRANCA DE UPÃ\u008dA      50110    2018
##       AREA LATITUD LONGITUD Shape_Leng Shape_Area
##      <num>   <num>    <num>      <num>      <num>
##  1.286e+09   4.092   -73.49       2.04     0.1047
##  1.124e+09    4.01   -73.72      2.074    0.09151
##  4.038e+08   4.519      -73      1.312    0.03289

se convierte el objeto raster a objeto terra con el siguiente código:

(nslope  <- as(slope, "SpatRaster"))
## class       : SpatRaster 
## dimensions  : 4866, 5871, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 473927.7, 943607.7, 155399.3, 544679.3  (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   : 87.38766

Ahora se calculan las categorias de pendientes agregadas por municipio. Este onjeto se llamará “munic.slope”

## S4 method for signature 'SpatRaster,SpatVector'
munic.slope <- zonal(nslope, nmunic, fun=mean, as.raster=TRUE)

El resultado del objeto terra es:

munic.slope
## class       : SpatRaster 
## dimensions  : 4866, 5871, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 473927.7, 943607.7, 155399.3, 544679.3  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :      slope 
## min value   :  0.6657628 
## max value   : 26.3263620

Finalmente este es el resultado del mapa interactivo:

tmap_mode("view")
  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="Pendiente media")+
tm_layout(legend.outside = TRUE)

10. Bibliografía.

Para la elaboración de este cuaderno se uso como guia el siguiente cuaderno:

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

sessionInfo()
## R version 4.2.1 Patched (2022-09-21 r82893 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
## 
## 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.2       
##  [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-3        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.1             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.1          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.1              parallel_4.2.1          crayon_1.5.2           
## [58] deldir_1.0-9            stars_0.6-1             hms_1.1.2              
## [61] knitr_1.42              pillar_1.9.0            markdown_1.5           
## [64] codetools_0.2-18        stats4_4.2.1            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.1           
## [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-20            viridisLite_0.4.1       units_0.8-2            
## [85] timechange_0.2.0        ellipsis_0.3.2