El caso de estudio esta centrado en realizar la predicción de la variable temperatura en un cultivo de aguacate en el Cauca usando el interpolador geoestadístico Kriging. Para esto, se definen los siguientes pasos:

1. Análisis exploratorio de los datos

Se inserta la base de datos que contiene la ubicación de cada árbol de aguacate en el departamento del Cauca.

library(readxl)
datos = read_excel("C:/Maria Camila/Maestría Ciencia de datos/2do semestre/Análisis de información geográfica/Modulo 3/Datos_Completos_Aguacate.xlsx")
head(datos)

En la tabla se puede observar las variables climáticas de los registros de cada árbol como son, la humedad relativa, temperatura, dirección del viento, entre otras.

Para el análisis de predicción de la variable temperatura usando el interpolador Kriging, se hace la exploración de datos sujetos a su ubicación espacial

library(leaflet)
leaflet() %>% addTiles() %>% addCircleMarkers(lng = datos$Longitude,lat = datos$Latitude,radius = 0.2,color = "orange")

De acuerdo con la visualización anterior, existen 4 lotes donde se encuentran los cultivos de aguacate cercanos al municipio de Timbio, Cauca

Seguidamente, se realiza el filtro para 01/10/2020

library(dplyr)
datos$fecha <- as.Date(datos$FORMATTED_DATE_TIME, format = "%d/%m/%Y")
# Filtrar los datos por la fecha específica
datos_2020 <- datos %>%
  filter(FORMATTED_DATE_TIME == "01/10/2020  10:11:12 a, m,")
# Muestra los datos filtrados
head(datos_2020)

Ahora visualizamos estos datos

library(leaflet)
leaflet() %>% addTiles() %>% addCircleMarkers(lng = datos_2020$Longitude,lat = datos_2020$Latitude,radius = 0.2,color = "blue")

Este conjunto de árboles de aguacate se encuentra cercano al sector conocido del Cauca “La Chorrera”

En consecuencia, se procede a visualizar los datos filtrados en relación a su ubicación y distribución de temperaturas correspondientes a cada árbol de aguacate en octubre del 2020.

library(geoR)
library(raster)
temp=as.geodata(datos_2020,coords.col = 3:2,data.col = 9)
plot(temp)

En el primer gráfico, se relacionan los niveles de temperaturas, siendo las zonas de color azul con temperaturas más bajas, las verdes y amarillas temperaturas medias y las rojas son áreas donde se encuentran las temperaturas más altas.

2. Semivariograma Experimental

Con la finalidad de analizar la variabilidad espacial de la variable temperatura en relación con su ubicación y distancia, se procede a ejecutar el semivariograma experimental o empírico; para ello, es relavante conocer los datos mínimos y máximos, así, como los cuartiles (25%, 50% y 75%). Dado que estos valores serán los utilizados para definir el sill y el range del semivariograma.

summary(dist(datos_2020[,3:2]))
##      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
library(geoR)
variograma=variog(temp,option = "bin",uvec=seq(0,0.000917,9.178e-05))
## variog: computing omnidirectional variogram
plot(variograma)
variograma_temp=variog.mc.env(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
lines(variograma_temp)

En el semivariograma experimental, se observa que los puntos negros corresponden a la variabilidad de cada ubicación de acuerdo a la distancias, que en este caso se observa en el semivariograma una autocorrelación espacial fuerte, lo que indica que en los puntos de temperatura para árbol que estén cercanos unos a otros presentan temperaturas similares en contraste con las distantes.

3. Ajuste del semivariograma experimental

Ahora, se hace el análisis de cuál es el mejor modelo matemático en relacióna la ubicación de los puntos frente a la variabilidad que tengan de la variable temperatura. Se evalua el semivariograma experimental con el modelo Esférico, Gaussiano y Exponencial.

ini.vals = expand.grid(seq(1,15,l=10), seq(6e-04,8e-05,l=10))

plot(variograma)
model_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 "4.11"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 20717.4624927274
model_gaus=variofit(variograma, ini=ini.vals, cov.model="gaussian", wei="npair", min="optim")
## 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 "2.56"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 25714.590255448
model_sph=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 "2.56"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 25766.449209192
lines(model_exp,col="blue")
lines(model_gaus,col="red")
lines(model_sph,col="purple")

model_gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.2189  2.7757  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002328613
## 
## variofit: minimised weighted sum of squares = 6537.706
model_exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.0042  3.0606  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0003115809
## 
## variofit: minimised weighted sum of squares = 3531.275
model_sph
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  3.0049  0.0002 
## Practical Range with cor=0.05 for asymptotic range: 0.0001955546
## 
## variofit: minimised weighted sum of squares = 5817.208

De acuerdo, al desempeño de los modelos teóricos con el comportamiento del semivariograma empírico, se define que el mejor modelo teórico que representa la variabilidad de la temperatura es el Exponencial de acuerdo con la suma de mínimos cuadrados 3531.275, en relación al modelo gaussiano y esférico.

4. Predicción espacial

Ahora vamos a realizar la interpolación (predicción espacial), esta se hace por medio de la definción de una cuadrícula o grilla de coordenadas que se define desde el punto mínimo al máximo de longitud y latitud, que representa toda el área de estudio de la variable temperatura en el cultivo de aguacate.

##Predicción Espacial Kriging
datos_grid=expand.grid( lon=seq(-76.710,-76.712,l=100),lat=seq(2.3920 ,2.3937 ,l=100))
plot(datos_grid)
points(datos_2020[,3:2],col="red")

5. Procesamiento con el interpolador kriging para la predicción de la variable temperatura

En el procesamiento para la ejecución de la predicción de la variable temperatura en el cultivo de aguacate en el Cauca, se procede a la determinación de los datos por medio de los componentes de su semivariograma teórico, que para este caso fue el exponencial, como el nugget, sigmasq y phi.

datos_ko=krige.conv(temp, loc=datos_grid,
      krige= krige.control(nugget=0,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=3.0606, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(datos_ko, main="kriging Predict", xlab="Longitud", ylab="Latitud")

En el gráfico anterior, se observa la predicción de Kriging en cuánto a la temperatura y su distribución espacial en el cultivo de aguacate.

image(datos_ko, main="kriging StDv Predicted",val=sqrt(datos_ko$krige.var), xlab="Longitud", ylab="Latitud")

En relación al gráfico de la desviación estándar en la prediccón se comporta de forma similar en cada uno de los puntos, lo que no afecta la predicción de la temperatura.

Finalmente, se transforma la imagen de la predicción con Kriging a raster, con la finalidad de tener una análisis más descriptivo.

library(raster)
library(rasterVis)
pred=cbind(datos_grid,datos_ko$predict)
temp_predict=rasterFromXYZ(cbind(datos_grid,datos_ko$predict))
colores <- colorRampPalette(c("blue", "green", "yellow", "red"))
plot(temp_predict, col.regions= colores, main = "Predicción de Temperatura con kriging")

De acuerdo al gráfico anterior, podemos evidenciar que en la zona suroeste hay una predicción de que la temperatura se comporte entre los 28°C a 29°C, en contraste, en la zona central del cultivo la probabilidad es que la temperatura tenga un comportamiento mucho menor entre los 23°C y 25°C, esto puede deberse a la topografía y geomorfología del lote.