Actividad 3 - Exploración Datos Geoestadísticos
Juan Camilo López
2024-11-23
Introducción
Este informe presenta el análisis de datos geoespaciales de una finca de aguacates en el departamento del Cauca, en la zona rural del municipio de Timbío. En particular, se analizan los datos de temperatura tomados por sensores dispuestos en los árboles de la finca. A través de este ejercicio se quiere construir un mapa donde se muestren las predicciones de la temperatura en la finca de aguacate utilizando la metodología “Kriging”. El documento está dividido en cuatro partes. Primero se exploran los datos de interés por medio de técnicas de estadística descriptiva, en segundo lugar, se presenta la exploración de datos por medio de la Geoestadística, donde se plantea la exploración de la correlación espacial de los datos por medio de la construcción de un semivariograma y se ajustan tres modelos para escoger el más apropiado. Finalmente, se utiliza el modelo más apropiado para hacer el mapa de la predicción de los valores de la temperatura y se presentan algunas conclusiones.
# Cargar datos
aguacates <- read_excel("Datos_Completos_Aguacate.xlsx")
aguacates_filtro = aguacates %>%
filter (FORMATTED_DATE_TIME == "01/10/2020 10:11:12 a, m,")
glimpse(aguacates_filtro) ## Rows: 534
## Columns: 21
## $ id_arbol <chr> "1", "2", "3", "4", "5", "6", "7", "8",…
## $ Latitude <dbl> 2.393549, 2.393573, 2.393541, 2.393503,…
## $ Longitude <dbl> -76.71124, -76.71120, -76.71113, -76.71…
## $ FORMATTED_DATE_TIME <chr> "01/10/2020 10:11:12 a, m,", "01/10/20…
## $ Psychro_Wet_Bulb_Temperature <dbl> 22.0, 21.4, 21.8, 22.8, 22.6, 21.5, 22.…
## $ Station_Pressure <dbl> 825.1, 825.3, 825.5, 825.4, 825.2, 825.…
## $ Relative_Humidity <dbl> 85.2, 84.0, 79.6, 77.6, 76.5, 77.7, 76.…
## $ Crosswind <dbl> 0.0, 0.0, 0.2, 0.4, 0.0, 0.3, 0.0, 0.0,…
## $ Temperature <dbl> 23.9, 23.5, 24.5, 25.9, 26.0, 24.5, 25.…
## $ Barometric_Pressure <dbl> 825.2, 825.2, 825.5, 825.4, 825.2, 825.…
## $ Headwind <dbl> 0.0, 0.0, 0.4, 0.2, 0.0, 0.0, 0.0, 0.0,…
## $ Direction_True <dbl> 313, 317, 338, 299, 265, 265, 210, 304,…
## $ Direction_Mag <dbl> 312, 317, 337, 299, 264, 264, 209, 303,…
## $ Wind_Speed <dbl> 0.0, 0.0, 0.5, 0.5, 0.0, 0.3, 0.0, 0.0,…
## $ Heat_Stress_Index <dbl> 25.3, 24.8, 25.7, 28.1, 28.0, 25.5, 27.…
## $ Altitude <dbl> 1696, 1696, 1694, 1694, 1696, 1698, 169…
## $ Dew_Point <dbl> 21.3, 20.7, 20.8, 21.7, 21.5, 20.4, 21.…
## $ Density_Altitude <dbl> 2.504, 2.485, 2.518, 2.572, 2.575, 2.52…
## $ Wind_Chill <dbl> 23.9, 23.5, 24.5, 25.9, 25.9, 24.5, 25.…
## $ Estado_Fenologico_Predominante <dbl> 717, 717, 717, 717, 717, 717, 717, 717,…
## $ Frutos_Afectados <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …Análisis exploratorio
Observemos la ubicación de los árboles en el mapa.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = aguacates_filtro$Longitude,lat = aguacates_filtro$Latitude,radius = 0.2,color = "blue")Podemos ver que los árboles están reunidos en un sólo espacio de la finca, cerca del municipio de Timbío. Ahora exploremos la distribución de la temperatura en esta muestra.
ggplot(aguacates_filtro, aes(x = Temperature)) +
geom_histogram(binwidth = 1, fill = "lightblue", color = "black") +
theme_minimal() +
labs(title = "Distribución de la Temperatura", x = "Temperatura", y = "Frecuencia")
Según el histograma, la temperatura tiene una distribución similar a la normal, con la mayoría de los puntos donde se tomó la temperatura alrededor de los 26 grádos, en un rango entre los 22 y 30 grados.
Análisis Geoestadístico
Para realizar los anális geoestadísticos se transforman los datos en una variable regionalizada, es decir de tipo Geodata.
geo_aguacates = as.geodata(aguacates_filtro,coords.col = 3:2,data.col = 9)
plot(geo_aguacates)
En el primer gráfico de la esquina superior izquierda podemos ver que los círculos azules corresponden a las temperaturas más bajas, mientras que las équis corresponden a las más altas, encontrando un indício de un patrón para la distribución de la temperatura. Adicionalmente, en el histograma vemos nuevamente que el rango de la temperatura está entre los 22 y los 30 grados.
Creación semivariograma
El semivariograma nos permitirá determinar si existe una covarianza entre la temperatura y la ubicación geográfica de los árboles, presentando una gráfica que nos muestre cómo cambia la semivarianza a medida que los puntos observados tiene mayor distancia. Para realizar este gráfico es necesario conocer la distribución de la distancia entre cada árbol:
summary(dist(geo_aguacates$coords))## 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-03Una vez se conoce la distribución de los valores de las distancias, podemos usar la distancia mínima entre ellos y el tercer cuartil como parámetros para construír el semivariograma. Adicionalmente, para ayudar en la interpretación del gráfico, se generan simulaciones de semivariogramas si las distribuciones fueran totalmente aleatorias. Las distribuciones se representan por dos líneas, si la distribución del semivariograma resulta estar en su mayoría en ese rango, podemos decir que no existe una correlación entre la distribución en el espacio y la variable temperatura.
semi_aguacates = variog(geo_aguacates,option = "bin", uvec=seq(0,0.001,0.00002))## variog: computing omnidirectional variogramdatos.env=variog.mc.env(geo_aguacates,obj=semi_aguacates)## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelopsplot(semi_aguacates)
lines(datos.env)
Como se puede observar en el gráfico, aunque hay algunos casos que caen dentro del rango determinado para distribuciones aleatorias, el semivariograma para nuestro caso puede considerarse válido para determinar que existe una correlación entre la distancia de los árboles y la temperatura. Además, la distribución de los puntos sigue el patrón que normalmente se espera, mayor similitud a menor distancia o, menor semivarianza a menor distancia.
Ajuste de modelos
Ahora, con el objetivo de encontrar un modelo que mejor se ajuste a la distribución encontrada y nos permita predecir la variable objetivo “temperatura”, se probarán tres modelos y se eligirán los más adecuados
Modelo exponencial
ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,0.0008,l=10))
model_agu_exp=variofit(semi_aguacates, 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: 296807.231725087model_agu_exp## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 0.7969 2.2860 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.000232806
##
## variofit: minimised weighted sum of squares = 9158.212Modelo Gausiano
model_agu_gaus=variofit(semi_aguacates, 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.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 288205.86784414model_agu_gaus## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
## tausq sigmasq phi
## 0.7902 2.2842 0.0001
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
##
## variofit: minimised weighted sum of squares = 8947.68Modelo esférico
model_agu_spe=variofit(semi_aguacates, 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.5" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 283743.808438713model_agu_spe## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq = 0
## parameter estimates:
## sigmasq phi
## 3.0186 0.0000
## Practical Range with cor=0.05 for asymptotic range: 0
##
## variofit: minimised weighted sum of squares = 17731.76Para determinar cuál de los tres modelos se ajusta mejor a nuestro caso debemos observar la diferencia de cuadrados, que para el caso corresponde al modelo gausiano, con 8,947. Veamos los tres modelos frente a la distribución de los datos en el semivariograma.
plot(semi_aguacates)
lines(model_agu_exp,col="blue")
lines(model_agu_gaus,col="red")
lines(model_agu_spe,col="green")
Al conocer el modelo que mejor se ajusta podemos proceder a realizar el mapa con la predicción de los valores, empleando el modelo Gausiano.
Predicción de valores de temperatura
Primero determinamos la malla o grid sobre la que se realizarán las predicciones, determinada por las coordenadas dentro del dataset.
c(min(aguacates_filtro[,3]),
max(aguacates_filtro[,3]),
min(aguacates_filtro[,2]),
max(aguacates_filtro[,2]))## [1] -76.711799 -76.710215 2.392101 2.393634geodatos_grid = expand.grid( lon=seq(-76.710215,-76.711799,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_aguacates$coords,col="red")
Mapa de predicción
Una vez conocemos el espacio donde se implementará el método Kriging podemos generar los mapas.
geodatos_kr=krige.conv(geo_aguacates, loc=geodatos_grid,
krige= krige.control(nugget=0,trend.d="cte",
trend.l="cte",cov.pars=c(sigmasq=2.2842, phi=0.0001 )))## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodpar(mfrow=c(1,2))
contour(geodatos_kr,main="kriging Predict", drawlabels=TRUE)
image(geodatos_kr, main="kriging Predict", xlab="East", ylab="North")
Conclusiones
Al observar los mapas generados es claro que las mayores varianzas están dentro del espacio que ocupan los árboles de aguacates, que es precizamente de dónde se extrajo la información para el ejercicio. De este método pueden sacarse varias utilizadaes, para el caso del aguacate por ejemplo, podría determinarse el terreno más adecuado para el cultivo de árboles que sean más productivo e incluso tengan mejores frutos.