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)

Semivariograma

summary(dist(geod_aguacate$coords))
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000000 0.001756 0.077784 0.063476 0.106252 0.109962
variograma=variog(geod_aguacate,option = "bin",uvec=seq(0,0.001,0.00002))
## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origin
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: 5351625789.89722

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: 5297606635.41097

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: 5325691774.74162

Graficando todos los modelos se observa:

plot(variograma)

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

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 
## 12.5756  0.9163  0.0013 
## Practical Range with cor=0.05 for asymptotic range: 0.004020816
## 
## variofit: minimised weighted sum of squares = 1234962

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 
## 12.8480  0.0011  0.0013 
## Practical Range with cor=0.05 for asymptotic range: 0.002165812
## 
## variofit: minimised weighted sum of squares = 1484166
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 
## 12.8448  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 23884586

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.606710,-76.711799,l=100),lat=seq(2.316405 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geod_aguacate$coords,col="red")

Predicción del modelo

# Supongamos que tus datos están en un objeto llamado geod_aguacate
# y que la columna coords contiene las coordenadas

# Determinar los rangos de las coordenadas x e y
x.range <- range(geod_aguacate$coords[, 1])
y.range <- range(geod_aguacate$coords[, 2])

# Ajustar la resolución del grid
# Un valor más alto de grid_resolution resultará en un grid menos denso
grid_resolution <- 1  # Ajusta este valor según tus necesidades

# Crear el grid con la resolución ajustada
geodatos_grid <- expand.grid(
  x = seq(x.range[1], x.range[2], by = grid_resolution),
  y = seq(y.range[1], y.range[2], by = grid_resolution)
)

# Verificar el tamaño del grid
cat("Tamaño del grid:", nrow(geodatos_grid), "puntos\n")
## Tamaño del grid: 1 puntos
# Ver un resumen del grid
head(geodatos_grid)
##          x        y
## 1 -76.7118 2.316405