Nota: Documento académico para la Maestría en ciencia de datos de la Universidad Javeriana de Cali

Introducción

El propósito de este análisis es estudiar la variabilidad espacial de variables regionalizadas de un cultivo, utilizando datos recolectados de sensores instalados en árboles de aguacate. Los datos serán filtrados para seleccionar un primer periodo de medición, específicamente el 1 de octubre de 2020. Este análisis se enfoca en la variable Temperature.

Para este estudio, se utilizará el video de YouTube titulado “[Video] M2U1- Clima en finca de aguacate”:

Ver video de YouTube aquí

En resumen, según lo comprendido en el vídeo, los pasos para realizar el análisis son los siguientes:

  1. Filtrar los datos climáticos: Se seleccionará el primer periodo de medición, específicamente el 1 de octubre de 2020, para realizar un análisis de la variable Temperature.
  2. Análisis de autocorrelación espacial: Utilizar un semivariograma para identificar la existencia de autocorrelación espacial en los datos de temperatura.
  3. Identificación del mejor modelo teórico: Seleccionar el modelo teórico de semivariograma que mejor se ajuste a los datos entre esférico, exponencial y gaussiano.
  4. Predicción espacial mediante kriging: Utilizar el modelo teórico seleccionado para realizar una interpolación espacial de la temperatura en el área de estudio mediante el método de kriging, generando un mapa de predicción espacial.

1 cargue de datos

# Leer el archivo de datos
datos_aguacate <- read_excel("Datos_Completos_Aguacate.xlsx")

# Filtrar los datos por la fecha específica
datos_aguacate <- datos_aguacate %>%
  filter(as.Date(FORMATTED_DATE_TIME, format="%d/%m/%Y") == as.Date("2020-10-01"))

# Ver el número de registros después de filtrar
num_registros <- nrow(datos_aguacate)
print(paste("Número de registros para el 01/10/2020:", num_registros))
## [1] "Número de registros para el 01/10/2020: 534"
leaflet() %>%
  addTiles() %>%
  addCircleMarkers(
    lng = datos_aguacate$Longitude, 
    lat = datos_aguacate$Latitude,
    radius = 0.2,
    color = "black"
  )

2 patrones de temperatura

# Cargar la librería geoR
library(geoR)

# Convertir los datos en un objeto geodata
geo_temp <- as.geodata(datos_aguacate, coords.col = c("Longitude", "Latitude"), data.col = "Temperature")

# Graficar el objeto geodata
plot(geo_temp)

Las zonas más cálidas están concentradas en el sur-occidente y en la franja oriental del área de estudio. Las áreas más frías se encuentran en el centro y el noroeste. Esto sugiere un patrón de variabilidad espacial de la temperatura, posiblemente influido por otros factores geográficos locales. Los datos muestran una clara autocorrelación espacial.

Se estudiarán registros de 534 árboles como menciona el vídeo.

# Obtener estadísticas descriptivas de la variable Temperature
# Opción 1: Usando summary()
summary(datos_aguacate$Temperature)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70

La temperatura varía entre 22.2 °C y 29.7 °C, con una media de 25.83 °C y una mediana de 25.8 °C, lo que indica una distribución bastante de 10°C entre los sensores de loárboles de aguacate.

3 Semivariograma, análisis de autocorrelación espacial

Limité el semivariograma al tercer cuartil de las distancias para enfocarme en el rango donde es más probable encontrar autocorrelación espacial, evitando el ruido de las distancias largas. Esto me permite obtener un semivariograma más claro y confiable.

# Calcular las distancias entre los puntos de muestreo (coordenadas de longitud y latitud)
distancias <- dist(datos_aguacate[, c("Longitude", "Latitude")])

# Resumen de las distancias
summary(distancias)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.712e-05 4.051e-04 6.408e-04 6.827e-04 9.178e-04 1.959e-03
# Calcular el semivariograma omnidireccional limitado al percentil 75 de las distancias
variograma <- variog(geo_temp, option = "bin", uvec = seq(0, 0.0009178, 9.178e-05))
## variog: computing omnidirectional variogram
# Graficar el semivariograma
plot(variograma)

# Calcular las envolventes mediante simulaciones
variograma_mc <- variog.mc.env(geo_temp, obj = variograma)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
# Añadir las envolventes al gráfico
lines(variograma_mc)

En este gráfico, estamos analizando cómo varía la temperatura entre puntos a medida que aumenta la distancia entre ellos. Vemos que, en distancias cortas, la diferencia de temperatura es baja (puntos cercanos tienen temperaturas similares), lo que muestra que hay una relación fuerte entre la temperatura y la ubicación. A medida que la distancia aumenta y se acerca al límite del tercer cuartil (nuestro rango de análisis), esta relación se pierde y las diferencias de temperatura se estabilizan, indicando que los puntos alejados tienen temperaturas más independientes. Las líneas superior e inferior muestran el rango de variación que esperaríamos si la temperatura no estuviera influenciada por la ubicación. En distancias cortas, varios puntos están fuera de este rango, lo que confirma que hay autocorrelación espacial: puntos cercanos tienen temperaturas similares. Sin embargo, en distancias más largas, los puntos caen dentro de las líneas, indicando que la influencia de la ubicación disminuye y las temperaturas se vuelven más independientes unas de otras.

En este semivariograma, la mayoría de la autocorrelación espacial se concentra en el rango de distancias de 0 a 4e-04, donde los puntos de semivarianza están por debajo de las líneas de confianza, indicando dependencia espacial significativa. A partir de 4e-04, los puntos se acercan a las envolventes, lo que sugiere que la autocorrelación disminuye y los valores de la variable se vuelven más independientes.

Este comportamiento indica que, a partir de estas distancias, los valores de temperatura entre puntos se vuelven independientes. Por lo tanto, voy a ajustar el modelo de semivarianza en el rango de 2e-04 a 4e-04-

4 Identificación del mejor modelo teórico

# Ajuste de modelos de semivariograma en el rango (2e-04 a 5e-04)
ini.vals = expand.grid(seq(8, 12, l=10), seq(2e-04, 4e-04, l=10))

# Ajuste de cada modelo
plot(variograma)
model_mco_exp = variofit(variograma, ini=ini.vals, cov.model="exponential", wei="npair", min="optim")
## 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 845476.826874663
## Warning in variofit(variograma, ini = ini.vals, cov.model = "exponential", :
## unreasonable initial value for sigmasq (too high)
model_mco_gaus = variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim", nugget=0)
## 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1283514.7904775
## Warning in variofit(variograma, ini = ini.vals, cov.model = "gaussian", :
## unreasonable initial value for sigmasq (too high)
model_mco_spe = variofit(variograma, ini=ini.vals, cov.model="spheric", fix.nug=TRUE, wei="npair", min="optim")
## 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 2307931.2112438
## Warning in variofit(variograma, ini = ini.vals, cov.model = "spheric", fix.nug
## = TRUE, : unreasonable initial value for sigmasq (too high)
lines(model_mco_exp, col="blue")
lines(model_mco_gaus, col="red")
lines(model_mco_spe, col="purple")

# Valores del semivariograma experimental
observed_semivariances <- variograma$v
observed_distances <- variograma$u

# Funciones de los modelos de semivarianza
exponential_model <- function(h, cov.pars) {
  sigmasq <- cov.pars[1]
  phi <- cov.pars[2]
  sigmasq * (1 - exp(-h / phi))
}

gaussian_model <- function(h, cov.pars) {
  sigmasq <- cov.pars[1]
  phi <- cov.pars[2]
  sigmasq * (1 - exp(-(h / phi)^2))
}

spherical_model <- function(h, cov.pars) {
  sigmasq <- cov.pars[1]
  phi <- cov.pars[2]
  ifelse(h <= phi, sigmasq * (1.5 * (h / phi) - 0.5 * (h / phi)^3), sigmasq)
}

# Calcular los valores predichos por cada modelo
predicted_exp <- exponential_model(observed_distances, model_mco_exp$cov.pars)
predicted_gaus <- gaussian_model(observed_distances, model_mco_gaus$cov.pars)
predicted_spe <- spherical_model(observed_distances, model_mco_spe$cov.pars)

# Calcular ECM para cada modelo
ecm_exp <- mean((observed_semivariances - predicted_exp)^2)
ecm_gaus <- mean((observed_semivariances - predicted_gaus)^2)
ecm_spe <- mean((observed_semivariances - predicted_spe)^2)

# Imprimir los ECM para comparar
print(paste("ECM Exponencial:", ecm_exp))
## [1] "ECM Exponencial: 0.78386198832747"
print(paste("ECM Gaussiano:", ecm_gaus))
## [1] "ECM Gaussiano: 6.09305529875679"
print(paste("ECM Esférico:", ecm_spe))
## [1] "ECM Esférico: 1.03288272083215"
lines(model_mco_exp, col="blue")
lines(model_mco_gaus, col="red")
lines(model_mco_spe, col="purple")

En el gráfico, observo los ajustes de tres modelos de semivariograma: exponencial (línea roja), gaussiano (línea azul) y esférico (línea púrpura). El modelo exponencial tiene el menor Error Cuadrático Medio (ECM) de 0.78 y un valor de pérdida bajo (845476.83), lo que indica un buen ajuste general y representa bien la estructura de autocorrelación espacial. El modelo esférico también se ajusta razonablemente bien, con un ECM de 1.03, pero su valor de pérdida es significativamente mayor.

El modelo gaussiano, sin embargo, tiene problemas. No se ajusta en las distancias cortas, pero representa bien las mayores distancias donde la autocorrelación disminuye (es decir, en las zonas donde los puntos se acercan a las líneas de confianza), lo que se refleja en su ECM más alto (6.09). Por lo tanto, el modelo exponencial es la mejor elección aquí para capturar la estructura de autocorrelación en los datos.

5 Predicción

# Verificar los límites de latitud y longitud en los datos
min(datos_aguacate$Longitude)
## [1] -76.7118
max(datos_aguacate$Longitude)
## [1] -76.71022
min(datos_aguacate$Latitude)
## [1] 2.392101
max(datos_aguacate$Latitude)
## [1] 2.393634
# Crear la cuadrícula de predicción basada en los límites de los datos
geodatos_grid <- expand.grid(
  lon = seq(-76.7118, -76.71022, length.out = 100),
  lat = seq(2.392101, 2.393634, length.out = 100)
)

# Visualizar la cuadrícula y los puntos de observación
plot(geodatos_grid, main = "Cuadricula de Prediccion")
points(datos_aguacate[, c("Longitude", "Latitude")], col = "red")

# Definir los parámetros obtenidos del modelo exponencial
sigmasq <- model_mco_exp$cov.pars[1]  # valor de sigmasq del modelo exponencial
phi <- model_mco_exp$cov.pars[2]      # valor de phi del modelo exponencial

# Realizar la predicción de kriging usando el modelo exponencial
geodatos_ko <- krige.conv(
  geo_temp,
  loc = geodatos_grid,
  krige = krige.control(
    nugget = 0,
    trend.d = "cte",
    trend.l = "cte",
    cov.pars = c(sigmasq = sigmasq, phi = phi)
  )
)
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
# Graficar la predicción de kriging y su desviación estándar
par(mfrow = c(1, 2))

# Mapa de predicción de kriging
image(geodatos_ko, main = "Kriging Prediction", xlab = "Longitude", ylab = "Latitude")
contour(geodatos_ko, add = TRUE, drawlabels = TRUE)

# Mapa de desviación estándar de la predicción de kriging
image(geodatos_ko, main = "Kriging StDv Predicted", val = sqrt(geodatos_ko$krige.var), xlab = "Longitude", ylab = "Latitude")
contour(geodatos_ko, val = sqrt(geodatos_ko$krige.var), add = TRUE, drawlabels = TRUE)

Los gráficos muestran que las zonas más calientes están concentradas en el suroccidente, a lo largo de una franja en el oriente, y en algunas áreas al norte. El modelo de kriging logró capturar bien estos patrones de temperatura. En cuanto a la incertidumbre (gráfico de desviación estándar), se observa que en el centro del área la predicción es más confiable, es decir de donde hay datos, ya que la desviación estándar es baja, mientras que en los bordes la incertidumbre aumenta (donde no hay datos), indicando menos precisión en esas zonas.

# Cargar la librería necesaria
require(rasterVis)
## Cargando paquete requerido: rasterVis
## Cargando paquete requerido: lattice
# Combinar la cuadrícula de predicción con los valores predichos de kriging
pred <- cbind(geodatos_grid, geodatos_ko$predict)

# Crear el raster de predicción a partir de los valores de kriging
temp_predict <- rasterFromXYZ(pred)

# Graficar el raster de predicción
plot(temp_predict, main = "Kriging Prediction")

Finalmente, convertí los valores predichos de kriging en un raster para visualizar la temperatura en la región con una escala continua de colores. Esto nos permite observar fácilmente las variaciones de temperatura en el área,

6 Conclusión

El análisis de estos mapas de predicción de Kriging en un cultivo de aguacates nos permite identificar zonas de variación térmica dentro del área estudiada. Las áreas en verde en el mapa de predicción de temperaturas son las más cálidas, mientras que las áreas en tonos rosados son más frías. Esto indica que existen patrones de temperatura que podrían influir en el crecimiento del cultivo. Las zonas más cálidas están hacia el suroeste y en una franja en el oriente, lo que sugiere que estas áreas pueden tener condiciones diferentes que deben considerarse en la gestión del cultivo.

Fin del documento

dibujar