Clima en finca de Aguacate

Teniendo en cuenta la base de datos de los registros de árboles de aguacates en el Cauca, para el primer periodo de la base (01-10-2020) que corresponde a 534 registros, realizar las siguientes tareas.

  1. Análisis exploratorio 1.1. Crear el geodata 1.2. Analizar si existe tendencia en la variable
  2. Evaluar la autocorrelación espacial con el semivariograma
  3. Realizar ajuste con el mejor modelo teórico (exponencial, gaussiano, esférico)
  4. Realizar el ajuste espacial con la metodología Kriging

1. Análisis exploratorio

Dado que el ejercicio requiere filtrar la base por una fecha específica, para realizar este proceso inicialmente se extrae la fecha de la columna FORMATTED_DATE_TIME y posteriormente se realiza el filtro solicitado

# Extraer fecha de la columna FORMATTED_DATE_TIME
data_completa <- data_completa %>%
  mutate(
    fecha = str_extract(FORMATTED_DATE_TIME, "^\\d{2}/\\d{2}/\\d{4}")
  )

# Filtrar los registros solicitados para el análisis
periodo1 <- data_completa %>% 
  filter(as.Date(fecha) == as.Date("01/10/2020"))

# Seleccionar las columnas que se utilizarán del dataframe
periodo1 <- periodo1 %>% 
  dplyr::select(Latitude, Longitude, Temperature)

# graficar la distrubución de los puntos de coordenadas
plot(periodo1[,2:3])

La distribución de los puntos de coordenadas indica que el muestreo se realizó por zonas específicas dado que los puntos se aglomeran en ciertos rangos de coordenadas y no están dispersos a lo largo de todo el plano.

1.1. Crear el geodata

Se crea el geodata para realizar el análisis exploratorio de datos. Dado que las coordenadas están en grados decimales, es necesario convertir estas coordenadas a una proyección métrica y presentarlas como un dataframe para que la función as.geodata pueda tomar los datos correctamente y posteriormente poder calcular el semivariograma.

Conversión de coordenadas de grados decimales a UTM

# Convertir a objeto sf
p1_sf <- st_as_sf(periodo1, coords = c("Longitude", "Latitude"), crs = 4326)

# Convertir sf a UTM - Zona 18S - Se toma este código como el que más usado en la zona 18S a la que pertenece Cauca
p1_utm <- st_transform(p1_sf, 32618)

Convertir en dataframe

df_geo <- data.frame(
  x = st_coordinates(p1_utm)[,1],
  y = st_coordinates(p1_utm)[,2],
  temp = periodo1$Temperature
)

Crear el geodata

p1_geod <- as.geodata(
  df_geo,
  coords.col = 1:2,
  data.col = 3
)

# graficar geodata
plot(p1_geod)

1.2. Analizar si existe tendencia en la variable

Los resultados gráficos del geodata indican que la temperatura no tiene diferencias en una dirección específica dentro del espacio geográfico, lo que indicaría que es estacionaria (Gráfico superior izquierdo).

Por otro lado los gráficos que representan la temperatura frente a la latitud y longitud no indican una clara tendencia ascendente o descendente en ninguna dirección, lo cual confirmaría que la temperatura no se comporta de manera específica en ciertas ubicaciones geográficas (Gráficos superior derecho e inferior izquierdo). Por último el gráfico inferior derecho no indica sesgos marcados en la distribución de las temperaturas, tampoco datos atípicos en esta variable.

2. Evaluar la autocorrelación espacial con el semivariograma

En esta etapa se calcula la matriz de distancias para evaluar los sitios más cercanos y lejanos de los puntos donde se tomó la información de temperatura y con ello definir el rango que se utilizará en el semivariograma

# Cálculo de la matriz de distancias
m <-dist(p1_geod$coords)
summary(m)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.893  44.935  71.060  75.688 101.743 217.351

A partir del summary se puede identificar que las distancias entre los puntos son en general oscilan entre 1.893 y 217.351. Para el semivariograma se graficará hasta el 3er cuartil.

# Semivariograma
semiv <- variog(
  geodata = p1_geod,
  uvec = seq(0, 101, length = 20)
)
## variog: computing omnidirectional variogram
plot(semiv)

Si se analiza el semivariograma, se encuentra que la temperatura tiene patrones espaciales claros, donde la semivarianza es baja al inicio (< 20mts) lo cual indica que la temperatura no presenta mayor cambio en distancias cortas, luego hay mayores y más rápidos aumentos de temperatura a medida que se incrementa la distancia (entre 20 y 30 mts). A partir de esa distancia, no hay un cambio significativo de la temperatura lo que indica que a distancias grandes no existen patrones espaciales y, gráficamente, se confirmaría la autocorrelación de esta variable.

Luego se realizan simulaciones de Montecarlo para analizar si al aleatorizar esta variable el comportamiento es similar y para indicar los rangos de valores donde se mueven los datos cuando no hay correlación espacial

# Se realizan 500 simulaciones de Montecarlo para establecer las bandas
bandas <- variog.mc.env(geodata = p1_geod, obj.variog = semiv, nsim = 500)
## variog.env: generating 500 simulations by permutating data values
## variog.env: computing the empirical variogram for the 500 simulations
## variog.env: computing the envelops
# Graficar variograma y bandas
plot(semiv)
lines(bandas, col="purple")

3. Realizar ajuste con el mejor modelo teórico (exponencial, gaussiano, esférico)

Se ajustan los datos a diferentes modelos teóricos para identificar cuál tiene mejor ajuste, utilizando la función variofit

# Modelo esférico
mod_esf <- variofit(
  semiv,
  cov.model = "spherical",
  wei = "equal"
)
## variofit: covariance model used is spherical 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
## Warning in variofit(semiv, cov.model = "spherical", wei = "equal"): initial
## values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi     tausq  kappa
## initial.value "2.66"  "32.32" "0.36" "0.5"
## status        "est"   "est"   "est"  "fix"
## loss value: 1.52959128345751
# Modelo exponencial
mod_exp <- variofit(
  semiv,
  cov.model = "exponential",
  wei = "equal"
)
## variofit: covariance model used is exponential 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
## Warning in variofit(semiv, cov.model = "exponential", wei = "equal"): initial
## values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi     tausq  kappa
## initial.value "2.66"  "32.32" "0.89" "0.5"
## status        "est"   "est"   "est"  "fix"
## loss value: 0.93433119931597
# Modelo gaussiano
mod_gauss <- variofit(
  semiv,
  cov.model = "gaussian",
  wei = "equal"
)
## variofit: covariance model used is gaussian 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
## Warning in variofit(semiv, cov.model = "gaussian", wei = "equal"): initial
## values not provided - running the default search
## variofit: searching for best initial value ... selected values:
##               sigmasq phi     tausq  kappa
## initial.value "2.66"  "16.16" "0.36" "0.5"
## status        "est"   "est"   "est"  "fix"
## loss value: 2.15612311976169
#Gráficos de los 3 modelos
plot(semiv)
lines(mod_esf, col = "blue", lwd = 2)
lines(mod_exp, col = "purple", lwd = 2)
lines(mod_gauss, col = "orange", lwd = 2)

legend(
  "bottomright",
  legend = c("Esférico", "Exponencial", "Gaussiano"),
  col = c("blue", "purple", "orange"),
  lwd = 2
)

Teniendo en cuenta los diferentes modelos teóricos, se encuentra que el que más se ajusta al semivariograma es el modelo exponencial.

Por otro lado, utilizando el paquete gstat y las funciones variogram, vgm y fit.variogram se presentan los siguientes resultados.

# Realizar variograma con función variogram
semiv2 <- variogram(Temperature ~ 1, data = p1_utm)

#Definición de modelos teóricos iniciales con vgm

# Se definen los valores para el variograma iniciales teniendo en cuenta el rango del variograma de los datos

# Modelo esférico
vt_esf <- vgm(psill = 3.2, model = "Sph",
                range = 60, nugget = 0.2)

# Modelo exponencial
vt_exp <- vgm(psill = 3.2, model = "Exp",
                range = 60, nugget = 0.2)

# Modelo Gaussiano
vt_gauss <- vgm(psill = 3.2, model = "Gau",
                range = 60, nugget = 0.2)

# Ajustar variogramas

# Exponencial
fit_exp <- fit.variogram(semiv2, vt_esf)

# Esférico
fit_esf <- fit.variogram(semiv2, vt_exp)

# Gaussiano
fit_gauss <- fit.variogram(semiv2, vt_gauss)

A continuación se grafican las funciones ajustadas a los modelos

# Gráfica del semivariograma y los modelos teóricos 

#Modelo esférico
plot(semiv2, fit_esf)

#Modelo exponencial
plot(semiv2, fit_exp)

#Modelo Gaussiano
plot(semiv2, fit_gauss)

De manera similar al gráfico que presentaba el variograma con los 3 modelos teóricos, se valida que bajo este método, igualmente el modelo que más se ajusta es el exponencial.

4. Realizar el ajuste espacial con la metodología Kriging

Para realizar la predicción espacial utilizando la metodología Kriging y generar la imagen de la predicción espacial, se utiliza la función gstat. La metodología Kriging utiliza los valores del variograma ajustado (a alguno de los modelos teóricos) para obtener predicciones en ubicaciones no muestreadas. En este caso, teniendo en cuenta que el modelo exponencial es el que se observa gráficamente con mejor ajuste respecto al modelo esférico y Gaussiano, se aplica esta metodología con el ajuste de este modelo.

# Construir grilla basada en los límites UTM (p1_utm está en 32618)
bb <- st_bbox(p1_utm)

x.range <- seq(bb["xmin"], bb["xmax"], length.out = 200)
y.range <- seq(bb["ymin"], bb["ymax"], length.out = 200)

grid <- expand.grid(x = x.range, y = y.range)

# Convertir grilla a sf (Asegurando el CRS de cálculo: 32618)
grid_sf <- st_as_sf(grid, coords = c("x","y"), crs = 32618)

# Construir el modelo de Kriging 
k <- gstat(formula = log(Temperature) ~ 1, data = p1_utm, model = fit_exp)

# Obtener predicciones (Objeto sf de PUNTOS en UTM)
kpred_puntos <- predict(k, grid_sf)
## [using ordinary kriging]
# Transformar los resultados a Lat/Lon (EPSG 4326), para que ggplot pueda manejar correctamente los ejes geográficos.
kpred_latlon <- st_transform(kpred_puntos, 32618)

# Transformar datos nuevamente a valores de temperatura (exp) y Extracción de Coordenadas (Data Frame)
coords <- st_coordinates(kpred_latlon)

df_kriging <- data.frame(
  lon = coords[,1],
  lat = coords[,2],
  
  #Conversión de Logaritmo a Real
  pred_real = exp(kpred_latlon$var1.pred),
  var_real = exp(kpred_latlon$var1.var)
)

# Calcular Límites de los rangos de predicciones para graficar adecuadamente
min_pred <- min(df_kriging$pred_real, na.rm = TRUE)
max_pred <- max(df_kriging$pred_real, na.rm = TRUE)
min_var <- min(df_kriging$var_real, na.rm = TRUE)
max_var <- max(df_kriging$var_real, na.rm = TRUE)

# Visualización

# Gráfico 1: Predicción de Temperatura
g1 <- ggplot(df_kriging, aes(x = lon, y = lat)) +
  # geom_tile es la solución robusta para rasters
  geom_tile(aes(fill = pred_real)) +
  # coord_sf define el sistema de coordenadas. Usamos 4326 que coincide con la data
  coord_sf(crs = 32618) + 
  scale_fill_viridis(option = "C", name = "Temp (°C)",
                     limits = c(min_pred, max_pred)) +
  theme_bw(base_size = 12) +
  labs(title = "Kriging: Predicción de Temperatura (Lat/Lon)",
       x = "Longitud", y = "Latitud")

# Gráfico 2: Varianza
g2 <- ggplot(df_kriging, aes(x = lon, y = lat)) +
  geom_tile(aes(fill = var_real)) +
  coord_sf(crs = 32618) +
  scale_fill_viridis(option = "A", name = "Varianza",
                     limits = c(min_var, max_var)) + 
  theme_bw(base_size = 12) +
  labs(title = "Kriging: Varianza de Predicción (Lat/Lon)",
       x = "Longitud", y = "Latitud")

# Mostrar gráficos
g1

g2

De acuerdo con los resultados de la predicción, que se observan como una superficie suavizada de los datos, se encuentra que en el noroeste y suroeste se encuentran las temperaturas más cálidas (colores amarillo y naranja), lo cual es coherente con las zonas identificadas en el gráfico descriptivo inicial que mostraba la distribución de los puntos de temperatura de los árboles en colores de acuerdo con los niveles de temperatura. Así mismo, coinciden las zonas con las temperaturas más bajas en el rango de datos (color morado, azul oscuro).

De otra parte, el gráfico de varianza indica que las predicciones son confiables, toda vez que la mayor parte del área de predicción se encuentra en color púrpura o morado, donde se encuentran las zonas de baja varianza y por, ende, se indica que el error en las predicciones es muy bajo, por lo que se puede confiar en la temperatura estimada para las áreas predichas.