Ejercicio: Clima Finca de Aguacate

Actividad

  • Hacer el geo data de una variable climática presente en la base de datos de la finca de aguacates.

  • Obtener el semivariograma.

  • Realizar el ajuste del modelo teórico

  • Predicción espacial con Kriging

Primero se realiza el cargue de las librerías a necesitar

#cargar librerias
library(readxl)
library(geoR)
library(raster)
library(ggplot2)
library(leaflet)
library(lubridate)
library(sp)

Como se va a filtrar por la fecha 01-10-2020, se decide crear una variable que contenga separada la fecha y la hora, y así obtener el dataframe filtrado.

También se revisa que se tenga las coordenadas geográficas y la variable temperatura por cada árbol en la finca.

#se lee el excel con los datos
aguacates = read_excel("C:/Users/dolli/Documents/ESTUDIO/Maestria ciencia de datos/Segundo_semestre/electiva/Datos_Completos_Aguacate.xlsx") 

#se convierte la columna de fecha en formato date
aguacates$fecha = as.Date(aguacates$FORMATTED_DATE_TIME, format = "%d-%m-%Y %H:%M")

#se crea la columna fecha_1 para separar fecha de hora
fecha_hora = dmy_hms(aguacates$FORMATTED_DATE_TIME)
hora = hour(fecha_hora)
fecha_1 = date(fecha_hora)
aguacates$fecha_1 = fecha_1

#se filtra por la fecha 01/10/2020
aguacates_2020 = aguacates[aguacates$fecha_1 == as.Date("2020-10-01"), ]
names(aguacates_2020)
##  [1] "id_arbol"                       "Latitude"                      
##  [3] "Longitude"                      "FORMATTED_DATE_TIME"           
##  [5] "Psychro_Wet_Bulb_Temperature"   "Station_Pressure"              
##  [7] "Relative_Humidity"              "Crosswind"                     
##  [9] "Temperature"                    "Barometric_Pressure"           
## [11] "Headwind"                       "Direction_True"                
## [13] "Direction_Mag"                  "Wind_Speed"                    
## [15] "Heat_Stress_Index"              "Altitude"                      
## [17] "Dew_Point"                      "Density_Altitude"              
## [19] "Wind_Chill"                     "Estado_Fenologico_Predominante"
## [21] "Frutos_Afectados"               "fecha"                         
## [23] "fecha_1"

Mapas de árboles ubicados en la finca de aguacate

La finca se encuentra en la zona rural del municipio de Popayan y su distribucción espacial es bastante regular.

#se grafican los puntos obtenidos
plot(aguacates_2020[,3:2])
title("Ubicación de árboles de la finca de aguacates")

leaflet() %>% addTiles() %>% addCircleMarkers(lng = aguacates_2020$Longitude, 
                                              lat = aguacates_2020$Latitude,
                                              radius = 0.2, color = "green") %>%
  addControl("Ubicación de árboles de la finca de aguacates", position = "topright")

Geodata obtenida para la variable temperatura de la finca aguacate

Como primer paso realizamos la creación de la variable regionalizada:

#se pasa aguacates_2020 a formato geo con variable temperatura
aguacates_geo = as.geodata(aguacates_2020, coords.col = 3:2, data.col=9)
plot(aguacates_geo)

Semivariograma para temperatura de la finca aguacate

Con el semivariograma se observa una autocorrelación espacial indicando que las mediciones cercanas, presentan temperaturas similares en contraste con las distantes.

También se genera un intervalo de confianza en donde garantizamos que nuestro semivariograma muestra un patrón y no son datos aleatorios.

# obtener la matriz de distancia
semi_aguacate = variog(geodata = aguacates_geo, option = "bin", uvec = seq(0, 0.00098, 0.00005))
## variog: computing omnidirectional variogram
plot(semi_aguacate)
# generar intervalo de confianza para el variograma
variograma_ag = variog.mc.env(aguacates_geo, obj.variog = semi_aguacate)
## 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(variograma_ag)

Ajuste del modelo teórico

Se realiza el ajuste del modelo teorico:

Aqui cada linea en el semivariograma representa un modelo,

  • linea azul es para el modelo exponencial

  • linea purpura es para el modelo gaussiano

  • linea roja es para el modelo spherical

#Ajuste del modelo de semivarianza

ini.vals = expand.grid(seq(2.5,3.5,l=10), seq(6e-04,8e-04,l=10)) # para temperatura
plot(semi_aguacate)

model_mco_exp = variofit(semi_aguacate, 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 "3.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 134058.673070192
model_mco_gaus = variofit(semi_aguacate, 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 "3.5"   "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 202567.144389584
model_mco_spe = variofit(semi_aguacate, ini = ini.vals, cov.model = "spherical", 
                         fix.nugget = 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 "3.39"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 22755.7548444676
plot(semi_aguacate)

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

Se observa gráficamente que el modelo exponencial es el que mejor se ajusta a los puntos del semivariograma muestral, y esto se confirma con los valores de SCE del modelo.

Predicción espacial con Kriging

Se realiza la interpolación (predicción espacial)

# prediccion Espacial Kriging

datos_grid = expand.grid(lon=seq(-76.7120, -76.7100, l=100), lat=seq(2.3920,2.3938, l=100))
plot(datos_grid)
points(aguacates_2020[,3:2], col="red")

#prediccion para temperatura con modelo exponencial
datos_ko = krige.conv(aguacates_geo, loc=datos_grid,
                      krige = krige.control(nugget = 0, trend.d = "cte",
                                            trend.l = "cte", cov.pars = c(sigmasq=4.6320, phi=0.0028)))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
image(datos_ko, main = "kriging predict", xlab="longitud", ylab="latitud")
contour(datos_ko, main="kriging predict", add=TRUE, drawlabels=TRUE)

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

Imagen ráster con la predicción de temperatura para la finca de aguacates

Se tranforma la imagen en ráster:

#transformar en imagen raster
pred = cbind(datos_grid, datos_ko$predict)
temp_predict = rasterFromXYZ(cbind(datos_grid, datos_ko$predict))
plot(temp_predict, main = "Ráster temperatura predict")

Utilizando el método convex hull, obtenemos un polígono que aproxima el contorno del conjunto de puntos correspondientes a los árboles de aguacate. Luego, aplicamos este polígono como máscara para recortar la imagen ráster de la predicción de temperatura.

#calcular los indices del convex hull
hull_indices = chull(aguacates_2020$Longitude, aguacates_2020$Latitude)

#crear un poligono a partir del convex
hull_coords = aguacates_2020[c(hull_indices, hull_indices[1]), c("Longitude", "Latitude") ]
poligono_hull = sf::st_polygon(list(as.matrix(hull_coords)))

#transformo el poligno en sf y luego a spatial para poder recortarlo
poligono_sf = sf::st_sfc(poligono_hull)
poligono_sp = as(poligono_sf, "Spatial")
raster_cut = mask(temp_predict, poligono_sp)

#visualizacion 
plot(raster_cut, main = "Ráster temperatura predict finca")
plot(poligono_sp, add = TRUE, border = "red", lwd = 2)

Conclusiones

En la predicción de temperatura se observa que la parte sur de la finca presenta un sector con temperaturas alrededor de los 29 grados Celsius. Esto podría atribuirse a la ubicación en ladera de la finca, que hace que esta área esté más cerca del nivel del mar, o bien a factores como un manejo efectivo de las podas o a que los árboles en esta zona sean aún jóvenes, lo que reduce su capacidad de influir en la creación de un microclima. En contraste, en el sector norte la temperatura estimada varía entre 24 y 26 grados Celsius, lo que podría indicar que esta parte de la finca está a mayor altitud y cuenta con árboles de copas más densas. Además, en el lado noreste se identifica una línea con temperaturas más bajas, cercanas a los 23 grados Celsius, lo cual podría deberse a una elevación aún mayor en el terreno o a un manejo inadecuado de las podas, que fomenta la creación de un microclima más fresco en esta área.

En conclusión, los métodos de interpolación espacial permiten predecir y analizar la dinámica entre variables dentro de un área específica, proporcionando una comprensión más detallada de la variación espacial de la temperatura en la finca.