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

Introducción

En este cuaderno se ilustra como obtener, procesar y visualizar modelos digitales de elevación (DEM) en R.

Configuración

Instalar los paquetes necesarios para elaborar el trabajo, estos son: “rgdal”,“raster”, “rgeos”, terra”,“elevatr”,“rasterVis”, “rgl”, “mapview”

Limpiar memoria

rm(list=ls())

Cargar librerias

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)

elevatr

Se usa el paquete elevatr para acceder a los datos de elevación desde las APl web

Obtención de datos raster para el departamento de interes

Primero se cargara un shapefile que represente los municipios de nuestro(s) departamento(S). En este caso se usara un shapefile que contiene los municipios de Caldas, Quindio y Risaralda.

(mun.eje =  st_read('C:/Users/LENOVO/Desktop/Geomatica/R/Datos/muni_eje.shp'))

## Reading layer `muni_eje' from data source 
##   `C:\Users\LENOVO\Desktop\Geomatica\R\Datos\muni_eje.shp' using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 14 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
## Geodetic CRS:  WGS 84

## Simple feature collection with 51 features and 14 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -76.3595 ymin: 4.135801 xmax: -74.6489 ymax: 5.7626
## Geodetic CRS:  WGS 84
## First 10 features:
##    fid ID_0 ISO   NAME_0 ID_1 NAME_1 ID_2     NAME_2    TYPE_2    ENGTYPE_2
## 1    1   53 COL Colombia    8 Caldas  348    Aguadas Municipio Municipality
## 2    2   53 COL Colombia    8 Caldas  349    Anserma Municipio Municipality
## 3    3   53 COL Colombia    8 Caldas  350   Aranzazú Municipio Municipality
## 4    4   53 COL Colombia    8 Caldas  351 Belalcázar Municipio Municipality
## 5    5   53 COL Colombia    8 Caldas  352  Chinchiná Municipio Municipality
## 6    6   53 COL Colombia    8 Caldas  353 Filadelfia Municipio Municipality
## 7    7   53 COL Colombia    8 Caldas  354  La Dorada Municipio Municipality
## 8    8   53 COL Colombia    8 Caldas  355  La Merced Municipio Municipality
## 9    9   53 COL Colombia    8 Caldas  356  Manizales Municipio Municipality
## 10  10   53 COL Colombia    8 Caldas  357 Manzanares Municipio Municipality
##    NL_NAME_2 VARNAME_2      Area codigo                       geometry
## 1       <NA>      <NA> 501.64063  17013 POLYGON ((-75.3399 5.458, -...
## 2       <NA>      <NA> 196.73857  17616 POLYGON ((-75.6908 5.152099...
## 3       <NA>  Aranzazé 140.63766  17050 POLYGON ((-75.3997 5.2202, ...
## 4       <NA>      <NA> 114.27085  17088 POLYGON ((-75.7841 5.032401...
## 5       <NA>      <NA> 125.39270  17174 POLYGON ((-75.6204 4.979101...
## 6       <NA>      <NA> 202.22182  17272 POLYGON ((-75.5888 5.207699...
## 7       <NA>      <NA> 583.68634  17380 POLYGON ((-74.8504 5.290101...
## 8       <NA>      <NA>  90.16428  17388 POLYGON ((-75.6437 5.362199...
## 9       <NA>      <NA> 443.87241  17001 POLYGON ((-75.636 5.0011, -...
## 10      <NA>      <NA> 171.87523  17433 POLYGON ((-75.1662 5.1453, ...

municipioseje <- data.frame(mun.eje)

Despues se comprueba la geometria, el cuadro delimitador, el CRS y los atributos del objeto municipioseje

head(municipioseje)

##   fid ID_0 ISO   NAME_0 ID_1 NAME_1 ID_2     NAME_2    TYPE_2    ENGTYPE_2
## 1   1   53 COL Colombia    8 Caldas  348    Aguadas Municipio Municipality
## 2   2   53 COL Colombia    8 Caldas  349    Anserma Municipio Municipality
## 3   3   53 COL Colombia    8 Caldas  350   Aranzazú Municipio Municipality
## 4   4   53 COL Colombia    8 Caldas  351 Belalcázar Municipio Municipality
## 5   5   53 COL Colombia    8 Caldas  352  Chinchiná Municipio Municipality
## 6   6   53 COL Colombia    8 Caldas  353 Filadelfia Municipio Municipality
##   NL_NAME_2 VARNAME_2     Area codigo                       geometry
## 1      <NA>      <NA> 501.6406  17013 POLYGON ((-75.3399 5.458, -...
## 2      <NA>      <NA> 196.7386  17616 POLYGON ((-75.6908 5.152099...
## 3      <NA>  Aranzazé 140.6377  17050 POLYGON ((-75.3997 5.2202, ...
## 4      <NA>      <NA> 114.2709  17088 POLYGON ((-75.7841 5.032401...
## 5      <NA>      <NA> 125.3927  17174 POLYGON ((-75.6204 4.979101...
## 6      <NA>      <NA> 202.2218  17272 POLYGON ((-75.5888 5.207699...

Ahora, crearemos un nuevo objeto con los centroides de los municipios

# La tarea consiste en encontar un punto central para cada región
#Convertimos las caracteristicas simples en poligonos espaciales
sp.mun.eje = as_Spatial(mun.eje)

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

head(centers)

##           x        y region      munic
## 1 -75.48755 5.590398      1    Aguadas
## 2 -75.79169 5.204271      2    Anserma
## 3 -75.50984 5.255189      3   Aranzazú
## 4 -75.85170 4.981161      4 Belalcázar
## 5 -75.70418 4.975228      5  Chinchiná
## 6 -75.62369 5.282661      6 Filadelfia

Descargaremos los datos de elevación para nuestro(s) departamento(s)

# z indica el nivel de zoom de los datos
# a mayor zoom mayor resolución espacial

elevation <- get_elev_raster(mun.eje, z = 10)

## Mosaicing & Projecting

## Note: Elevation units are in meters.

elevation

## class      : RasterLayer 
## dimensions : 3066, 3078, 9437148  (nrow, ncol, ncell)
## resolution : 0.0006853648, 0.0006853648  (x, y)
## extent     : -76.64062, -74.53107, 3.864425, 5.965754  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : file175c72b8762c.tif 
## names      : file175c72b8762c

Se tiene que el tamaño de celda es: 0.0006853648, 0.0006853648 (x, y)

La elevacion es una capa raster temporal, sin embargo, es posible escribir el DEM en su directorio local usando la función writeRaster.

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

Visualizar los datos de elevación

Para visualizar, una resolución mas baja agiliza la tarea.

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

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

Una buena paleta es clave para una visulización adecuada. Probaremos con la siguiente:

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

Ahora, utilizaremos la biblioteca leaflet para ver los datos de elevación.

leaflet(mun.eje) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", mun.eje$NAME_2, "<br>",
                          "Km2: ", mun.eje$Area, "<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 eje cafetero (mts)")

Haz clic en el ID de cada región para obtener el nombre y la extensión del municipio.

Reproyeccón de los datos de elevación

Cuando se trabaja con un DEM, se recomienda utilizar coordenadas cartograficas en lugar de coordenadas geograficas.

Para reproyectar los datos de elevación en coordenadas cartograficas realizaremos dos pasos.

Primero. Definir el CRS del objeto.

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

Despues, reproyectar el DEM usando la función projectRaster

# Primero creamos un objeto raster
pr3 <- projectExtent(elevation, crs=spatialref)

# Ajusta el tamaño de celda al tuyo
res(pr3) <- 80

#proyectalo
rep_elev <- projectRaster(elevation, pr3)

¿Que se obtiene?

rep_elev

## class      : RasterLayer 
## dimensions : 2907, 2928, 8511696  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317827, 552067, 427129.2, 659689.2  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : file175c72b8762c 
## values     : -267.5508, 5280.105  (min, max)

Ahora, reproyectaremos el SpatialPolygons DataFrame, el cual representa los municipios de nuestro(s) departamento(s)

(rep_mun.eje = spTransform(sp.mun.eje,spatialref))

## class       : SpatialPolygonsDataFrame 
## features    : 51 
## extent      : 349337, 538878.8, 457177.6, 636970.4  (xmin, xmax, ymin, ymax)
## crs         : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## variables   : 14
## names       : fid, ID_0, ISO,   NAME_0, ID_1,    NAME_1, ID_2,  NAME_2,    TYPE_2,    ENGTYPE_2, NL_NAME_2,   VARNAME_2,             Area, codigo 
## min values  :   1,   53, COL, Colombia,    8,    Caldas,  348, Aguadas, Municipio, Municipality,        NA,    Aranzazé, 30.5922403667703,  17001 
## max values  :  51,   53, COL, Colombia,   25, Risaralda,  851, Viterbo, Municipio, Municipality,        NA, Mantequilla,  1033.6076464384,  66687

Para evitar problemas, guardaremos nuestro DEM.

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

Estadisticas basicas de los datos de elevación

Debemos limpiar la elevación inferior a 0.0

rep_elev[rep_elev < 0.0] <- 0.0

Vamos a crear histograma con los valores de elevación

hist(rep_elev, main = "Histograma elevacion", col= "green")

Resumen de las estadisticas DEM

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

Ahora crearemos un vector unidimensional con las estadisticas

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

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

##   metricas   valores
## 1     mean 1358.0057
## 2      min    0.0000
## 3      max 5280.1046
## 4      std  962.5001

Obtención de variables geomorfometricas

Calcularemos la pendiente, el aspecto y el relieve

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

Resultados:

Pendiente

slope

## class      : RasterLayer 
## dimensions : 2907, 2928, 8511696  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317827, 552067, 427129.2, 659689.2  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0, 85.36249  (min, max)

Aspecto

aspect

## class      : RasterLayer 
## dimensions : 2907, 2928, 8511696  (nrow, ncol, ncell)
## resolution : 80, 80  (x, y)
## extent     : 317827, 552067, 427129.2, 659689.2  (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 hacer un plot de la pendiente debemos reducir el tamaño de la capa raster

# Usalo cuando el zoom sea muy alto
slope2 <- aggregate(slope, fact=4, fun=mean)

Es conveniente proyectar los datos de elevación

slope3 <- projectRasterForLeaflet(slope2, "ngb")

Tambien crear una paleta de colores

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

Ahora, realizamos el plot

leaflet(mun.eje) %>% addTiles() %>%
    addPolygons(color = "white", weight = 1.0, smoothFactor = 0.5,
    opacity = 0.25, fillOpacity = 0.15,
    popup = paste("Region: ", mun.eje$NAME_2, "<br>",
                          "Km2: ", mun.eje$Area, "<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 = "Pendiente Eje cafetero (%)")

¿Que es el objeto slope2?

slope2

## class      : RasterLayer 
## dimensions : 727, 732, 532164  (nrow, ncol, ncell)
## resolution : 320, 320  (x, y)
## extent     : 317827, 552067, 427049.2, 659689.2  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : slope 
## values     : 0.09009511, 83.27359  (min, max)

Ahora, Convertiremos un Rasterlayer (objeto raster) en un SpatRaster (objeto terra)

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

## class       : SpatRaster 
## dimensions  : 727, 732, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 317827, 552067, 427049.2, 659689.2  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :       slope 
## min value   :  0.09009511 
## max value   : 83.27359377

La clasificación de pendientes es importante para entender la variabilidad del relieve. Es necesario introducir una tabla con los intervalos de pendiente, los nombres de las clases y la descripción, tal como establece el IGAC.

Se realizará una clasificación de pendientes basada en intervalos comunes.

Primero, creamos 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, 75, 100, 7)
(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
## [7,]   75  100    7

La matriz anterior clasifica en diferentes categorias la pendiente que se presenta en un terreno a partir de su valor expresado en porcentaje. Esto muestra que hay 7 clasificaciones categoricas para la pendiente, tal como se puede observar en la ultima columna. La interpretacion para las categorias de pendiente 1, 2, 3, 4, 5, 6, 7, obtenidas en la ultima columna, son respectivamente: A nivel, Ligeramente inclinada, Moderadamente inclinada, Fuertemente inclinada, Ligeramente empinada, Moderadamente empinada, Fuertemente empinada.

Reclasificación:

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

slope_rec

## class       : SpatRaster 
## dimensions  : 727, 732, 1  (nrow, ncol, nlyr)
## resolution  : 320, 320  (x, y)
## extent      : 317827, 552067, 427049.2, 659689.2  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        : slope 
## min value   :     1 
## max value   :     7

Visualización estática

Usaremos tmap para crear y guardar un mapa estatico de clasificación de pendientes.

slope.map <- tm_shape(slope_rec) +
  tm_raster(alpha = 0.6, style= "cat", 
            title="Pendiente") +  
  # Borders only
  tm_shape(mun.eje) + tm_borders(col = "darkgray", lwd = 0.5, alpha=0.8) +
  tm_text(text = "NAME_2", 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= 'clases de pendiente', title.position = c('right', 'top')) +
  tm_credits("Data source: Mapzen & DANE", size=0.3)

Visualizar

slope.map

Guardaremos el mapa

tmap_save(slope.map, "./eje_slope.png", width = 8.5, height =11, dpi=600)

## Map saved to C:\Users\LENOVO\Desktop\Geomatica\R\eje_slope.png

## Resolution: 5100 by 6600 pixels

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

Visualización interactiva

Ahora haremos un mapa de pendiente interactivo.

Convertir el objeto sf en un objeto SpatVector

(nmun.eje  <- as(rep_mun.eje, "SpatVector"))

##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 51, 14  (geometries, attributes)
##  extent      : 349337, 538878.8, 457177.6, 636970.4  (xmin, xmax, ymin, ymax)
##  coord. ref. : +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs 
##  names       :   fid  ID_0   ISO   NAME_0  ID_1 NAME_1  ID_2   NAME_2    TYPE_2
##  type        : <num> <num> <chr>    <chr> <num>  <chr> <num>    <chr>     <chr>
##  values      :     1    53   COL Colombia     8 Caldas   348  Aguadas Municipio
##                    2    53   COL Colombia     8 Caldas   349  Anserma Municipio
##                    3    53   COL Colombia     8 Caldas   350 Aranzazú Municipio
##     ENGTYPE_2 (and 4 more)
##         <chr>             
##  Municipality             
##  Municipality             
##  Municipality

Convertir el objeto raster en un objeto terra

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

## class       : SpatRaster 
## dimensions  : 2907, 2928, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 317827, 552067, 427129.2, 659689.2  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :    slope 
## min value   :  0.00000 
## max value   : 85.36249

Calculemos las categorias de pendientes agregadas por municipio:

mun.eje.slope <- zonal(nslope, nmun.eje, fun=mean, as.raster=TRUE)

Objeto terra:

mun.eje.slope

## class       : SpatRaster 
## dimensions  : 2907, 2928, 1  (nrow, ncol, nlyr)
## resolution  : 80, 80  (x, y)
## extent      : 317827, 552067, 427129.2, 659689.2  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 18N (EPSG:32618) 
## source(s)   : memory
## name        :     slope 
## min value   :  2.281582 
## max value   : 28.282159

Visualización del mapa

tmap_mode("view")

## tmap mode set to interactive viewing

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

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

Bibliografia

Este trabajo se realizó tomando como guía el cuaderno: Lizarazo, I. 2023, Elevation data processing and analysis in R.

Cite este trabajo así: Cely, N., 2023. Procesamiento y análisis de datos de elevación en R. Disponible en https://rpubs.com/ncely/1048046

Reproducibilidad

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=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.1      
## [13] ggplot2_3.4.2      tidyverse_2.0.0    sf_1.0-13          elevatr_0.4.2     
## [17] rgdal_1.6-6        rgeos_0.6-3        raster_3.6-20      sp_1.6-0          
## [21] rasterVis_0.51.5   lattice_0.20-45   
## 
## loaded via a namespace (and not attached):
##  [1] satellite_1.0.4         httr_1.4.6              progress_1.2.2         
##  [4] webshot_0.5.4           tools_4.2.2             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.2          progressr_0.13.0        leafem_0.2.0           
## [19] cli_3.6.1               sass_0.4.6              scales_1.2.1           
## [22] classInt_0.4-9          hexbin_1.28.3           proxy_0.4-27           
## [25] digest_0.6.31           rmarkdown_2.21          base64enc_0.1-3        
## [28] dichromat_2.0-0.1       jpeg_0.1-10             pkgconfig_2.0.3        
## [31] htmltools_0.5.5         highr_0.10              fastmap_1.1.1          
## [34] htmlwidgets_1.6.2       rlang_1.1.1             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.4                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         terra_1.7-29           
## [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.3               knitr_1.43              pillar_1.9.0           
## [64] markdown_1.7            codetools_0.2-18        stats4_4.2.2           
## [67] wk_0.7.3                XML_3.99-0.14           glue_1.6.2             
## [70] evaluate_0.21           leaflet.providers_1.9.0 latticeExtra_0.6-30    
## [73] png_0.1-8               vctrs_0.6.2             tzdb_0.4.0             
## [76] gtable_0.3.3            slippymath_0.3.1        cachem_1.0.8           
## [79] xfun_0.39               mime_0.12               lwgeom_0.2-13          
## [82] e1071_1.7-13            class_7.3-20            viridisLite_0.4.2      
## [85] units_0.8-2             timechange_0.2.0        ellipsis_0.3.2