1. Carga de datos y visualización inicial

Se incorpora una base de datos que contiene registros climáticos asociados a árboles de aguacate, recopilados en el departamento del Cauca. Este gráfico indica la ubicación espacial de los arboles.

geo_datos =read_excel("Datos_Completos_Aguacate.xlsx")
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "green")

2. Geoestadística

El propósito del análisis estadístico es determinar si existe correlación espacial en las mediciones de temperatura, utilizando técnicas geoestadísticas en R. Para ello, se comienza creando el objeto geodata, que es esencial para llevar a cabo el análisis

2.1. Crear geodata

El objeto se construye a partir de la base de datos previamente cargada, especificando las columnas correspondientes a la longitud y latitud de las mediciones, así como la columna que contiene la variable de interés, que en este caso es la temperatura (columna 5).

geo_temp=as.geodata(geo_datos,coords.col = 3:2,data.col = 5)
plot(geo_temp)

En este análisis preliminar se observa que las temperaturas registradas oscilan entre 19 y 25°C, con ciertas zonas donde tienden a concentrarse predominantemente valores altos o bajos.

2.2. Variograma

Inicialmente, se realiza un análisis de las distancias entre los puntos con el fin de definir el rango que se utilizará en la construcción del variograma. En este caso, las distancias consideradas varían entre 0 y 0.001, con incrementos de 0.00002.

options(scipen=999)
summary(dist(geo_datos[,3:2]))
##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## 0.00001712 0.00040512 0.00064078 0.00068267 0.00091776 0.00195913
variograma=variog(geo_temp,option = "bin",uvec=seq(0,0.001,0.00002))
## variog: computing omnidirectional variogram
variograma_mc=variog.mc.env(geo_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
plot(variograma)
lines(variograma_mc)

El semivariograma sugiere la presencia de correlación espacial, dado que a menores distancias se observa una semivarianza más baja. Específicamente, la semivarianza comienza alrededor de 0.7 para las distancias más cortas y se estabiliza en aproximadamente 1.3 a una distancia cercana a 0.00025.

2.3. Ajuste del modelo teórico

Se procedió a ajustar tres modelos teóricos de semivariograma: el exponencial, el gaussiano y el esférico.

ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0003,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: 6624.4190268492
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.47"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 11320.824524992
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.33"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 1216.01931350849
plot(variograma)

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

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

De acuerdo con el gráfico del semivariograma y el criterio de mínima suma de cuadrados, el modelo que mejor se ajusta a los datos es el exponencial. A partir de este ajuste, se estimaron los parámetros: sigmasq =1.3138 y phi=0.0001.

2.4. Interpolación (Predicción espacial)

Para realizar la interpolación, se construye una malla que abarca la extensión espacial de los puntos de medición. Esta malla se utiliza como insumo en el proceso de kriging, junto con los datos observados y los parámetros estimados del modelo exponencial.

c(min(geo_datos[,3]),
max(geo_datos[,3]),
min(geo_datos[,2]),
max(geo_datos[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634
geodatos_grid=expand.grid( Este=seq(-76.7120,-76.7100,l=100),Norte=seq(2.3918 ,2.3940 ,l=100))

plot(geodatos_grid)
points(geo_datos[,3:2],col="red")

geodatos_ko=krige.conv(geo_temp, loc=geodatos_grid,
                       krige= krige.control(nugget=0,trend.d="cte", 
                                            trend.l="cte",cov.pars=c(sigmasq=1.3138, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

Predicción espacial:

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

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

Desviación estandar:

par(mfrow=c(1,2))

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

2.5. Transformar imagen a raster

Finalmente se transforma el resultado anterior en un objeto raster. Para esto primero se realiza un corte donde es de interes revisar la predicción.

pts_sf <- st_as_sf(geo_datos, coords =3:2, crs = 4326)
convex_hull <- st_convex_hull(st_union(pts_sf))
convex_sp <- as(convex_hull, "Spatial")
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))

r_crop <- mask(temp_predict, convex_sp)

levelplot(r_crop,par.settings =BuRdTheme)

Este resultado permite obtener estimaciones de temperatura en ubicaciones donde no se realizaron mediciones directas.