Carga de librerías requeridas

require(geoR)
## Loading required package: geoR
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-2 (built on 2022-08-09) is now loaded
## --------------------------------------------------------------
require(MASS)
## Loading required package: MASS
require(readxl)
## Loading required package: readxl
require(rasterVis)
## Loading required package: rasterVis
## Loading required package: lattice
require(lattice)
require(raster)
## Loading required package: raster
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:MASS':
## 
##     select

Carga de Datos de árboles

geo_datos <- read_excel("./Datos_Completos_Aguacate.xlsx", sheet = "Hoja2")

Exploración de datos de árboles

head(geo_datos)
## # A tibble: 6 × 21
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Temp…¹
##      <dbl>    <dbl>     <dbl> <chr>                                        <dbl>
## 1        1     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22  
## 2        2     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.4
## 3        3     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.8
## 4        4     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.8
## 5        5     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.6
## 6        6     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.5
## # ℹ 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>

require(leaflet)
## Loading required package: leaflet
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.15,color = "blue")

Creación de la variable regionalizada

# Se selecciona la medición de humedad relativa y las coordenadas de latitud y longitud (2 y 3)
geo_hum=as.geodata(geo_datos,coords.col=2:3, data.col = 7)

plot(geo_hum)

### Matriz de distancias

D=dist(geo_hum$coords)
summary(D)
##      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(geo_hum,option = "bin",uvec=seq(0,9.178e-04,0.000018))
## variog: computing omnidirectional variogram
plot(variograma)

variograma_mc=variog.mc.env(geo_hum,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)

Ajuste del Modelo de Semivarianza

# Se realiza para el eje Y y eje X el seq.
ini.vals = expand.grid(seq(5,25,l=100), seq(0,2e-04,l=100))
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 "25"    "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1611649.48721746
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 "25"    "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1702390.47809032
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 "25"    "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 1705567.87228394
plot(variograma)
lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="purple")

Al graficar el comportamiento de cada modelo, es posible observar que el modelo exponencial es el modelo que mejor se ajusta. Adicionalmente, este tiene un menor valor de pérdida respecto a los modelos gaussiano y esférico.

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

Predicción

min(geo_datos[,3])
## [1] -76.7118
geodatos_grid=expand.grid( lon=seq(-76.7118,-76.71022,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")

Para realizar la interpolación usando Kriging, utilizaremos el Phi =26.3256 y sigmasq=0.0001 del modelo que mejor se ajusta a los puntos del semivariograma

geodatos_ko=krige.conv(geo_hum, loc=geodatos_grid,
      krige= krige.control(nugget=0,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=26.3256, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
image(geodatos_ko, main="Interpolacion con Kriging", 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")

#contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var), add=TRUE, drawlabels=TRUE)

Raster de la predicción

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