Cultivo de Aguacate
Juan Fernando Gutiérrez Ayala
# Instalación y cargue de librerías
library(DT)
library(kableExtra)
library(readxl)
library(dplyr)
library(ggplot2)
library(leaflet)
library(tidyverse)
library(MASS)
library(geoR)
library(rasterVis)
library(raster)
library(knitr)1 Introducción
El presente ejercicio tiene como objetivo la identificación y predicción de zonas que, según su tempertaura, son más aptas para el cultivo de aguacate, partiendo de la ubicación de los árboles de manera específica.
2 Preparación de datos
Se cargan y presenta la estructura de los datos a utilizar en el análisis, la cual contiene un total de 21 atributos que caracterizan a 20271 entradas
# Cargue de base de datos
datos_aguacate <- read_excel("C:/Users/Juan/Desktop/Maestría/Casos - Información Geográfica/Datos_Completos_Aguacate.xlsx")
head(datos_aguacate)Así mismo, dado que el horizonte temporal de interés responde al periodo de recolección de datos de 01/10/2020, se realiza el filtro, derivando en un total de 534 registros.
#Filtro a fecha requerida
# Eliminar la hora y los segundos del formato "21/08/2019 9:36:36 a. m."
datos_aguacate$FORMATTED_DATE_TIME <- substr(datos_aguacate$FORMATTED_DATE_TIME, 1, 10)
# Convertir a fecha
datos_aguacate$FORMATTED_DATE_TIME <- as.Date(datos_aguacate$FORMATTED_DATE_TIME, format = "%d/%m/%Y")
datos_aguacate_filt <- filter(datos_aguacate, FORMATTED_DATE_TIME == as.Date("2020-10-1"))
df<-datos_aguacate_filt3 Ubicación de los árboles de aguacate.
4 Análisis geoestadístico
Con los datos filtrados, se realiza un análisis de la temperatura en función de la ubicación de los árboles de la finca estudiada.
#Análisis geoestadístico de la variable temperatura
df_geo <- as.geodata(df, coords.col = 3:2, data.col = 9) # Variable regionalizada
plot(df_geo)
Se aprecia un patrón en el que la temperatura a la que está expuesto cada árbol responde a su ubicación, ya que presentan homogeneidad en el sentido que se aprecian agrupaciones de árboles con proximidad espacial y exposición a temperaturas iguales.
5 Variograma
En principio, se genera un resumen del comportamiento de las distancias para definir parámetros específicos en la construcción de los variogramas.
## 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 esto, se generan los variogramas a partir de la información identificada.
## variog: computing omnidirectional variogram## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
Como se puede ver, existe una correlación entre la temperatura y la distancia a la que se encuentran los árboles entre ellos, ya que la mayoría de puntos se encuentran fuera de la franja del intervalo de confianza, revelando que responden a un comportamiento que no es aleatorio.
6 Ajuste del modelo.
A continuación de realiza un ajuste de modelos teóricos al modelo derivado de los datos estudiados, haciendo enfasis en el exponencial, gaussiano y esférico.
#Ajuste del modelo
ini.vals <- expand.grid( seq(2,4,l=10), seq(0.0001,0.00025, l = 10))
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.11" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 3555.92876863577## 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.11" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 6633.5155506537## 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.11" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 6005.01775448873plot(variograma, pch=16, envelop = datos_env)
lines(mco_exp, col = "blue")
lines(mco_gau, col = "red")
lines(mco_sph, col = "green")
legend(0.00055, 6.9, legend=c("Gaussiano", "Exponencial", "Esférico"),
col=c("red", "blue","green"), lty=1, cex=0.7)
Cómo se puede ver gráficamente, el modelo exponencial se ajusta de manera un poco más fiel al comportamiento de los datos evidenciado a través del variograma realizado. Sin embargo, se confirma dicha decisión a partir del cálculo de la suma de los cuadrados de los errores para definir el más pertinente para este caso, confirmando que será el modelo exponencial el más eficiente.
# Cálculo de mejor modelo (suma de los cuadrados de los errores)
SCE_exp <- mco_exp$value
SCE_gaus <- mco_gau$value
SCE_sph <- mco_sph$value
SCE <- c(Exponencial = SCE_exp, Gaussiano = SCE_gaus, Esférico = SCE_sph)
print(SCE)## Exponencial Gaussiano Esférico
## 3314.176 6339.978 5837.5897 Predicción espacial
En un primer momento, se definen los límites de la zona a estudiar en esta predicción.
min_lat <- min(df_geo$coords[, 1]) # La primera columna es latitud
max_lat <- max(df_geo$coords[, 1])
min_long <- min(df_geo$coords[, 2]) # La segunda columna es longitud
max_long <- max(df_geo$coords[, 2])Con esto, se contruye la grilla sobre la cual se va a realizar la predicción a partir del modelo exponencial.
geodatos_grid<-expand.grid( lon=seq(-76.710215,-76.711799,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(df[,3:2],col="red")
Así, con el sigma (entendida como la varianza del proceso) y phi (rango de correlación espacial) arrojados por el modelo exponencial, se hace el kriging para predeccir el comportamiento de la temperatura en el sector priorizado.
geodatos_ko<-krige.conv(df_geo, loc=geodatos_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 neighbourhoodpar(mar = c(2, 2, 3, 1))
image(geodatos_ko, main = "Kriging Predict - Temperatura", xlab ="Latitude", ylab ="Longitude")
Finalmente, se visualiza el resultado de la predicción en imagen raster para una lectura más apropiada de las temperaturas esperadas en cada zona de la finca
pred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
levelplot(temp_predict,
par.settings=BuRdTheme,
scales=list(x=list(rot=45),
cex=1.2),
margin=FALSE,
asp=1) 