1 Introducción

Este notebook ilustra dos técnicas de interpolación espacial: Distancia Inversa Ponderada (IDW, por sus siglas en inglés) y Kriging Ordinario (OK). IDW es una técnica determinística; OK es una técnica probabilística basada en el análisis del variograma. Ambas técnicas se usan aquí para obtener una superficie continua del carbono orgánico del suelo (SOC, Soil Organic Carbon) a 15-30 cm de profundidad, a partir de muestras obtenidas de SoilGrids 250 m (ISRIC), para el departamento del Valle del Cauca, Colombia.

Este documento sigue la estructura metodológica propuesta por Lizarazo (2023) para la interpolación de SOC, adaptada a la zona de estudio del Valle del Cauca.

2 Configuración

Primero, se limpia la memoria de trabajo:

rm(list = ls())

Asegúrese de tener instaladas previamente las siguientes librerías:

install.packages(c("terra", "sf", "stars", "gstat", "automap",
                    "sp", "leaflet", "leafem", "dplyr", "ggplot2",
                    "geodata", "RColorBrewer"))

Ahora se cargan las librerías:

library(terra)
library(sf)
library(stars)
library(gstat)
library(automap)
library(sp)
library(leaflet)
library(leafem)
library(dplyr)
library(ggplot2)
library(geodata)
library(RColorBrewer)

Se define una paleta azul-amarillo personalizada, reutilizada en todos los mapas del documento. La paleta "YlGnBu" de RColorBrewer arranca en un amarillo casi blanco; aquí se fuerza un amarillo más intenso y saturado en el extremo inferior:

azul_amarillo <- colorRampPalette(c("#FFEA00", "#ADDD8E", "#41B6C4", "#225EA8", "#0C2C84"))

3 Lectura de los datos de entrada

3.1 Límite administrativo del Valle del Cauca

Se descarga la división administrativa de Colombia (nivel 1 — departamentos) usando el paquete geodata, y se filtra el departamento del Valle del Cauca:

col_dep <- gadm(country = "COL", level = 1, path = tempdir())
valle <- col_dep[col_dep$NAME_1 == "Valle del Cauca", ]
valle_sf <- st_as_sf(valle)
valle_sf
## Simple feature collection with 1 feature and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -77.54486 ymin: 3.060589 xmax: -75.70487 ymax: 5.015217
## Geodetic CRS:  WGS 84
##      GID_1 GID_0  COUNTRY          NAME_1 VARNAME_1 NL_NAME_1       TYPE_1
## 1 COL.31_2   COL Colombia Valle del Cauca     Valle      <NA> Departamento
##    ENGTYPE_1 CC_1 HASC_1  ISO_1                       geometry
## 1 Department <NA>  CO.VC CO-VAC MULTIPOLYGON (((-76.55499 3...

Se calcula el bounding box del departamento, que se usará para centrar todos los mapas leaflet en la zona de estudio (evitando que queden con la vista mundial por defecto):

valle_bbox <- st_bbox(valle_sf)
valle_bbox
##       xmin       ymin       xmax       ymax 
## -77.544861   3.060589 -75.704872   5.015217

3.2 Capa de SOC (SoilGrids 250 m)

Se lee la capa de SOC (15-30 cm) directamente desde el servidor de ISRIC, en formato cloud-optimized GeoTIFF (COG), y se recorta al extent del Valle del Cauca:

url_soc <- "/vsicurl/https://files.isric.org/soilgrids/latest/data/soc/soc_15-30cm_mean.vrt"
soc_col <- rast(url_soc)

valle_wgs84 <- st_transform(valle_sf, crs(soc_col))
soc_valle <- crop(soc_col, vect(valle_wgs84))
soc_valle <- mask(soc_valle, vect(valle_wgs84))
soc_valle
## class       : SpatRaster
## size        : 870, 800, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8636250, -8436250, 340750, 558250  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s)   : memory
## varname     : soc_15-30cm_mean
## name        : soc_15-30cm_mean
## min value   :              132
## max value   :             2529

SoilGrids expresa el SOC en dg/kg con un factor de escala de 10. Se convierte a porcentaje (%):

soc_valle <- soc_valle / 100
names(soc_valle) <- "soc_15_30"
soc_valle
## class       : SpatRaster
## size        : 870, 800, 1  (nrow, ncol, nlyr)
## resolution  : 250, 250  (x, y)
## extent      : -8636250, -8436250, 340750, 558250  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## source(s)   : memory
## varname     : soc_15-30cm_mean
## name        : soc_15_30
## min value   :      1.32
## max value   :     25.29

¿Cuál es el CRS de los datos? SoilGrids se distribuye nativamente en la proyección Homolosine Interrumpida (en metros), y la lectura vía COG no reproyecta nada — solo permite leer el archivo remoto por streaming sin descargarlo completo. Por lo tanto soc_valle sigue en Homolosine, no en WGS84. Verifique el CRS resultante:

crs(soc_valle, describe = TRUE)
##                           name authority code area         extent
## 1 Interrupted_Goode_Homolosine      <NA> <NA> <NA> NA, NA, NA, NA

Se convierte el SpatRaster a un objeto stars para su posterior visualización:

soc_stars <- st_as_stars(soc_valle)
soc_stars
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##            Min. 1st Qu. Median     Mean 3rd Qu.  Max.    NAs
## soc_15_30  1.32    2.98    4.5 6.389164    7.27 25.29 365200
## dimension(s):
##   from  to   offset delta                       refsys x/y
## x    1 800 -8636250   250 Interrupted_Goode_Homolosine [x]
## y    1 870   558250  -250 Interrupted_Goode_Homolosine [y]

3.3 Visualización de la capa de entrada

pal_soc <- colorBin(azul_amarillo(9), domain = c(0, 90), bins = seq(0, 90, 10),
                     na.color = "transparent")

leaflet() %>%
  addProviderTiles("OpenStreetMap") %>%
  addStarsImage(soc_stars, colors = pal_soc, opacity = 0.8) %>%
  addPolygons(data = valle_sf, fill = FALSE, color = "black", weight = 2) %>%
  addLegend(pal = pal_soc, values = c(0, 90), title = "SOC 15-30 cm [%]") %>%
  fitBounds(valle_bbox[["xmin"]], valle_bbox[["ymin"]],
            valle_bbox[["xmax"]], valle_bbox[["ymax"]])

4 Muestreo del área de estudio

Se obtiene una muestra aleatoria de aproximadamente 5000 sitios a partir de los datos “reales” (la capa de SoilGrids), simulando un muestreo de campo:

set.seed(2026)
muestras_v <- spatSample(soc_valle, size = 5000, method = "random",
                          na.rm = TRUE, as.points = TRUE)
muestras_v
## class       : SpatVector
## geometry    : points
## dimensions  : 5000, 1  (geometries, attributes)
## extent      : -8634125, -8436875, 341125, 557875  (xmin, xmax, ymin, ymax)
## coord. ref. : Interrupted_Goode_Homolosine
## names       : soc_15_30
## type        :     <num>
## values      :      6.27
##                    2.75
##                    7.31
##               ...

Se convierte el objeto SpatVector en un objeto sf:

muestras_sf <- st_as_sf(muestras_v)

# IMPORTANTE: soc_valle (y por tanto muestras_v) hereda el CRS nativo de
# SoilGrids (Homolosine Interrumpida, en metros), NO WGS84. Ese CRS nativo
# debe conservarse en la geometría de muestras_sf porque más adelante se usa
# para construir el objeto gstat, que debe coincidir con el CRS de la
# plantilla rrr_stars (también en Homolosine) — de lo contrario predict()
# falla por CRS incompatibles.
#
# Para el mapa leaflet, en cambio, se necesitan coordenadas en grados
# (WGS84), así que se calculan aparte, a partir de una copia reproyectada,
# sin modificar la geometría de muestras_sf:
coords_wgs84 <- st_coordinates(st_transform(muestras_sf, 4326))

muestras_sf <- muestras_sf %>%
  mutate(id = row_number(),
         longit = coords_wgs84[, 1],
         latit  = coords_wgs84[, 2]) %>%
  rename(soc = soc_15_30) %>%
  relocate(id, longit, latit, soc)

muestras_df <- st_drop_geometry(muestras_sf)
head(muestras_df, 10)
##    id    longit    latit   soc
## 1   1 -76.09177 3.286711  6.27
## 2   2 -75.84436 4.667871  2.75
## 3   3 -76.99300 3.416967  7.31
## 4   4 -76.35715 3.733623  2.26
## 5   5 -76.96985 3.535994  6.98
## 6   6 -77.40408 3.248533 13.48
## 7   7 -76.38209 3.466374  1.72
## 8   8 -76.22932 4.113161  2.07
## 9   9 -75.95266 3.951464  5.13
## 10 10 -76.96595 3.419213  6.06

Se eliminan los valores NA, si los hay:

muestras_sf <- muestras_sf %>% filter(!is.na(soc))
nrow(muestras_sf)
## [1] 5000

Visualización de las muestras sobre el mapa:

Con 5000 puntos, círculos pequeños sin borde tienden a mimetizarse con el mapa base. Además, en lugar de dejar que leaflet detecte automáticamente las coordenadas del objeto sf (lo que en algunos entornos no llega a dibujar los puntos), se le pasan explícitamente las columnas longit/latit que ya se calcularon antes:

pal_pts <- colorBin(azul_amarillo(9), domain = c(0, 90), bins = seq(0, 90, 10))

leaflet() %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(data = valle_sf, fill = FALSE, color = "black", weight = 2) %>%
  addCircleMarkers(data = muestras_df, lng = ~longit, lat = ~latit,
                    radius = 4, stroke = TRUE, color = "black", weight = 0.4,
                    fillOpacity = 1, fillColor = ~pal_pts(soc)) %>%
  addLegend(pal = pal_pts, values = c(0, 90), title = "SOC muestreado [%]") %>%
  fitBounds(valle_bbox[["xmin"]], valle_bbox[["ymin"]],
            valle_bbox[["xmax"]], valle_bbox[["ymax"]])

Con las muestras listas, se procede con las tareas de interpolación.

5 Interpolación

5.1 Creación del objeto gstat

Para interpolar se crea primero un objeto de clase gstat, que contiene:

  • La definición del modelo (fórmula de predicción)
  • Los datos de calibración
  • El modelo de variograma (si aplica)

Sin modelo de variograma → IDW. Con modelo de variograma y sin covariables → Kriging Ordinario.

5.2 Interpolación IDW

Se crea el objeto gstat para IDW, especificando solo la fórmula y los datos:

muestras_sp <- as(muestras_sf, "Spatial")

g1 <- gstat(formula = soc ~ 1, data = muestras_sp, nmax = 50)
g1
## data:
## var1 : formula = soc`~`1 ; data dim = 5000 x 4 nmax = 50

nmax = 50 indica que cada predicción usa solo los 50 vecinos más cercanos (vecindario local), en lugar de las 5000 muestras completas. Es una práctica estándar en interpolación local y acelera notablemente el cálculo sin afectar de forma perceptible el resultado, ya que los puntos lejanos aportan pesos prácticamente nulos en IDW.

Se crea una plantilla ráster con celdas de valor 1, sobre la cual se harán las predicciones. La capa soc_valle viene a la resolución nativa de SoilGrids (250 m), lo cual para un departamento entero implica cientos de miles de celdas; predecir en cada una de ellas —sobre todo con Kriging, que además de un promedio ponderado calcula la varianza resolviendo un sistema de ecuaciones— es innecesariamente lento para un mapa a escala departamental. Por eso se agrega (aggregate) el ráster a una resolución más gruesa (~1.25 km, factor 5) antes de predecir:

rrr <- aggregate(soc_valle, fact = 5, na.rm = TRUE)
values(rrr) <- 1
names(rrr) <- "valor"
rrr_stars <- st_as_stars(rrr)
rrr_stars
## stars object with 2 dimensions and 1 attribute
## attribute(s):
##        Min. 1st Qu. Median Mean 3rd Qu. Max.
## valor     1       1      1    1       1    1
## dimension(s):
##   from  to   offset delta                       refsys x/y
## x    1 160 -8636250  1250 Interrupted_Goode_Homolosine [x]
## y    1 174   558250 -1250 Interrupted_Goode_Homolosine [y]

Si necesita un mapa final a mayor detalle, puede reducir fact (p. ej. fact = 2) una vez que confirme que el resto del documento corre en un tiempo razonable.

Se interpola el SOC según el modelo g1 sobre la plantilla rrr_stars:

z1 <- predict(g1, rrr_stars)
## [inverse distance weighted interpolation]
z1
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                Min.  1st Qu.   Median    Mean  3rd Qu.     Max.   NAs     NA's
## var1.pred  1.526434 3.920044 6.132754 8.20181 14.01582 21.56844     0 1.526434
## var1.var         NA       NA       NA     NaN       NA       NA 27840 0.000000
## dimension(s):
##   from  to   offset delta                       refsys x/y
## x    1 160 -8636250  1250 Interrupted_Goode_Homolosine [x]
## y    1 174   558250 -1250 Interrupted_Goode_Homolosine [y]

Se extrae el primer atributo (var1.pred) y se renombra a soc:

idw_soc <- z1["var1.pred"]
names(idw_soc) <- "soc"

La plantilla rrr_stars se construyó sobre el bounding box rectangular del departamento (no sobre su forma real), por lo que la predicción idw_soc también sale como un rectángulo. Se enmascara el resultado con el polígono del Valle del Cauca para recortarlo a su forma correcta:

idw_soc <- st_as_stars(mask(rast(idw_soc), vect(valle_wgs84)))
names(idw_soc) <- "soc"

5.3 Interpolación OK (Kriging Ordinario)

Los métodos de Kriging requieren un modelo de variograma. El variograma cuantifica el patrón de autocorrelación espacial de los datos y asigna pesos en función de él.

Cálculo del variograma empírico (sin covariables):

v_emp <- variogram(soc ~ 1, data = muestras_sf)
plot(v_emp)

Ajuste automático del modelo de variograma con automap::autofitVariogram. Esta función no trabaja directamente con objetos sf, por lo que se convierte a Spatial:

v_fit <- autofitVariogram(soc ~ 1, input_data = muestras_sp)
plot(v_fit)

Parámetros del modelo ajustado:

v_fit$var_model
##   model      psill    range kappa
## 1   Nug   1.938418      0.0   0.0
## 2   Ste 190.723845 668125.6   0.6
  • model: tipo de modelo teórico ajustado (p. ej. Nugget + Stein/Matérn, esférico, exponencial).
  • psill: meseta parcial (varianza explicada por cada componente del modelo).
  • range: rango o distancia a partir de la cual las observaciones dejan de estar autocorrelacionadas.
  • kappa: parámetro de forma (solo aplica a modelos tipo Matérn/Stein).
  • ang1/ang2/ang3, anis1/anis2: parámetros de anisotropía (dirección e intensidad), iguales a 0/1 cuando el fenómeno se asume isotrópico.

Con el modelo de variograma ajustado, se realiza el Kriging Ordinario. Al igual que con IDW, se limita el vecindario a nmax = 50 — esto es especialmente importante en Kriging, ya que sin este límite gstat resolvería, para cada una de las celdas del ráster de salida, un sistema de ecuaciones con las 5000 muestras completas, lo cual sería extremadamente lento (es la causa más probable de que el knit se quede colgado en este chunk). Con vecindario local, cada predicción solo resuelve un sistema de 50×50:

g2 <- gstat(formula = soc ~ 1, data = muestras_sp, model = v_fit$var_model, nmax = 50)
z2 <- predict(g2, rrr_stars)
## [using ordinary kriging]
z2
## stars object with 2 dimensions and 2 attributes
## attribute(s):
##                     Min.  1st Qu.   Median     Mean   3rd Qu.     Max.
## var1.pred   1.661121e+00 3.986098 6.664645 8.214417 13.156002 19.27633
## var1.var   -1.989520e-13 2.342227 2.672893 7.825055  8.579866 54.29842
## dimension(s):
##   from  to   offset delta                       refsys x/y
## x    1 160 -8636250  1250 Interrupted_Goode_Homolosine [x]
## y    1 174   558250 -1250 Interrupted_Goode_Homolosine [y]

Se extrae y renombra el atributo de predicción:

ok_soc <- z2["var1.pred"]
names(ok_soc) <- "soc"

Igual que con IDW, se enmascara el resultado a la forma real del departamento:

ok_soc <- st_as_stars(mask(rast(ok_soc), vect(valle_wgs84)))
names(ok_soc) <- "soc"

5.4 Mapas de la interpolación

En el mapa anterior el rango de color se había fijado de antemano en 0-90 %. El problema es que IDW y, sobre todo, Kriging pueden extrapolar valores fuera del rango observado en los bordes del área (incluso negativos), y esas celdas quedaban fuera del dominio fijo — se pintaban transparentes, por lo que el mapa mostraba el borde del departamento pero sin relleno de color. Para evitarlo, el dominio de la paleta se calcula ahora a partir del rango real combinado de los datos “reales”, IDW y OK:

rango <- range(c(values(soc_valle), idw_soc[[1]], ok_soc[[1]]), na.rm = TRUE)
cortes <- pretty(rango, n = 8)
pal_res <- colorBin(azul_amarillo(length(cortes) - 1), domain = rango, bins = cortes, na.color = "transparent")
rango
## [1]  1.32 25.29

Mapa de la interpolación IDW:

leaflet() %>%
  addProviderTiles("OpenStreetMap") %>%
  addStarsImage(idw_soc, colors = pal_res, opacity = 0.85) %>%
  addPolygons(data = valle_sf, fill = FALSE, color = "black", weight = 2) %>%
  addLegend(pal = pal_res, values = rango, title = "IDW SOC [%]") %>%
  fitBounds(valle_bbox[["xmin"]], valle_bbox[["ymin"]],
            valle_bbox[["xmax"]], valle_bbox[["ymax"]])

Mapa de la interpolación OK:

leaflet() %>%
  addProviderTiles("OpenStreetMap") %>%
  addStarsImage(ok_soc, colors = pal_res, opacity = 0.85) %>%
  addPolygons(data = valle_sf, fill = FALSE, color = "black", weight = 2) %>%
  addLegend(pal = pal_res, values = rango, title = "OK SOC [%]") %>%
  fitBounds(valle_bbox[["xmin"]], valle_bbox[["ymin"]],
            valle_bbox[["xmax"]], valle_bbox[["ymax"]])

6 Evaluación de resultados

6.1 Evaluación cualitativa

Comparación visual de los tres productos: datos “reales” (SoilGrids), IDW y OK.

real_stars <- soc_stars
names(real_stars) <- "soc"

leaflet() %>%
  addProviderTiles("OpenStreetMap") %>%
  addStarsImage(real_stars, colors = pal_res, opacity = 0.85, group = "Real") %>%
  addStarsImage(idw_soc, colors = pal_res, opacity = 0.85, group = "IDW") %>%
  addStarsImage(ok_soc, colors = pal_res, opacity = 0.85, group = "OK") %>%
  addPolygons(data = valle_sf, fill = FALSE, color = "black", weight = 2) %>%
  addLayersControl(baseGroups = c("Real", "IDW", "OK"),
                    options = layersControlOptions(collapsed = FALSE)) %>%
  addLegend(pal = pal_res, values = rango, title = "SOC [%]") %>%
  fitBounds(valle_bbox[["xmin"]], valle_bbox[["ymin"]],
            valle_bbox[["xmax"]], valle_bbox[["ymax"]])

6.2 Validación cruzada

Con 5000 muestras, la validación cruzada dejando-uno-fuera (Leave-One-Out, el valor por defecto de gstat.cv) exige 5000 predicciones individuales — cada una con su propia búsqueda de vecinos — y es la causa más probable de que el knit se quede pegado en este punto. En su lugar se usa validación cruzada de 10 particiones (nfold = 10): los datos se dividen en 10 grupos, y por turnos cada grupo se predice de una sola vez (vectorizado) usando los otros 9 como calibración. Es la práctica estándar para este tamaño de muestra y da una estimación del error igual de confiable, en una fracción del tiempo:

set.seed(2026)
cv1 <- gstat.cv(g1, nfold = 10, verbose = FALSE)
bubble(cv1, "residual", main = "Residuales IDW")

rmse_idw <- sqrt(mean(cv1$residual^2, na.rm = TRUE))
rmse_idw
## [1] 1.405091

Se repite el proceso para los resultados de OK:

set.seed(2026)
cv2 <- gstat.cv(g2, nfold = 10, verbose = FALSE)
bubble(cv2, "residual", main = "Residuales OK")

rmse_ok <- sqrt(mean(cv2$residual^2, na.rm = TRUE))
rmse_ok
## [1] 1.511265
data.frame(metodo = c("IDW", "Kriging Ordinario"),
           RMSE = c(rmse_idw, rmse_ok))
##              metodo     RMSE
## 1               IDW 1.405091
## 2 Kriging Ordinario 1.511265

Resumen de resultados: compare los valores de RMSE obtenidos para IDW y OK. El método con menor RMSE es, en términos de esta validación cruzada, el más preciso para representar la variabilidad espacial del carbono orgánico del suelo en el Valle del Cauca. Comente además si el modelo de variograma ajustado (rango, meseta) es coherente con la escala de variabilidad edáfica esperada para la zona (valle interandino, piedemonte y zona montañosa del departamento), y qué limitaciones tiene usar datos de SoilGrids como “verdad de campo” en lugar de mediciones directas de laboratorio.

7 Referencias

  • Hengl, T., de Jesus, J. M., Heuvelink, G. B. M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., … & Kempen, B. (2017). SoilGrids250m: Global gridded soil information based on machine learning. PLOS ONE, 12(2), e0169748. https://doi.org/10.1371/journal.pone.0169748
  • ISRIC — World Soil Information. SoilGrids REST API y datos raster. https://www.isric.org/explore/soilgrids
  • Hijmans, R. J. (2023). terra: Spatial Data Analysis. R package. https://CRAN.R-project.org/package=terra
  • Pebesma, E. (2018). Simple Features for R: Standardized Support for Spatial Vector Data. The R Journal, 10(1), 439-446.
  • Pebesma, E., & Graeler, B. (2023). gstat: Spatial and Spatio-Temporal Geostatistical Modelling, Prediction and Simulation. R package.
  • Hiemstra, P. H., Pebesma, E. J., Twenhöfel, C. J. W., & Heuvelink, G. B. M. (2009). Real-time automatic interpolation of ambient gamma dose rates from the Dutch radioactivity monitoring network. Computers & Geosciences, 35(8), 1711-1721. (paquete automap)
  • Hijmans, R. J. (2023). geodata: Download Geographic Data. R package. https://CRAN.R-project.org/package=geodata
  • Lizarazo, I. (2023). Spatial interpolation of soil organic carbon. Disponible en https://rpubs.com/ials2un/soc_interp
sessionInfo()
## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Spain.utf8  LC_CTYPE=Spanish_Spain.utf8   
## [3] LC_MONETARY=Spanish_Spain.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Spain.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3 geodata_0.6-9      ggplot2_4.0.3      dplyr_1.2.1       
##  [5] leafem_0.2.5       leaflet_2.2.3      sp_2.2-1           automap_1.1-20    
##  [9] gstat_2.1-6        stars_0.7-2        abind_1.4-8        sf_1.1-1          
## [13] terra_1.9-27      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6            xfun_0.57               bslib_0.10.0           
##  [4] raster_3.6-32           htmlwidgets_1.6.4       lattice_0.22-9         
##  [7] leaflet.providers_3.0.0 vctrs_0.7.3             tools_4.6.0            
## [10] crosstalk_1.2.2         generics_0.1.4          parallel_4.6.0         
## [13] tibble_3.3.1            proxy_0.4-29            spacetime_1.3-3        
## [16] xts_0.14.2              pkgconfig_2.0.3         KernSmooth_2.23-26     
## [19] S7_0.2.2                lifecycle_1.0.5         compiler_4.6.0         
## [22] farver_2.1.2            FNN_1.1.4.1             codetools_0.2-20       
## [25] htmltools_0.5.9         class_7.3-23            sass_0.4.10            
## [28] yaml_2.3.12             pillar_1.11.1           jquerylib_0.1.4        
## [31] classInt_0.4-11         cachem_1.1.0            tidyselect_1.2.1       
## [34] digest_0.6.39           fastmap_1.2.0           grid_4.6.0             
## [37] cli_3.6.6               magrittr_2.0.5          base64enc_0.1-6        
## [40] e1071_1.7-17            withr_3.0.2             scales_1.4.0           
## [43] rappdirs_0.3.4          rmarkdown_2.31          otel_0.2.0             
## [46] zoo_1.8-15              png_0.1-9               evaluate_1.0.5         
## [49] knitr_1.51              rlang_1.2.0             Rcpp_1.1.1-1.1         
## [52] glue_1.8.1              DBI_1.3.0               rstudioapi_0.18.0      
## [55] reshape_0.8.10          jsonlite_2.0.0          R6_2.6.1               
## [58] plyr_1.8.9              intervals_0.15.5        units_1.0-1