ACTIVIDAD 3 - Análisis de Geoestadistica

Introducción

El presente trabajo tiene como objetivo principal analizar y modelar la distribución espacial de la variable temperatura en una zona del departamento del Cauca, Colombia, utilizando técnicas de análisis estadístico y geoestadística. La base de datos comprende un total de 534 registros correspondientes a mediciones realizadas en árboles de producción de aguacate. Cada árbol cuenta con un identificador único y las coordenadas geográficas de su ubicación, lo que permite realizar un análisis espacial detallado.

Además, la información inicial incluye variables climáticas recolectadas mediante sensores portátiles, tales como humedad relativa, velocidad del viento, temperatura y altura, entre otras. Sin embargo, este estudio se enfocará exclusivamente en la variable de temperatura, dejando de lado las demás variables climáticas para centrar los recursos analíticos en su modelado y predicción.

La metodología propuesta incluye:

Construcción del semivariograma experimental: Se calcularán las diferencias espaciales de la temperatura para representar su variación y dependencia espacial.

Ajuste del mejor modelo teórico: Se evaluarán distintos modelos teóricos, como el exponencial, gaussiano o esférico, para seleccionar el que mejor describa el comportamiento espacial de los datos.

Interpolación espacial: A partir del modelo ajustado, se aplicará la técnica de krigeado para predecir la temperatura en ubicaciones no medidas y generar una representación gráfica de la distribución espacial de la variable.

Este análisis permitirá identificar patrones espaciales de la temperatura en el área de estudio, lo que resulta clave para la toma de decisiones en la gestión del cultivo de aguacate, optimizando recursos y entendiendo mejor las condiciones climáticas locales que afectan la producción.

Exploración de los datos

datos=read.xlsx("D:/Docs C/Desktop/Sergio/Maestria/Segundo Semestre/Georeferenciación/Actividad 3/Datos_Completos_Aguacate.xlsx")

# Filtrar los datos para tomar solo las filas con la fecha específica
datos_filtrados <- subset(datos, FORMATTED_DATE_TIME == "01/10/2020  10:11:12 a, m,")

# Verificar los datos filtrados
head(datos_filtrados)
##      id_arbol Latitude Longitude        FORMATTED_DATE_TIME
## 9140        1 2.393549 -76.71124 01/10/2020  10:11:12 a, m,
## 9141        2 2.393573 -76.71120 01/10/2020  10:11:12 a, m,
## 9142        3 2.393541 -76.71113 01/10/2020  10:11:12 a, m,
## 9143        4 2.393503 -76.71119 01/10/2020  10:11:12 a, m,
## 9144        5 2.393486 -76.71121 01/10/2020  10:11:12 a, m,
## 9145        6 2.393441 -76.71126 01/10/2020  10:11:12 a, m,
##      Psychro_Wet_Bulb_Temperature Station_Pressure Relative_Humidity Crosswind
## 9140                         22.0            825.1              85.2       0.0
## 9141                         21.4            825.3              84.0       0.0
## 9142                         21.8            825.5              79.6       0.2
## 9143                         22.8            825.4              77.6       0.4
## 9144                         22.6            825.2              76.5       0.0
## 9145                         21.5            825.0              77.7       0.3
##      Temperature Barometric_Pressure Headwind Direction_True Direction_Mag
## 9140        23.9               825.2      0.0            313           312
## 9141        23.5               825.2      0.0            317           317
## 9142        24.5               825.5      0.4            338           337
## 9143        25.9               825.4      0.2            299           299
## 9144        26.0               825.2      0.0            265           264
## 9145        24.5               825.0      0.0            265           264
##      Wind_Speed Heat_Stress_Index Altitude Dew_Point Density_Altitude
## 9140        0.0              25.3     1696      21.3            2.504
## 9141        0.0              24.8     1696      20.7            2.485
## 9142        0.5              25.7     1694      20.8            2.518
## 9143        0.5              28.1     1694      21.7            2.572
## 9144        0.0              28.0     1696      21.5            2.575
## 9145        0.3              25.5     1698      20.4            2.522
##      Wind_Chill Estado_Fenologico_Predominante Frutos_Afectados
## 9140       23.9                            717                0
## 9141       23.5                            717                0
## 9142       24.5                            717                0
## 9143       25.9                            717                0
## 9144       25.9                            717                0
## 9145       24.5                            717                0
total_filas_filtradas <- nrow(datos_filtrados)
cat("El DataFrame filtrado tiene", total_filas_filtradas, "filas.\n")
## El DataFrame filtrado tiene 534 filas.

geo_datos=datos_filtrados
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")

Creación Variable Regionalizada

geodatos=as.geodata(datos_filtrados,coords.col=3:2,datos_filtrados.col=9)
plot(geodatos)

summary(dist(geodatos$coords))
##      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
variograma=variog(geodatos,option = "bin",uvec=seq(0,9.178e-04,9.178e-05))
## variog: computing omnidirectional variogram
plot(variograma)
variograma_mc=variog.mc.env(geodatos,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_mc)

Se odserva que existe correlación espacial, al observar la temperatura se puede ver que a distacias cortas la similitud es baja, los valores de temparatura son parecidos y a mayores distancias estos valores cambian.

Ajuste al modelo de semivarianza

ini.vals = expand.grid(seq(5e-07,2e-07,l=10), seq(8e-04,10e-04,l=10))
plot(variograma)

options(digits = 10)

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 "0"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 8.73794905221097e-11
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 "0"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1.74142138294591e-12
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 "0"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 3.99417562052865e-10
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")

# Acceder a los parámetros ajustados
sigmasq <- model_mco_gaus$cov.pars[1]  # Este es el valor de sigmasq
phi <- model_mco_gaus$cov.pars[2]     # Este es el valor de phi

# Ver los valores ajustados
print(paste("Sigmasq:", sigmasq))
## [1] "Sigmasq: 3e-07"
print(paste("Phi:", phi))
## [1] "Phi: 0.001"

Como se observa en la figura, el modelo gaussiano es el que mejor se ajusta a los puntos, lo anterior se puede comprobar con el valor del calculo de la suma de cuadrados de los errores (SCE), este modelo es que presenta un valor mas bajo indicando que el modelo predice bien los valores observados.

Predicción espacial Krigind

geodatos_grid=expand.grid( lon=seq(-76.7120,-76.7100,l=100),lat=seq(2.3920 ,2.3940 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")

Se usa los datos del modelo gaussanio Phi =26.3256 y sigmasq=0.0001

geodatos_ko=krige.conv(geodatos, loc=geodatos_grid,
      krige= krige.control(nugget=0,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=3e-07, phi=0.001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")
contour(geodatos_ko,main="kriging Predict", add=TRUE, drawlabels=TRUE)

image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")

pred=cbind(geodatos_grid,geodatos_ko$predict)
hum_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
plot(hum_predict)