M2U1 - Caso: Análisis Geoestadístico de Clima en Finca de Aguacate (Popayán)

options(repos = c(CRAN = "https://cloud.r-project.org"))

paquetes <- c("readxl", "dplyr", "geoR", "raster", "rasterVis",
              "ggplot2", "leaflet", "sp", "knitr", "kableExtra")

for (pkg in paquetes) {
  if (!require(pkg, character.only = TRUE, quietly = TRUE)) {
    install.packages(pkg, dependencies = TRUE,
                     repos = "https://cloud.r-project.org")
    library(pkg, character.only = TRUE)
  }
}

paleta <- c("#2C7BB6", "#D7191C", "#1A9641", "#FDAE61", "#ABD9E9")

1 Introducción

La temperatura incide directamente en la productividad del aguacate. Conocer su distribución espacial dentro de la finca permite decisiones de manejo más precisas, pero los sensores solo capturan valores en puntos discretos, haciendo necesaria la interpolación para obtener una superficie continua.

Este informe aplica geoestadística a 534 mediciones de temperatura georeferenciadas por árbol, tomadas el 1 de octubre de 2020 en una finca aguacatera de Popayán, Cauca. El objetivo es estimar la distribución continua de temperatura mediante Kriging Ordinario, aprovechando la autocorrelación espacial de los datos para producir predicciones insesgadas y de varianza mínima. El análisis sigue el flujo estándar: exploración de datos, semivariograma empírico, ajuste del modelo teórico, interpolación espacial, generación de rasters y validación cruzada.

Los resultados obtenidos muestran la magnitud real de la variación térmica entre árboles en una misma jornada, identifican los sectores más fríos y más cálidos dentro del lote y señalan concretamente dónde la red de sensores actual deja zonas con predicciones poco confiables. Cabe anotar que el análisis cubre un único día de medición, por lo que las conclusiones reflejan las condiciones de ese momento y no necesariamente el patrón térmico permanente de la finca.

2 Análisis Exploratorio de Datos

geo_datos_raw <- read_excel(
  "C:/Pontificia_Un_Javeriana/Analisis_Geoespacial/Datos_Completos_Aguacate.xlsx"
)

cat("Dimensiones totales:", dim(geo_datos_raw), "\n")

## Dimensiones totales: 20271 21

cat("Columnas:", paste(names(geo_datos_raw), collapse = ", "), "\n")

## Columnas: id_arbol, Latitude, Longitude, FORMATTED_DATE_TIME, Psychro_Wet_Bulb_Temperature, Station_Pressure, Relative_Humidity, Crosswind, Temperature, Barometric_Pressure, Headwind, Direction_True, Direction_Mag, Wind_Speed, Heat_Stress_Index, Altitude, Dew_Point, Density_Altitude, Wind_Chill, Estado_Fenologico_Predominante, Frutos_Afectados

cat("Tipo FORMATTED_DATE_TIME:", class(geo_datos_raw$FORMATTED_DATE_TIME), "\n\n")

## Tipo FORMATTED_DATE_TIME: character

# Filtrado robusto: prueba 4 estrategias en cascada según tipo de columna
fecha_obj <- as.Date("2020-10-01")
vals <- geo_datos_raw$FORMATTED_DATE_TIME

if (inherits(vals, c("Date","POSIXct","POSIXlt"))) {
  idx <- !is.na(as.Date(vals)) & as.Date(vals) == fecha_obj
} else {
  txt <- as.character(vals)
  # intento dd/mm/yyyy
  p1 <- suppressWarnings(as.Date(txt, format = "%d/%m/%Y"))
  if (any(!is.na(p1))) {
    idx <- !is.na(p1) & p1 == fecha_obj
  } else {
    # intento yyyy-mm-dd (primeros 10 chars)
    p2 <- suppressWarnings(as.Date(substr(txt, 1, 10)))
    if (any(!is.na(p2))) {
      idx <- !is.na(p2) & p2 == fecha_obj
    } else {
      # búsqueda de texto
      idx <- grepl("2020-10-01|01/10/2020|01-10-2020", txt)
    }
  }
}

geo_datos <- geo_datos_raw[idx, ]

if (nrow(geo_datos) == 0) {
  warning("Filtro de fecha sin resultados. Se usarán todos los datos.")
  geo_datos <- geo_datos_raw
}

cat("Registros filtrados:", nrow(geo_datos), "\n")

## Registros filtrados: 534

cat("Rango Longitude:", round(range(as.double(geo_datos$Longitude), na.rm=TRUE), 6), "\n")

## Rango Longitude: -76.7118 -76.71022

cat("Rango Latitude :", round(range(as.double(geo_datos$Latitude),  na.rm=TRUE), 6), "\n")

## Rango Latitude : 2.392101 2.393634

# Variables inline del dataset
res_temp_min   <- round(min(geo_datos$Temperature,  na.rm = TRUE), 1)
res_temp_max   <- round(max(geo_datos$Temperature,  na.rm = TRUE), 1)
res_temp_media <- round(mean(geo_datos$Temperature, na.rm = TRUE), 1)
res_temp_sd    <- round(sd(geo_datos$Temperature,   na.rm = TRUE), 1)
res_temp_cv    <- round(sd(geo_datos$Temperature,   na.rm = TRUE) /
                         mean(geo_datos$Temperature, na.rm = TRUE) * 100, 1)
res_temp_rango <- round(res_temp_max - res_temp_min, 1)
res_lon_rng_km <- round(diff(range(as.double(geo_datos$Longitude), na.rm = TRUE)) * 111, 2)
res_lat_rng_km <- round(diff(range(as.double(geo_datos$Latitude),  na.rm = TRUE)) * 111, 2)
res_moda_temp  <- round(as.numeric(names(which.max(table(round(geo_datos$Temperature, 0))))), 0)

kable(head(geo_datos[, c("id_arbol","Latitude","Longitude",
                          "Temperature","Relative_Humidity",
                          "Wind_Speed","Altitude")], 10),
      caption = "Tabla 1. Primeros 10 registros del conjunto de datos filtrado al período 01/10/2020",
      digits = 3, align = "c") %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = FALSE)

Tabla 1. Primeros 10 registros del conjunto de datos filtrado al período 01/10/2020

id_arbol	Latitude	Longitude	Temperature	Relative_Humidity	Wind_Speed	Altitude
1	2.394	-76.711	23.9	85.2	0.0	1696
2	2.394	-76.711	23.5	84.0	0.0	1696
3	2.394	-76.711	24.5	79.6	0.5	1694
4	2.394	-76.711	25.9	77.6	0.5	1694
5	2.393	-76.711	26.0	76.5	0.0	1696
6	2.393	-76.711	24.5	77.7	0.3	1698
7	2.393	-76.711	25.5	76.5	0.0	1698
8	2.393	-76.711	25.7	77.7	0.0	1698
9	2.393	-76.711	26.0	78.3	0.0	1696
10	2.393	-76.711	25.9	77.8	0.4	1694

En esta sección se realiza el análisis exploratorio de la Temperatura mediante estadísticas descriptivas y visualizaciones.

2.1 Estadísticas Descriptivas

resumen_temp <- data.frame(
  Estadístico = c("Mínimo","Percentil 25","Mediana","Media",
                  "Percentil 75","Máximo","Desv. Estándar","CV (%)"),
  Valor = c(
    round(min(geo_datos$Temperature,                     na.rm=TRUE), 3),
    round(quantile(geo_datos$Temperature, 0.25,          na.rm=TRUE), 3),
    round(median(geo_datos$Temperature,                  na.rm=TRUE), 3),
    round(mean(geo_datos$Temperature,                    na.rm=TRUE), 3),
    round(quantile(geo_datos$Temperature, 0.75,          na.rm=TRUE), 3),
    round(max(geo_datos$Temperature,                     na.rm=TRUE), 3),
    round(sd(geo_datos$Temperature,                      na.rm=TRUE), 3),
    round(sd(geo_datos$Temperature,   na.rm=TRUE) /
            mean(geo_datos$Temperature, na.rm=TRUE) * 100, 2)
  )
)

kable(resumen_temp,
      caption = "Tabla 2. Estadísticas descriptivas de la variable Temperatura (°C)",
      align = c("l","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Tabla 2. Estadísticas descriptivas de la variable Temperatura (°C)

Estadístico	Valor
Mínimo	22.200
Percentil 25	24.500
Mediana	25.800
Media	25.829
Percentil 75	27.175
Máximo	29.700
Desv. Estándar	1.771
CV (%)	6.860

Interpretación: La tabla muestra que la temperatura registrada el 1 de octubre de 2020 osciló entre 22.2 °C y 29.7 °C, con una media de 25.8 °C y una desviación estándar de 1.8 °C. El coeficiente de variación de 6.9 % indica una dispersión moderada que sugiere variabilidad espacial detectable. Ese rango de aproximadamente 7.5 °C entre el árbol más frío y el más caliente es agronómicamente relevante, ya que diferencias superiores a 2 °C dentro de un lote pueden afectar la sincronía de floración del aguacate Hass.

2.2 Distribución de Temperatura

ggplot(geo_datos, aes(x = Temperature)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 25, fill = paleta[1], color = "white", alpha = 0.85) +
  geom_density(color = paleta[2], linewidth = 1.2) +
  labs(x = "Temperatura (°C)", y = "Densidad") +
  theme_minimal(base_size = 13) +
  theme(plot.title    = element_text(face = "bold", color = paleta[1]))

Figura 1. Distribución de frecuencias de la temperatura registrada en los árboles de aguacate durante el período 01/10/2020.

Interpretación: El histograma muestra una distribución con dos modos aproximados alrededor de 26 °C y en el rango de 26-27 °C, lo que sugiere que en la finca coexisten al menos dos microambientes térmicos diferenciados. La asimetría leve hacia valores altos indica que un grupo de árboles registró temperaturas notoriamente superiores a la mayoría. Dado que la distribución no se aleja demasiado de la normalidad, el Kriging Ordinario puede aplicarse sin transformación previa de la variable, aunque sería prudente verificar este supuesto formalmente con una prueba de normalidad si el análisis se extiende a otros períodos.

2.3 Distribución en Árboles

pal_temp <- colorNumeric(palette = c(paleta[5], paleta[4], paleta[2]),
                         domain  = geo_datos$Temperature)

leaflet(geo_datos) %>%
  addTiles() %>%
  addCircleMarkers(lng = ~Longitude, lat = ~Latitude, radius = 3,
                   color = ~pal_temp(Temperature),
                   stroke = FALSE, fillOpacity = 0.85,
                   popup = ~paste0("<b>Árbol:</b> ", id_arbol,
                                   "<br><b>Temp:</b> ", round(Temperature,2), " °C")) %>%
  addLegend(position = "bottomright", pal = pal_temp,
            values = ~Temperature, title = "Temperatura (°C)", opacity = 0.85)

Figura 2. Mapa interactivo de distribución espacial de los árboles de aguacate con código de color por temperatura (Popayán, Cauca).

Interpretación: El mapa confirma que los 534 árboles monitoreados cubren de forma continua la extensión de la finca, con una distribución que abarca aproximadamente 0.18 km en longitud y 0.17 km en latitud. Se observa visualmente que los árboles con temperatura más alta se concentran en un sector particular del lote, lo que podría estar relacionado con la exposición solar, el drenaje del suelo o la densidad del dosel en esa zona. Esta concentración espacial de valores extremos es precisamente lo que el semivariograma deberá capturar.

2.4 Temperatura y Coordenadas

ggplot(geo_datos, aes(x = Latitude, y = Temperature)) +
  geom_point(color = paleta[1], alpha = 0.5, size = 1.8) +
  geom_smooth(method = "lm", formula = y ~ x,
              color = paleta[2], se = TRUE, linewidth = 1.2,
              fill = paleta[5], alpha = 0.3) +
  labs(x = "Latitud", y = "Temperatura (°C)") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", color = paleta[1]))

Figura 3. Dispersión de la temperatura en función de la latitud. La línea de regresión permite evaluar tendencia espacial norte-sur.

Interpretación: La nube de puntos contra la latitud muestra una dispersión amplia sin un patrón lineal evidente, lo que sugiere que no existe una tendencia térmica sistemática en dirección norte-sur dentro del lote. La banda de confianza de la regresión lineal es relativamente ancha, reforzando que la pendiente no es estadísticamente robusta. Esto es consistente con la topografía relativamente plana de la zona baja de Popayán donde se ubica la finca. La ausencia de tendencia norte-sur es un indicio favorable para asumir estacionariedad intrínseca y proceder con Kriging Ordinario.

ggplot(geo_datos, aes(x = Longitude, y = Temperature)) +
  geom_point(color = paleta[3], alpha = 0.5, size = 1.8) +
  geom_smooth(method = "lm", formula = y ~ x,
              color = paleta[2], se = TRUE, linewidth = 1.2,
              fill = paleta[5], alpha = 0.3) +
  labs(x = "Longitud", y = "Temperatura (°C)") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(face = "bold", color = paleta[3]))

Figura 4. Dispersión de la temperatura en función de la longitud. La línea de regresión permite evaluar tendencia espacial este-oeste.

Interpretación: La dispersión contra la longitud tampoco evidencia una tendencia lineal clara en dirección este-oeste. Sin embargo, se aprecia una mayor densidad de puntos con temperatura elevada en el extremo occidental del lote, lo que podría corresponder a la parte de la finca con mayor exposición a la tarde. Si este patrón se confirma con mediciones de otros días, sería un argumento para explorar Kriging Universal con una tendencia lineal en longitud. Por ahora, con un único día de datos, el Kriging Ordinario es la elección más parsimoniosa.

3 Evaluación con Geodata

En esta sección se construye el objeto geodata de geoR y se visualizan los cuatro paneles diagnóstico equivalentes.

# Extraer vectores double puros y eliminar NA/Inf
lon_vec  <- as.double(geo_datos[["Longitude"]])
lat_vec  <- as.double(geo_datos[["Latitude"]])
temp_vec <- as.double(geo_datos[["Temperature"]])

idx_ok   <- is.finite(lon_vec) & is.finite(lat_vec) & is.finite(temp_vec)
lon_vec  <- lon_vec[idx_ok]
lat_vec  <- lat_vec[idx_ok]
temp_vec <- temp_vec[idx_ok]

cat("Observaciones limpias:", length(temp_vec), "\n")

## Observaciones limpias: 534

cat("Rango Longitud       :", round(range(lon_vec), 6), "\n")

## Rango Longitud       : -76.7118 -76.71022

cat("Rango Latitud        :", round(range(lat_vec), 6), "\n")

## Rango Latitud        : 2.392101 2.393634

cat("Rango Temperatura    :", round(range(temp_vec), 3), "\n")

## Rango Temperatura    : 22.2 29.7

# Construir geodata manualmente (evita problemas con tibbles de readxl)
geo_temp <- list(
  coords = matrix(c(lon_vec, lat_vec), ncol = 2,
                  dimnames = list(NULL, c("lon","lat"))),
  data   = temp_vec
)
class(geo_temp) <- "geodata"
res_n <- length(geo_temp$data)

# Panel 1: distribución espacial por cuartil (equivalente a plot.geodata panel 1)
q_num <- dplyr::ntile(geo_temp$data, 4)
lbl_map <- c("1"="Q1 Bajo","2"="Q2","3"="Q3","4"="Q4 Alto")

df_geo <- data.frame(
  lon     = geo_temp$coords[,"lon"],
  lat     = geo_temp$coords[,"lat"],
  temp    = geo_temp$data,
  cuartil = factor(lbl_map[as.character(q_num)], levels = lbl_map)
)

ggplot(df_geo, aes(x = lon, y = lat, color = cuartil)) +
  geom_point(size = 1.5, alpha = 0.8) +
  scale_color_manual(values = c("#2C7BB6","#1A9641","#FDAE61","#D7191C"),
                     name = "Cuartil", drop = FALSE) +
  labs(x = "Longitud", y = "Latitud") +
  theme_minimal(base_size = 12) + coord_fixed()

Figura 5. Distribución espacial de los puntos de muestreo clasificados por cuartil de temperatura.

ggplot(df_geo, aes(x = lat, y = temp)) +
  geom_point(color = paleta[1], alpha = 0.4, size = 1.5) +
  geom_smooth(method = "loess", formula = y ~ x,
              color = paleta[2], se = FALSE, linewidth = 1) +
  labs(x = "Latitud", y = "Temperatura (°C)") +
  theme_minimal(base_size = 12)

Figura 6. Temperatura en función de la latitud con suavizado LOESS. Panel diagnóstico de tendencia norte-sur del geodata.

ggplot(df_geo, aes(x = lon, y = temp)) +
  geom_point(color = paleta[3], alpha = 0.4, size = 1.5) +
  geom_smooth(method = "loess", formula = y ~ x,
              color = paleta[2], se = FALSE, linewidth = 1) +
  labs(x = "Longitud", y = "Temperatura (°C)") +
  theme_minimal(base_size = 12)

Figura 7. Temperatura en función de la longitud con suavizado LOESS. Panel diagnóstico de tendencia este-oeste del geodata.

ggplot(df_geo, aes(x = temp)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 20, fill = paleta[1], color = "white", alpha = 0.8) +
  geom_density(color = paleta[2], linewidth = 1.2) +
  labs(x = "Temperatura (°C)", y = "Densidad") +
  theme_minimal(base_size = 12)

Figura 8. Distribución marginal de la temperatura con histograma y curva de densidad kernel.

Interpretación: El panel de cuartiles (Figura 5) muestra que los árboles del Q4 (temperaturas más altas) no están distribuidos aleatoriamente: forman un cluster visible en el sector noroccidental del lote. Los árboles del Q1 (más fríos) tienden a ubicarse en el sector suroriental. Esta agrupación espacial de valores similares es la señal más directa de que existe autocorrelación espacial positiva en la temperatura, lo que justifica el uso de geoestadística frente a métodos de interpolación determinísticos. Los paneles LOESS contra latitud y longitud confirman lo observado en las regresiones lineales: no hay una tendencia global dominante, aunque la curva contra latitud muestra una ligera curvatura que merece atención. El histograma final (Figura 8) reitera la bimodalidad observada previamente.

4 Semivariograma Empírico

En esta sección se calcula el semivariograma empírico y los envolventes de Monte Carlo para evaluar la autocorrelación espacial de la temperatura.

cat("Observaciones en geodata:", length(geo_temp$data), "\n")

## Observaciones en geodata: 534

if (length(geo_temp$data) < 10)
  stop("Geodata con menos de 10 obs. Revisar filtrado de fecha.")

# max.dist = mitad de la diagonal del bounding box
delta_lon <- diff(range(geo_temp$coords[,"lon"]))
delta_lat <- diff(range(geo_temp$coords[,"lat"]))
max_dist  <- sqrt(delta_lon^2 + delta_lat^2) * 0.5

# Protección si la finca tiene extensión casi nula en alguna dirección
if (!is.finite(max_dist) || max_dist < 1e-10)
  max_dist <- max(delta_lon, delta_lat, na.rm = TRUE) * 0.5

cat("max_dist para variog:", round(max_dist, 8), "| finito:", is.finite(max_dist), "\n")

## max_dist para variog: 0.00110217 | finito: TRUE

variograma <- variog(geo_temp, option = "bin", max.dist = max_dist)

## variog: computing omnidirectional variogram

# dist_50 reutilizado en ajuste-modelos para ini.vals
dist_50 <- max_dist

cat("Bins calculados:", length(variograma$u), "\n")

## Bins calculados: 13

set.seed(params$seed)
variograma_mc <- variog.mc.env(geo_temp, obj.variog = variograma, nsim = 99)

## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops

cat("Envolventes Monte Carlo: 99 simulaciones completadas.\n")

## Envolventes Monte Carlo: 99 simulaciones completadas.

plot(variograma,
     xlab = "Distancia (grados)", ylab = "Semivarianza",
     pch = 16, col = paleta[1])
lines(variograma_mc, col = "grey60", lty = 2)

Figura 9. Semivariograma empírico omnidireccional de la temperatura con envolventes de Monte Carlo (99 simulaciones).

Interpretación: El semivariograma empírico calculado sobre 534 observaciones con una distancia máxima de 0.0011 grados muestra un incremento claro de la semivarianza en las primeras clases de distancia, lo que confirma que los árboles cercanos entre sí registraron temperaturas más parecidas que los distantes. Los puntos del semivariograma que superan el envolvente superior de Monte Carlo en las distancias cortas indican que esta autocorrelación no es producto del azar. La meseta visual del semivariograma se alcanza aproximadamente a la mitad del rango analizado, lo que implica que más allá de esa distancia los árboles ya no se influyen térmicamente entre sí. Una limitación importante es que el semivariograma se calculó en modo omnidireccional; si la autocorrelación fuera anisotrópica (diferente según la dirección), este análisis la subestimaría.

5 Ajuste de Modelos

En esta sección se ajustan tres modelos teóricos al semivariograma empírico y se selecciona el de menor SCE.

# El argumento correcto de variofit() es ini.cov.pars (NO ini.vals).
# ini.cov.pars acepta: vector c(sigmasq, phi), o matriz/data.frame con
# 2 columnas donde col1=sigmasq y col2=phi.
# expand.grid() produce data.frame, que variofit convierte internamente a matrix.
sigmasq_seq <- seq(var(geo_temp$data) * 0.3, var(geo_temp$data) * 2, length.out = 8)
phi_seq     <- seq(dist_50 * 0.1, dist_50 * 0.8, length.out = 8)
ini_grid    <- expand.grid(sigmasq_seq, phi_seq)   # data.frame: col1=sigmasq, col2=phi

model_exp  <- variofit(variograma,
                       ini.cov.pars = ini_grid,
                       cov.model    = "exponential",
                       weights      = "npairs")

## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.23"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 5720.67737872684

cat("Exponential SCE:", round(model_exp$value, 4), "\n")

## Exponential SCE: 5678.797

model_gaus <- variofit(variograma,
                       ini.cov.pars = ini_grid,
                       cov.model    = "gaussian",
                       weights      = "npairs")

## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.23"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 11370.2209167688

cat("Gaussian    SCE:", round(model_gaus$value, 4), "\n")

## Gaussian    SCE: 10284.05

model_spe  <- variofit(variograma,
                       ini.cov.pars = ini_grid,
                       cov.model    = "spherical",
                       weights      = "npairs")

## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "3.23"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 10020.0397796141

cat("Spherical   SCE:", round(model_spe$value, 4), "\n")

## Spherical   SCE: 8927.608

tabla_modelos <- data.frame(
  Modelo    = c("Exponential","Gaussian","Spherical"),
  SCE       = round(c(model_exp$value,       model_gaus$value,       model_spe$value),  2),
  Nugget    = round(c(model_exp$nugget,       model_gaus$nugget,      model_spe$nugget), 4),
  Meseta    = round(c(model_exp$cov.pars[1],  model_gaus$cov.pars[1], model_spe$cov.pars[1]), 4),
  Rango_phi = round(c(model_exp$cov.pars[2],  model_gaus$cov.pars[2], model_spe$cov.pars[2]), 6)
)

kable(tabla_modelos,
      caption   = "Tabla 3. Comparación de modelos teóricos de semivariograma según Suma de Cuadrados del Error (SCE)",
      col.names = c("Modelo","SCE","Pepita","Meseta parcial","Rango (phi)"),
      align     = c("l","r","r","r","r")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE) %>%
  row_spec(which.min(tabla_modelos$SCE), bold = TRUE, background = "#e8f5e9")

Tabla 3. Comparación de modelos teóricos de semivariograma según Suma de Cuadrados del Error (SCE)

Modelo	SCE	Pepita	Meseta parcial	Rango (phi)
Exponential	5678.80	0.0000	3.2017	0.000109
Gaussian	10284.05	0.0493	3.0762	0.000105
Spherical	8927.61	0.1062	3.0333	0.000258

plot(variograma,
     xlab = "Distancia (grados)", ylab = "Semivarianza",
     pch = 16, col = "grey30")
lines(model_exp,  col = paleta[1], lwd = 2)
lines(model_gaus, col = paleta[2], lwd = 2)
lines(model_spe,  col = paleta[3], lwd = 2)
legend("bottomright",
       legend = c("Exponential","Gaussian","Spherical"),
       col = paleta[1:3], lwd = 2, bty = "n")

Figura 10. Semivariograma empírico con los tres modelos teóricos ajustados: exponencial, gaussiano y esférico.

sce_vec    <- c(model_exp$value, model_gaus$value, model_spe$value)
nom_geoR   <- c("exponential","gaussian","spherical")
nom_disp   <- c("Exponential","Gaussian","Spherical")
mejor_idx  <- which.min(sce_vec)

mejor_modelo_nombre <- nom_geoR[mejor_idx]
mejor_modelo        <- list(model_exp, model_gaus, model_spe)[[mejor_idx]]

cat("Modelo seleccionado:", nom_disp[mejor_idx], "\n")

## Modelo seleccionado: Exponential

cat("SCE  :", round(mejor_modelo$value,       4), "\n")

## SCE  : 5678.797

cat("Nugget:", round(mejor_modelo$nugget,      4), "\n")

## Nugget: 0

cat("Sigma2:", round(mejor_modelo$cov.pars[1], 4), "\n")

## Sigma2: 3.2017

cat("Phi   :", round(mejor_modelo$cov.pars[2], 6), "\n")

## Phi   : 0.000109

# Variables inline del modelo seleccionado
res_modelo_nombre <- nom_disp[mejor_idx]
res_sce           <- round(mejor_modelo$value, 2)
res_nugget        <- round(mejor_modelo$nugget, 3)
res_sigma2        <- round(mejor_modelo$cov.pars[1], 3)
res_phi_grados    <- mejor_modelo$cov.pars[2]
res_phi_metros    <- round(mejor_modelo$cov.pars[2] * 111000, 0)

Interpretación: De los tres modelos ajustados, el Exponential obtuvo el menor SCE (5678.80), frente a 10284.05 del peor ajuste, lo que indica una diferencia sustancial en la capacidad de cada modelo para reproducir el semivariograma empírico. El rango efectivo estimado de 0.000109 grados equivale aproximadamente a 12 metros, que es la distancia máxima a la que dos árboles en esta finca presentan temperaturas correlacionadas. La meseta parcial de 3.202 °C² representa la varianza explicada por la estructura espacial, mientras que la pepita de 0.000 °C² refleja la variabilidad a escala menor que la separación mínima entre sensores o el error de medición propio del instrumento portátil. Una pepita elevada en relación a la meseta sería indicativa de que el sensor introduce ruido considerable.

6 Interpolación con Kriging

En esta sección se predice la temperatura sobre una grilla de 100 x 100 puntos mediante Kriging Ordinario.

lon_min <- min(geo_temp$coords[,"lon"])
lon_max <- max(geo_temp$coords[,"lon"])
lat_min <- min(geo_temp$coords[,"lat"])
lat_max <- max(geo_temp$coords[,"lat"])

cat("Longitud:", round(lon_min,6), "a", round(lon_max,6), "\n")

## Longitud: -76.7118 a -76.71022

cat("Latitud :", round(lat_min,6), "a", round(lat_max,6), "\n")

## Latitud : 2.392101 a 2.393634

# Columnas x e y requeridas por rasterFromXYZ
geodatos_grid <- expand.grid(
  x = seq(lon_min, lon_max, length.out = 100),
  y = seq(lat_min, lat_max, length.out = 100)
)

plot(geodatos_grid, pch = 1, cex = 0.2, col = "grey70",
     xlab = "Longitud", ylab = "Latitud")
points(geo_temp$coords, col = paleta[2], pch = 16, cex = 0.5)

Figura 11. Grilla de predicción de 100 x 100 puntos (gris) superpuesta sobre los puntos de muestreo (verde).

geodatos_ko <- krige.conv(
  geo_temp,
  loc   = geodatos_grid,
  krige = krige.control(
    nugget    = mejor_modelo$nugget,
    trend.d   = "cte",
    trend.l   = "cte",
    cov.pars  = mejor_modelo$cov.pars,
    cov.model = mejor_modelo_nombre
  )
)

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

cat("Rango predicciones:", round(range(geodatos_ko$predict),   3), "°C\n")

## Rango predicciones: 22.182 29.536 °C

cat("Rango varianza KO :", round(range(geodatos_ko$krige.var), 6), "\n")

## Rango varianza KO : 0.007075 3.307027

# Variables inline del Kriging
res_pred_min     <- round(min(geodatos_ko$predict), 1)
res_pred_max     <- round(max(geodatos_ko$predict), 1)
res_sd_error_max <- round(max(sqrt(geodatos_ko$krige.var)), 2)
res_sd_error_min <- round(min(sqrt(geodatos_ko$krige.var)), 3)
res_sd_umbral    <- round(max(sqrt(geodatos_ko$krige.var)) * 0.7, 2)

par(mar = c(4,4,3,2))
image(geodatos_ko,
      xlab = "Longitud", ylab = "Latitud",
      col  = heat.colors(100, rev = TRUE))
contour(geodatos_ko, add = TRUE, drawlabels = TRUE, col = "grey20", lwd = 0.8)
points(geo_temp$coords, pch = 16, cex = 0.4, col = "black")

Figura 12. Mapa de predicción espacial de la temperatura (°C) mediante Kriging Ordinario con isolíneas de contorno.

Interpretación: El mapa de predicción muestra que la temperatura estimada varía entre 22.2 °C y 29.5 °C dentro de los límites de la finca. Las isotermas revelan un gradiente espacial con el sector de mayores temperaturas ubicado en la zona de mayor concentración de árboles del Q4 identificados previamente. Desde el punto de vista agronómico, las zonas que superan los 27 °C de forma consistente merecen atención especial en cuanto a riego y sombrío, ya que el aguacate Hass presenta estrés calórico por encima de ese umbral durante períodos sostenidos. Es importante recordar que estas predicciones corresponden a un único día; para recomendaciones de manejo permanente se requeriría un análisis multitemporal.

# Graficar error estándar con ggplot2 (image.kriging no acepta val=)
df_error <- data.frame(
  x    = geodatos_grid$x,
  y    = geodatos_grid$y,
  sdKO = sqrt(geodatos_ko$krige.var)
)
df_pts <- data.frame(
  x = geo_temp$coords[,"lon"],
  y = geo_temp$coords[,"lat"]
)

ggplot(df_error, aes(x = x, y = y, fill = sdKO)) +
  geom_raster(interpolate = TRUE) +
  scale_fill_gradientn(colours = rev(heat.colors(100)), name = "SD error (°C)") +
  geom_point(data = df_pts, aes(x = x, y = y),
             inherit.aes = FALSE, pch = 16, size = 0.4, col = "black") +
  labs(x = "Longitud", y = "Latitud") +
  theme_minimal(base_size = 12) + coord_fixed()

Figura 13. Mapa de desviación estándar del error de predicción Kriging. Las zonas de mayor error coinciden con menor densidad de muestreo.

Interpretación: La desviación estándar del error de predicción oscila entre 0.084 °C y 1.82 °C. Los valores más bajos coinciden con el interior del lote donde la densidad de árboles monitoreados es mayor, mientras que los bordes y esquinas de la finca presentan los errores más altos, como es esperable en Kriging Ordinario cuando se extrapola más allá de la nube de datos. En términos prácticos, las predicciones en el borde occidental tienen una incertidumbre cerca del doble que las del centro del lote, lo que debería tenerse en cuenta antes de tomar decisiones de manejo diferencial basadas en esas zonas periféricas.

7 Predicción en RASTERS

En esta sección se transforma la predicción en rasters para visualización avanzada e integración con plataformas SIG.

require(rasterVis)

temp_predict <- rasterFromXYZ(
  cbind(geodatos_grid, z = geodatos_ko$predict)
)
temp_error_r <- rasterFromXYZ(
  cbind(geodatos_grid, z = sqrt(geodatos_ko$krige.var))
)

suppressWarnings({
  proj4string(temp_predict)  <- CRS("+proj=longlat +datum=WGS84")
  proj4string(temp_error_r)  <- CRS("+proj=longlat +datum=WGS84")
})

cat("Raster predicción creado.\n")

## Raster predicción creado.

print(temp_predict)

## class      : RasterLayer 
## dimensions : 100, 100, 10000  (nrow, ncol, ncell)
## resolution : 1.6e-05, 1.548485e-05  (x, y)
## extent     : -76.71181, -76.71021, 2.392093, 2.393642  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : z 
## values     : 22.18185, 29.53588  (min, max)

levelplot(temp_predict,
          par.settings = BuRdTheme,
          xlab = "Longitud", ylab = "Latitud")

Figura 14. Raster de temperatura predicha mediante Kriging Ordinario, visualizado con la paleta BuRd.

Interpretación: La representación raster, con resolución de celda determinada por la grilla de 100 x 100 puntos sobre la extensión de la finca, permite apreciar con mayor suavidad los gradientes térmicos que el mapa de isolíneas. La paleta BuRd, que va de azul (frío) a rojo (cálido), hace evidente que el núcleo más cálido de la finca no coincide con el centroide del lote sino que está desplazado hacia un sector específico. Este raster puede integrarse directamente en QGIS para cruzarlo con capas de pendiente, orientación o uso del suelo y así profundizar en las causas de la heterogeneidad térmica observada.

levelplot(temp_error_r,
          par.settings = BuRdTheme,
          xlab = "Longitud", ylab = "Latitud")

Figura 15. Raster de desviación estándar del error de predicción Kriging.

Interpretación: El raster de error estándar muestra que la incertidumbre no es homogénea en el lote: crece notoriamente en los bordes y en particular en la esquina con menor densidad de puntos. Con base en este mapa, si se quisiera reducir el error máximo a menos de 1.27 °C, sería suficiente con instalar sensores adicionales en las dos o tres zonas periféricas de mayor error, sin necesidad de aumentar la densidad de monitoreo en el centro del lote donde el modelo ya es preciso.

leaflet() %>%
  addTiles() %>%
  addRasterImage(temp_predict, opacity = 0.65) %>%
  addLegend(position = "bottomright",
            colors   = c("#2C7BB6","#FDAE61","#D7191C"),
            labels   = c("Baja","Media","Alta"),
            opacity  = 0.8)

Figura 16. Mapa interactivo con el raster de temperatura predicha superpuesto sobre la cartografía base.

Interpretación: Al superponer el raster de predicción sobre la cartografía base se puede confirmar que la finca se ubica en zona rural de Popayán, rodeada de otros cultivos y caminos que podrían influir en la dinámica de temperatura a través de variaciones en la cobertura del suelo y el flujo de aire. El sector más cálido predicho por el modelo coincide visualmente con la parte del lote con menor protección natural de colindantes, lo que es consistente con una mayor exposición a la radiación directa.

8 Validación Cruzada

En esta sección se evalúa la calidad predictiva del modelo mediante validación cruzada leave-one-out con xvalid de geoR.

valida <- xvalid(geodata = geo_temp, model = mejor_modelo)

## xvalid: number of data locations       = 534
## xvalid: number of validation locations = 534
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 
## xvalid: end of cross-validation

cat("Columnas xvalid:", paste(names(valida), collapse = ", "), "\n")

## Columnas xvalid: data, predicted, krige.var, error, std.error, prob

MAE  <- mean(abs(valida$error))
RMSE <- sqrt(mean(valida$error^2))
ME   <- mean(valida$error)
MSDR <- mean((valida$error^2) / valida$krige.var)

cat("MAE :", round(MAE,  4), "°C\n")

## MAE : 0.7861 °C

cat("RMSE:", round(RMSE, 4), "°C\n")

## RMSE: 1.0141 °C

cat("ME  :", round(ME,   4), "°C (sesgo)\n")

## ME  : -0.011 °C (sesgo)

cat("MSDR:", round(MSDR, 4), "(~1 = bien calibrado)\n")

## MSDR: 0.9304 (~1 = bien calibrado)

# Variables inline de la validación cruzada
res_MAE        <- round(MAE, 3)
res_RMSE       <- round(RMSE, 3)
res_ME         <- round(ME, 4)
res_MSDR       <- round(MSDR, 3)
res_MAE_redondo <- round(MAE, 1)
res_msdr_texto <- ifelse(abs(MSDR - 1) < 0.3,
  "se acerca razonablemente a 1, lo que indica que las varianzas Kriging son estimaciones realistas del error",
  "se aleja de 1, lo que sugiere que las varianzas Kriging subestiman o sobreestiman el error real y el modelo teórico podría afinarse")

tabla_val <- data.frame(
  Métrica    = c("MAE","RMSE","Error Medio (sesgo)","MSDR"),
  Valor      = round(c(MAE, RMSE, ME, MSDR), 4),
  Unidad     = c("°C","°C","°C","adimensional"),
  Referencia = c("Menor es mejor","Menor es mejor","Cercano a 0","Cercano a 1")
)

kable(tabla_val,
      caption = "Tabla 4. Métricas de desempeño del modelo Kriging en validación cruzada leave-one-out",
      align   = c("l","r","l","l")) %>%
  kable_styling(bootstrap_options = c("striped","hover"), full_width = FALSE)

Tabla 4. Métricas de desempeño del modelo Kriging en validación cruzada leave-one-out

Métrica	Valor	Unidad	Referencia
MAE	0.7861	°C	Menor es mejor
RMSE	1.0141	°C	Menor es mejor
Error Medio (sesgo)	-0.0110	°C	Cercano a 0
MSDR	0.9304	adimensional	Cercano a 1

errores_df <- data.frame(error = valida$error)

ggplot(errores_df, aes(x = error)) +
  geom_histogram(aes(y = after_stat(density)),
                 bins = 25, fill = paleta[1], color = "white", alpha = 0.85) +
  geom_density(color = paleta[2], linewidth = 1.2) +
  geom_vline(xintercept = 0, color = paleta[3], linetype = "dashed", linewidth = 1) +
  labs(x = "Error de predicción (°C)", y = "Densidad") +
  theme_minimal(base_size = 13)

Figura 17. Histograma de errores de predicción en la validación cruzada leave-one-out.

Interpretación: La validación cruzada leave-one-out sobre los 534 árboles arrojó un MAE de 0.786 °C y un RMSE de 1.014 °C. Estos valores indican que en promedio el modelo se equivoca en menos de 0.8 °C al predecir la temperatura de un árbol a partir de sus vecinos, lo que se considera aceptable dado el rango total observado de 7.5 °C. El error medio de -0.011 °C confirma que el modelo no tiene sesgo sistemático de sobreestimación ni subestimación. El MSDR de 0.93 se acerca razonablemente a 1, lo que indica que las varianzas Kriging son estimaciones realistas del error. Una limitación de esta validación es que al ser leave-one-out, no evalúa la capacidad del modelo para predecir en zonas sin ningún sensor cercano, que es precisamente donde el error del raster es mayor.

obs_pred_df <- data.frame(
  observado = valida$data,
  predicho  = valida$predict
)

ggplot(obs_pred_df, aes(x = observado, y = predicho)) +
  geom_point(color = paleta[1], alpha = 0.4, size = 1.8) +
  geom_abline(slope = 1, intercept = 0,
              color = paleta[2], linewidth = 1.2, linetype = "dashed") +
  geom_smooth(method = "lm", formula = y ~ x,
              color = paleta[3], se = TRUE, linewidth = 1,
              fill = paleta[5], alpha = 0.3) +
  labs(x = "Temperatura observada (°C)", y = "Temperatura predicha (°C)") +
  theme_minimal(base_size = 13) +
  theme(plot.title    = element_text(face = "bold", color = paleta[1]))

Figura 18. Diagrama de dispersión entre valores observados y predichos en la validación cruzada.

Interpretación: El diagrama observado-predicho muestra que la nube de puntos sigue razonablemente la diagonal de predicción perfecta, aunque con mayor dispersión en los extremos del rango de temperatura. Esto es típico del Kriging Ordinario, que tiende a suavizar los valores extremos porque los estima como promedios ponderados de vecinos: los árboles más calientes quedan ligeramente subestimados y los más fríos ligeramente sobreestimados. Este efecto de regresión a la media es una limitación inherente al método y debe considerarse si el objetivo fuera identificar con precisión los árboles individuales en condición de estrés térmico extremo.

9 Conclusiones

En esta sección se sintetizan los principales hallazgos del análisis.

El análisis confirmó que la temperatura en la finca no es espacialmente homogénea: el rango de 7.5 °C entre el árbol más frío y el más caliente, registrado en una sola jornada, es suficiente para generar condiciones microclimáticas diferenciadas que afectan el desarrollo del aguacate. El semivariograma empírico capturó esta estructura con un rango efectivo de aproximadamente 12 metros, que representa la escala espacial relevante para el manejo diferencial dentro del lote.

El modelo Exponential resultó el mejor ajuste con un SCE de 5678.8, y el Kriging Ordinario construido sobre él produjo predicciones con un MAE de 0.786 °C, error razonable para decisiones de manejo de precisión. El sector noroccidental de la finca concentra consistentemente las temperaturas más altas, lo que sugiere revisar las condiciones de sombrío, riego y densidad de siembra en esa zona.

Como limitaciones del estudio se destacan tres: el análisis cubre un único día de medición, lo que impide generalizar los patrones; el semivariograma se calculó sin analizar anisotropía, que podría existir dada la posible influencia de pendientes o corrientes de aire; y la validación cruzada leave-one-out no evalúa el desempeño en zonas sin sensores cercanos, donde el error real puede superar el 1.82 °C estimado en los bordes. Para fortalecer el análisis se recomienda incorporar mediciones de al menos tres períodos distintos y explorar semivariogramas direccionales.

---

sessionInfo()

## R version 4.6.0 (2026-04-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## 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    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] kableExtra_1.4.0 knitr_1.51       leaflet_2.2.3    ggplot2_4.0.3   
##  [5] rasterVis_0.51.7 lattice_0.22-9   raster_3.6-32    sp_2.2-1        
##  [9] geoR_1.9-6       dplyr_1.2.1      readxl_1.5.0    
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6        xfun_0.58           bslib_0.11.0       
##  [4] htmlwidgets_1.6.4   latticeExtra_0.6-31 vctrs_0.7.3        
##  [7] tools_4.6.0         crosstalk_1.2.2     generics_0.1.4     
## [10] parallel_4.6.0      proxy_0.4-29        tibble_3.3.1       
## [13] pkgconfig_2.0.3     KernSmooth_2.23-26  Matrix_1.7-5       
## [16] RColorBrewer_1.1-3  S7_0.2.2            lifecycle_1.0.5    
## [19] compiler_4.6.0      farver_2.1.2        stringr_1.6.0      
## [22] deldir_2.0-4        textshaping_1.0.5   terra_1.9-27       
## [25] codetools_0.2-20    class_7.3-23        htmltools_0.5.9    
## [28] sass_0.4.10         yaml_2.3.12         hexbin_1.28.5      
## [31] pillar_1.11.1       jquerylib_0.1.4     MASS_7.3-65        
## [34] classInt_0.4-11     cachem_1.1.0        nlme_3.1-169       
## [37] tidyselect_1.2.1    digest_0.6.39       stringi_1.8.7      
## [40] sf_1.1-1            labeling_0.4.3      splines_4.6.0      
## [43] fastmap_1.2.0       grid_4.6.0          cli_3.6.6          
## [46] magrittr_2.0.5      e1071_1.7-17        withr_3.0.2        
## [49] scales_1.4.0        rmarkdown_2.31      jpeg_0.1-11        
## [52] interp_1.1-6        otel_0.2.0          cellranger_1.1.0   
## [55] zoo_1.8-15          png_0.1-9           evaluate_1.0.5     
## [58] splancs_2.01-45     tcltk_4.6.0         viridisLite_0.4.3  
## [61] mgcv_1.9-4          rlang_1.2.0         Rcpp_1.1.1-1.1     
## [64] DBI_1.3.0           glue_1.8.1          xml2_1.5.2         
## [67] svglite_2.2.2       rstudioapi_0.18.0   jsonlite_2.0.0     
## [70] R6_2.6.1            units_1.0-1         systemfonts_1.3.2