Mediante el uso de la geoestadística, se requiere realizar un análisis detallado sobre el clima y mejores condiciones geográficas para la simebra y cultivo del mismo.

Se procede a cargar la información de la base de datos correspondiente a los datos completos del aguate:

Aguacates <- read_excel("C:/Users/User/Documents/Unidad3/Datos_Completos_Aguacate.xlsx")
head(Aguacates)
## # A tibble: 6 × 21
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME       Psychro_Wet_Bulb_Tempe…¹
##   <chr>       <dbl>     <dbl> <chr>                                        <dbl>
## 1 1            2.38     -76.6 21/08/2019  9:22:57 a, m,                     14.8
## 2 2            2.38     -76.6 21/08/2019  9:27:13 a, m,                     11.6
## 3 3            2.38     -76.6 21/08/2019  9:36:36 a, m,                     12.9
## 4 4            2.38     -76.6 21/08/2019  9:38:02 a, m,                     14.1
## 5 5            2.38     -76.6 21/08/2019  9:39:38 a, m,                     14.3
## 6 6            2.38     -76.6 21/08/2019  9:42:02 a, m,                     14.2
## # ℹ abbreviated name: ¹​Psychro_Wet_Bulb_Temperature
## # ℹ 16 more variables: Station_Pressure <dbl>, Relative_Humidity <dbl>,
## #   Crosswind <dbl>, Temperature <dbl>, Barometric_Pressure <dbl>,
## #   Headwind <dbl>, Direction_True <dbl>, Direction_Mag <dbl>,
## #   Wind_Speed <dbl>, Heat_Stress_Index <dbl>, Altitude <dbl>, Dew_Point <dbl>,
## #   Density_Altitude <dbl>, Wind_Chill <dbl>,
## #   Estado_Fenologico_Predominante <dbl>, Frutos_Afectados <dbl>

Posterior al cargue de la información de la base de datos es conveniente graficar en mapa la plantación de los aguacates:

require(leaflet)

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

Objeto Geodata

Se convierten los datos en una variable regionalizada.

require(geoR)
geod_aguacate=as.geodata(Aguacates,coords.col = 3:2,data.col = 9)
## as.geodata: 18586 replicated data locations found. 
##  Consider using jitterDupCoords() for jittering replicated locations. 
## WARNING: there are data at coincident or very closed locations, some of the geoR's functions may not work.
##  Use function dup.coords() to locate duplicated coordinates.
##  Consider using jitterDupCoords() for jittering replicated locations
plot(geod_aguacate)

Al encontrarse la situación de la gráfica geodata anterior, se procede a realizar el análisis con la variable temperatura, y la fecha 21 de agosto de 2019.

Aguacates$fecha <- as.Date(Aguacates$FORMATTED_DATE_TIME, format = "%d/%m/%Y")
Data_20 <- filter(Aguacates, fecha == "2020-10-01")

tablaAguacate <- (head(Data_20,5))

tablaAguacate %>%
  kbl() %>%
  kable_paper("hover", 
              full_width = F)
id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Temperature Station_Pressure Relative_Humidity Crosswind Temperature Barometric_Pressure Headwind Direction_True Direction_Mag Wind_Speed Heat_Stress_Index Altitude Dew_Point Density_Altitude Wind_Chill Estado_Fenologico_Predominante Frutos_Afectados fecha
1 2.393549 -76.71124 01/10/2020  10:11:12 a, m, 22.0 825.1 85.2 0.0 23.9 825.2 0.0 313 312 0.0 25.3 1696 21.3 2.504 23.9 717 0 2020-10-01
2 2.393573 -76.71120 01/10/2020  10:11:12 a, m, 21.4 825.3 84.0 0.0 23.5 825.2 0.0 317 317 0.0 24.8 1696 20.7 2.485 23.5 717 0 2020-10-01
3 2.393541 -76.71113 01/10/2020  10:11:12 a, m, 21.8 825.5 79.6 0.2 24.5 825.5 0.4 338 337 0.5 25.7 1694 20.8 2.518 24.5 717 0 2020-10-01
4 2.393503 -76.71119 01/10/2020  10:11:12 a, m, 22.8 825.4 77.6 0.4 25.9 825.4 0.2 299 299 0.5 28.1 1694 21.7 2.572 25.9 717 0 2020-10-01
5 2.393486 -76.71121 01/10/2020  10:11:12 a, m, 22.6 825.2 76.5 0.0 26.0 825.2 0.0 265 264 0.0 28.0 1696 21.5 2.575 25.9 717 0 2020-10-01 
geod_aguacate=as.geodata(Data_20,coords.col = 3:2,data.col = 9)
plot(geod_aguacate)

### Semivariograma

library(leaflet)
leaflet() %>% addTiles() %>% addCircleMarkers(lng = Data_20$Longitude,lat = Data_20$Latitude,radius = 0.2,color = "blue")
#summary(dist(Data_20[,3:2]))
summary(dist(Data_20[,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
#variograma=variog(geod_aguacate,option = "bin",uvec=seq(0,0.001,0.00002))

variograma=variog(geod_aguacate,option = "bin",uvec=seq(0,0.000917,9.178e-05))
## variog: computing omnidirectional variogram
datos.env=variog.mc.env(geod_aguacate,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

Graficando el variograma se obtiene lo siguiente:

plot(variograma)
lines(datos.env)

AJuste del modelo

Se realiza el ajuste del modelo de semivarianza de la siguiente forma:

Usando el modelo exponencial:

ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,0.0008,l=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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 240018.669934895

Usando el modelo esférico:

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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 226841.340327613

Usando el modelo Gaussiano:

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 "1.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 231440.514329651

Graficando todos los modelos se observa:

# Gráficando los módelos.
plot(variograma)

lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="green")

Mejorando la gráfica:

# Suponiendo que variograma, model_mco_exp, model_mco_gaus, y model_mco_spe están definidos

# Crear el gráfico del variograma
plot(variograma, main="Variograma y Modelos Ajustados", 
     xlab="Distancia", ylab="Semivarianza", 
     pch=19, col="black", cex=1.5)

# Agregar las líneas de los modelos ajustados
lines(model_mco_exp, col="blue", lwd=2, lty=1)
lines(model_mco_gaus, col="red", lwd=2, lty=2)
lines(model_mco_spe, col="green", lwd=2, lty=3)

# Agregar una leyenda
legend("topright", legend=c("Modelo Exponencial", "Modelo Gaussiano", "Modelo Esférico"), 
       col=c("blue", "red", "green"), lwd=2, lty=1:3, cex=0.8, bty="n")

# Mejorar la apariencia general del gráfico
par(mar=c(5, 5, 4, 2) + 0.1) # Ajustar márgenes

Los resultados de los modelos precedentes son:

model_mco_exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7849  2.2672  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0003390645
## 
## variofit: minimised weighted sum of squares = 3907.716

El exponencial es el que presenta menor suma de cuadrados.

model_mco_gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7594  2.2508  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
## 
## variofit: minimised weighted sum of squares = 5241.893
model_mco_spe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  2.9538  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 12203.46

Predicción espacial

c(min(Aguacates[,3]),
max(Aguacates[,3]),
min(Aguacates[,2]),
max(Aguacates[,2]))
## [1] -76.711799 -76.606710   2.316405   2.393634

Se continua con la creación de la malla:

#geodatos_grid=expand.grid( lon=seq(-76.711799,-76.606710,l=100),lat=seq(2.316405 ,2.393634 ,l=100))
#plot(geodatos_grid)
#points(geod_aguacate$coords,col="red")


datos_grid=expand.grid( lon=seq(-76.710,-76.712,l=100),lat=seq(2.3920 ,2.3937 ,l=100))
plot(datos_grid)
points(Data_20[,3:2],col="red")

Predicción del modelo

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

par(mfrow=c(1,2))
contour(geodatos_ko,main="kriging Predict", drawlabels=TRUE)

image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")