1. Importación de base de datos

Se importa la base de datos y las librerias que se utilizarán en el ejercicio. A su vez, se filtra para que solo se utilice la fecha de 2020-10-01, la cual es la que compete al análisis. También se hace algunos tratamientos a la variable temperatura, como transformarla a numérica.

library(sp)
library(gstat)
library(sf)
library(raster)
library(readxl)
library(dplyr)
library(ggplot2)
library(viridis)


datos <- read_excel("Datos_Completos_Aguacate.xlsx")


# Filtrar el primer período (01/10/2020)
datos$Fecha <- as.Date(sub(" .*", "", datos$FORMATTED_DATE_TIME), format = "%d/%m/%Y")

datos <- subset(datos, Fecha == as.Date("2020-10-01"))


datos <- datos %>%
  mutate(Temperature = gsub(",", ".", Temperature)) %>%
  mutate(Temperature = as.numeric(Temperature))



datos_sf <- st_as_sf(datos, coords = c("Longitude", "Latitude"), crs = 4326)

datos_sf <- datos_sf[, c("id_arbol", "Temperature")]

2. Análisis exploratorio

Se realiza el análisis exploratorio de datos en busca de patrones que puedan llegar a contribuir a una mejor comprensión de la variable temperatura. Para esto, se calculan algunas estadísticas de tendencia central y un histrograma para observar su distribución. Se observa que la temperatura sigue una distribución normal, una media de 25,83° centigrados y oscila entre 22° a 29° grados centigrados.

summary(datos$Temperature)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70
hist(datos$Temperature, main="Distribución de la Temperatura", xlab="Temperature", breaks=30)

3. Variograma

Se realiza el variograma para analizar la autocorrelación espacial de la variable temperatura y, a su vez, encontrar el mejor modelo para hacer predicciones.

datos_sp <- as(datos_sf, "Spatial")

g <- gstat(id = "Temperature", formula = Temperature ~ 1, data = datos_sp)

# Cálculo del variograma
vgm_exp <- variogram(g)

plot(vgm_exp)

Dada la forma de su distribución, se busca el modelo que mejor se ajusta. Para esto, se despliegan una serie de opciones con la funcion vgm() y show.vgms. Como se puede observar, el modelo esférico es el que más se parece a la forma mostrada por el variograma, por lo que se emplea este para ajustar el modelo.

vgm()
##    short                                      long
## 1    Nug                              Nug (nugget)
## 2    Exp                         Exp (exponential)
## 3    Sph                           Sph (spherical)
## 4    Gau                            Gau (gaussian)
## 5    Exc        Exclass (Exponential class/stable)
## 6    Mat                              Mat (Matern)
## 7    Ste Mat (Matern, M. Stein's parameterization)
## 8    Cir                            Cir (circular)
## 9    Lin                              Lin (linear)
## 10   Bes                              Bes (bessel)
## 11   Pen                      Pen (pentaspherical)
## 12   Per                            Per (periodic)
## 13   Wav                                Wav (wave)
## 14   Hol                                Hol (hole)
## 15   Log                         Log (logarithmic)
## 16   Pow                               Pow (power)
## 17   Spl                              Spl (spline)
## 18   Leg                            Leg (Legendre)
## 19   Err                   Err (Measurement error)
## 20   Int                           Int (Intercept)
show.vgms(par.strip.text = list(cex = 0.75))

# Ajustar un modelo teórico esférico
vgm_model <- fit.variogram(vgm_exp, model = vgm("Sph"))

# Grafca del semivariograma
plot(vgm_exp, model = vgm_model)

5. Predicción kriging

Se realiza la predicción de la temperatura mediante el método kirging a los lugares donde no hay información capturada. Para esto, se emplea el modelo anteriormente calculado y se proyecta las predicciones sobre un espacio de malla regular. Para facilitar la visualización de los datos, una vez realizada la predicción se transforma la data en un shape file y se grafica tanto la predicción como la varianza, esta última con el propósito de entender los lugares donde la información existente no es suficiente para realizar predicciones precisas y, por lo tanto, una mayor incertidumbre en la predicción.

# Creación de malla
grd <- st_make_grid(datos_sf, n = 100, what = "centers")
grd_sf <- st_as_sf(grd, crs = 4326)
grd_sp <- as(grd_sf, "Spatial")

# Predicción
k <- gstat(formula = Temperature ~ 1, data = datos_sp, model = vgm_model)
kpred <- predict(k, grd_sp)
## [using ordinary kriging]
# Convertirsión a shape file
kpred_sf <- st_as_sf(kpred)

# Gráfica predicciones
ggplot() + 
  geom_sf(data = kpred_sf, aes(color = var1.pred)) +
  geom_sf(data = datos_sf) +
  scale_color_viridis(name = "Temperature") + 
  theme_bw()

# Gráfica varianza
ggplot() + 
  geom_sf(data = kpred_sf, aes(color = var1.var)) +
  geom_sf(data = datos_sf) +
  scale_color_viridis(name = "Variance") + 
  theme_bw()


Como se observa, cada punto representa un arbol y los colores la temperatura en cada punto tras la predicción. Se obtienen algunas zonas interesantes con temperaturas bien establecidas. En la zona este y sur oeste se observan las temperaturas más altas alrededor de 28° centigrados (zona amarilla), La zona norte es una parte fria, con temperaturas alrededor de los 24° centigrados. Por último, la zona más fria corresponde a la zona que atraviesa de forma diagonal el sur, centro y este del mapa con temperaturas de 23° centigrados.

A su vez, la varianza de la predicción empieza a aumentar a medida que se va alejando de los arboles, es decir, de los registros que se tenian en la base de datos con respecto a la temperatura. Sin embargo, alcanza a abarcar una zona de alto interes para el análisis.


6. Conclusiones