Taller 3-Exploración de datos con Geoestadística
Presentado por Juan Sebastián Henao
A continuación, se presentan los resultados de un análisis de exploración de datos con geoestadística. Para esto, vamos a realizar una predicción de temperatura alrededor de un set de datos de una finca de aguacate. Para esto, lo primero que vamos a hacer es cargar las distintas librerías, luego haremos una exploración de los datos, por último se calculan 3 modelos para la predicción y luego se realiza la predicción a partir del modelo de que mejor se ajusta a partir del MSE.
library(readxl)
library(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
## --------------------------------------------------------------library(raster)## Loading required package: splibrary(ggplot2)
library(leaflet)
require(rasterVis)## Loading required package: rasterVis## Loading required package: latticegeo_datos =read_excel("~/Desktop/datos_aguacate.xlsx")
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>Exploración
Para iniciar, vamos a realizar una exploración de los datos que se contienen en el set de datos. Para esto, el primer paso es visualizar con leaflet los puntos que tiene la base. Esto muestra que los datos están agrupados y hacen referencia a una misma unidad geográfica, que en este caso, es una finca. Sobre este conjunto, es que se realizará el ejercicio.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = geo_datos$Longitude,lat = geo_datos$Latitude,radius = 0.2,color = "black")Luego de visualizar los puntos, se procede a transformar la base en un archivo como geodata y luego visualizar las características que este tiene.
geo_ag=as.geodata(geo_datos, coords.col = 2:3, data.col = 9)
plot(geo_ag)
Luego, calculamos las distancias en los datos para poder hacer el cálculo del semivariograma.
m=dist(geo_ag$coords)
summary(m)## 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-03Con los datos que nos arroja el summary de m, podemos realizar el cálculo del semivariograma tomando como base los datos del tercer cuartil y el mínimo.
variograma=variog(geodata = geo_ag, uvec= seq(0,9.9e-04,1.99e-05))## variog: computing omnidirectional variogramplot(variograma)
Con esto en mente, realizamos el plot del variograma mostrando que no se encuentra entre las líneas, lo cual indica que si existe correlación. Luego de esta exploración, procedemos a realizar el cálculo de los modelos y selección para la predicción.
datos_env = variog.mc.env(geo_ag, obj.variog =variograma)## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelopsplot(variograma, pch=16, envelop = datos_env)
Cálculo de los modelos y su selección
Para este caso, con los datos que nos arroja el modelo, vamos a calcular un modelo exponencial, gaussiano y esférico, y a partir de los cuales se hará la selección.
ini.vals = expand.grid(seq(1, 3, l = 10), seq(6e-04, 9e-04, l = 10))
plot(variograma, pch=16, envelop = datos_env)
modelo_mco_exp = variofit(variograma, ini.vals, cov.model = "exponential") ## 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 "3" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 204754.631532749modelo_mco_gaus = variofit(variograma, ini.vals, cov.model = "gaussian") ## 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 "3" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 261629.832919343modelo_mco_sph = variofit(variograma, ini.vals, cov.model = "spheric") ## 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 "3" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 35901.582017216lines(modelo_mco_exp, col = "blue")
lines(modelo_mco_gaus, col = "red")
lines(modelo_mco_sph, col = "green")
legend(0.0007, 6, legend=c("Gaussiano", "Exponencial", "Esférico"),
col=c("red", "blue","green"), lty=1, cex=0.8)
Lo que nos muestra la gráfica, es que aparentemente el modelo que mejor se ajusta a los datos es el exponencial. Pero antes de esto, vamos a ver el resultado de MSE para los tres modelos y confirmar esta selección.
MSE_exp <- modelo_mco_exp$value
MSE_gaus <- modelo_mco_gaus$value
MSE_sph <- modelo_mco_sph$value
MSE_vector <- c(Exponencial = MSE_exp, Gaussiano = MSE_gaus, Esférico = MSE_sph)
print(MSE_vector)## Exponencial Gaussiano Esférico
## 3865.548 40674.146 15245.417mejor_modelo <- names(MSE_vector)[which.min(MSE_vector)]
cat("El mejor modelo es:", mejor_modelo, "con MSE =", min(MSE_vector))## El mejor modelo es: Exponencial con MSE = 3865.548Luego de comparar el MSE de los tres modelos, se confirma que el exponencial es el que mejor se ajusta con un valor de 3865.55.
Predicción con el modelo seleccionado
Teniendo en cuenta que el modelo que mejor se ajusta a los datos es el exponencial, se inicia imprimiendo los parámetros para el modelo ajustado.
parametros = variofit(variograma, ini.vals, cov.model = "exponential")## 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 "3" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 204754.631532749print(parametros)## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 2.1786 4.1706 0.0024
## Practical Range with cor=0.05 for asymptotic range: 0.007257417
##
## variofit: minimised weighted sum of squares = 3865.548Ahora, para poder hacer la predicción, obtengamos las coordenadas del plano sobre las cuales se hará la predicción
# Obtener los mínimos y máximos de latitud y longitud
min_latitud <- min(geo_ag$coords[, 1]) # La primera columna es latitud
max_latitud <- max(geo_ag$coords[, 1])
min_longitud <- min(geo_ag$coords[, 2]) # La segunda columna es longitud
max_longitud <- max(geo_ag$coords[, 2])
# Imprimir los valores mínimos y máximos
cat("Mínimo de latitud:", min_latitud, "\n")## Mínimo de latitud: 2.392101cat("Máximo de latitud:", max_latitud, "\n")## Máximo de latitud: 2.393634cat("Mínimo de longitud:", min_longitud, "\n")## Mínimo de longitud: -76.7118cat("Máximo de longitud:", max_longitud, "\n")## Máximo de longitud: -76.71022Con estos datos, se imprime el plano sobre el cual se realizará y se toman los parámetros del modelo ajustado.
datos_grid = expand.grid(Latitude = seq(2.392,2.394, l=100),Longitud = seq(-76.7118,-76.71022, l=100))
datos_ko = krige.conv(geo_ag, loc = datos_grid, krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", cov.pars = c(sigmasq = 4.1706, phi = 0.0024)))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodplot(datos_grid)
Por último, se presenta la predicción de temperatura en el plano. Esto lo realizamos con Kriging y donde se muestra que hay consistencia con los datos de la exploración.
image(datos_ko, main = "Kriging Predict - Temperatura", xlab ="Latitude", ylab ="Longitude")
Ahora, trasformamos los datos en raster y realizamos una visualización con levelplot, que permite una visualización más amena y permite observar las variaciones.
pred=cbind(datos_grid,datos_ko$predict)
temp_predict=rasterFromXYZ(cbind(datos_grid,datos_ko$predict))
levelplot(temp_predict,par.settings =BuRdTheme)
En conclusión, lo que muestra esta predicción, es que las áreas resaltadas en color rojo representan las zonas donde se observan las temperaturas más altas durante el análisis exploratorio.