Condiciones Climáticas en finca de Aguacate
Carlos Valbuena Acosta
2023-09-27
Análisis de Información Geográfica y Espacial
Maestría en Ciencia de Datos - Javeriana Cali
Caso de Estudio.
Se toma como referencia para este análisis las mediciones obtenidas para unos árboles de Aguacate en el departamento del Cauca.
Para este caso, se utilizará la información capturada en el año 2020.
Desarrollo del Caso
Importación de Librerías
Se realiza la importación de las librerías para el análisis geoestadístico como es el caso de geoR.
# Importar librerías
require(readxl)
require(geoR)
require(leaflet)
require(lea)## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'lea'require(gstat)Cargue de la Base de Datos
Se realiza el cargue de la base de datos, donde se observa que el set de datos contiene 21 columnas con características incluidas para la georreferenciación que contienen los registros de 534 árboles.
# Se realiza el cargue de la base de datos de la Finca de Aguacate
BD_Aguacate_2020 <- read_excel("C:/Users/User/OneDrive - PUJ Cali/Maestría en Ciencia de Datos - Javeriana Cali/3. Información Geográfica y Espaciales/BD Aguacate_2020.xlsx")
# Se muestran las variables de la base de datos
str(BD_Aguacate_2020)## tibble [534 × 21] (S3: tbl_df/tbl/data.frame)
## $ id_arbol : num [1:534] 1 2 3 4 5 6 7 8 9 10 ...
## $ Latitude : num [1:534] 2.39 2.39 2.39 2.39 2.39 ...
## $ Longitude : num [1:534] -76.7 -76.7 -76.7 -76.7 -76.7 ...
## $ FORMATTED_DATE_TIME : chr [1:534] "01/10/2020 10:11:12 a, m," "01/10/2020 10:11:12 a, m," "01/10/2020 10:11:12 a, m," "01/10/2020 10:11:12 a, m," ...
## $ Psychro_Wet_Bulb_Temperature : num [1:534] 22 21.4 21.8 22.8 22.6 21.5 22.2 22.6 23 22.7 ...
## $ Station_Pressure : num [1:534] 825 825 826 825 825 ...
## $ Relative_Humidity : num [1:534] 85.2 84 79.6 77.6 76.5 77.7 76.5 77.7 78.3 77.8 ...
## $ Crosswind : num [1:534] 0 0 0.2 0.4 0 0.3 0 0 0 0.3 ...
## $ Temperature : num [1:534] 23.9 23.5 24.5 25.9 26 24.5 25.5 25.7 26 25.9 ...
## $ Barometric_Pressure : num [1:534] 825 825 826 825 825 ...
## $ Headwind : num [1:534] 0 0 0.4 0.2 0 0 0 0 0 0.3 ...
## $ Direction_True : num [1:534] 313 317 338 299 265 265 210 304 284 315 ...
## $ Direction_Mag : num [1:534] 312 317 337 299 264 264 209 303 284 315 ...
## $ Wind_Speed : num [1:534] 0 0 0.5 0.5 0 0.3 0 0 0 0.4 ...
## $ Heat_Stress_Index : num [1:534] 25.3 24.8 25.7 28.1 28 25.5 27.3 27.5 28.4 28.1 ...
## $ Altitude : num [1:534] 1696 1696 1694 1694 1696 ...
## $ Dew_Point : num [1:534] 21.3 20.7 20.8 21.7 21.5 20.4 21.1 21.5 22 21.7 ...
## $ Density_Altitude : num [1:534] 2.5 2.48 2.52 2.57 2.58 ...
## $ Wind_Chill : num [1:534] 23.9 23.5 24.5 25.9 25.9 24.5 25.5 25.6 26 25.9 ...
## $ Estado_Fenologico_Predominante: num [1:534] 717 717 717 717 717 717 717 717 717 717 ...
## $ Frutos_Afectados : num [1:534] 0 0 0 0 0 0 0 0 0 0 ...Se construye un mapa a través de Leaflet donde se logra determinar la ubicación exacta de los árboles en el municipio de Tibío.
# Se construye el mapa para la ubicación geoespacial de los Árboles.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = BD_Aguacate_2020$Longitude,lat = BD_Aguacate_2020$Latitude,radius = 0.2, color = "blue")Análisis Geoestadístico
# Se crea la variable regionalizada
df_aguacate_geo <- as.geodata(BD_Aguacate_2020, coords.col = 2:3, data.col = 9)
plot(df_aguacate_geo)
Con la creación de la variable regionalizada, a partir de la información correspondiente a la temperatura de los árboles de Aguacate, como se observa en el primer gráfico los puntos en rojo determinan aquellos con mayor temperatura, seguidos de los amarillos y verdes, mientras que los azules obtienen los niveles más bajos de este indicador. En toda la región existen areas con temperaturas altas y bajas, por tanto, no se observa una concentración marcada de esta característica.
# Se construye la matriz de distancias
summary(dist(BD_Aguacate_2020[,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-03Con la matriz se obtienen los rangos de la distancia los cuales se usarán como referente para la construcción del Semivariograma, tomando un valor superior al mínimo y otro por inferior al tercer cuartil.
# Construcción del Semivariograma.
variograma=variog(df_aguacate_geo, option = "bin",uvec=seq(0,0.0008978,0.00001912))## variog: computing omnidirectional variogramplot(variograma)
datos_env = variog.mc.env(df_aguacate_geo, 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 envelops# lines(datos_env)
plot(variograma, pch=18, envelop = datos_env)
Como se aprecia en el gráfico, los árboles tienen una alta similitud en sus niveles de temperatura, dentro de la banda se observa que hay un número importante de semivariogramas que no presentan correlación, pero a su vez existe tambien un grupo importante de semivariogramas que presentan correlación espacial.
Modelo de Semivarianza
Se procede a la construcción y ajuste del Modelo.
##Ajuste del Modelo de Semivarianza
ini.vals = expand.grid(seq(6,12,l=10), seq(5e-04,8e-04,l=10))
plot(variograma)
# Se construye el modelo exponencial
model_mco_exp=variofit(variograma, ini=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 "6" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 44365.7832119024# Se construye el modelo gaussiano
model_mco_gaus=variofit(variograma, ini=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 "6" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 156706.740941865# Se construye el modelo esférico
model_mco_spe=variofit(variograma, ini=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 "6" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 391553.149802868lines(model_mco_exp,col="blue")
lines(model_mco_gaus,col="red")
lines(model_mco_spe,col="green")
legend(0.0007, 7, legend=c("Gausiano", "Exponencial", "Esferico"),
col=c("red", "blue","green"), lty=1, cex=0.8)
En la gráfica se observa que el modelo que presenta el mejor ajuste es el exponencial, lo anterior, se comprueba calculando el MSE:
print(paste("MSE modelo Exponencial =", as.character(model_mco_exp$value)))## [1] "MSE modelo Exponencial = 44334.6829205354"print(paste("MSE modelo Gaussiano =", as.character(model_mco_gaus$value)))## [1] "MSE modelo Gaussiano = 156377.750582178"print(paste("MSE Modelo Esférico =", as.character(model_mco_spe$value)))## [1] "MSE Modelo Esférico = 83111.7273276977"El modelo con el valor más bajo de MSE es el modelo exponencial con un valor de 44.334,68, por tanto sobre este se realizará la predicción espacial.
Predicción Espacial
max(BD_Aguacate_2020[,3])## [1] -76.71022# Predicción espacial Kriging
datos_grid = expand.grid(Latitude = seq(2.392,2.394, l=100),Longitud = seq(-76.7118,-76.71022, l=100))datos_ko = krige.conv(df_aguacate_geo, loc = datos_grid, krige = krige.control(nugget = 0, trend.d = "cte", trend.l = "cte", cov.pars = c(sigmasq = 3.1344, phi = 0.0001)))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodimage(datos_ko, main = "Kriging Predict - Temperatura", xlab ="Latitude", ylab ="Longitude")
Como se puede observar en el gráfico geoestadístico, las regiones con el color rojo más oscuro representan la ubicación de los árboles de aguacate que registraron la mayor temperatura.