#Objetivo

Analizar la variabilidad espacial de la temperatura en un cultivo de aguacates usando métodos geoestadísticos para generar un mapa de predicción

#Preparación del Espacio de Trabajo

En este primer capitulo, nuestro objetivo será preparar el espacio de trabajo, para ello usaremos las siguientes librerías con sus respectivos usos:

#install.packages(c("sf", "gstat", "ggplot2", "viridis", 
                  # "dplyr", "readr", "rmarkdown", "rsconnect"))

#Carga y Preparación de Datos

Dentro de nuestro proceso para la carga y preparación de los datos, solo usaremos la información de la fecha 01/10/2020, para ello procederemos a leer las respectivas librerías, lectura del archivo base de la data, con la selección de la variable de interés de temperatura y finalmente un análisis exploratorio de datos incial.

Para ello transformaremos el sistema de coordenadas proyectadas UTM18N, EPSG:32618 buscando lograr medir las distancias en metros

# --- Carga de Librerías ---

library(readr)       
library(dplyr)       
library(sf)          
library(gstat)       
library(ggplot2)    
library(viridis)  

# --- Carga de Datos ---


file_path <- "C:/Maestria/Semestre 2/Analisis de Datos Geográficos y Espaciales/Módulo 3/Actividad/Insumo a Evaluar.csv"

#file_path <- "https://drive.google.com/file/d/1azvCBFuOBaMXLvOCd-QSDYbEfoWM0nT-/view?usp=drive_link
#datos <- read_csv(file_path)"

datos <- read_csv(file_path)

names(datos)
##  [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"
# --- Limpieza y Selección ---

# Seleccionamos las columnas de interés

datos_limpios <- datos %>%
  
  select(id_arbol, Latitude, Longitude, Temperature)

# --- Creación de Objeto Espacial (sf) ---

# 1. Convertimos el data.frame a un objeto 'sf'

# 2. Asignamos el CRS original (WGS84, EPSG:4326), que es el estándar para Lat/Lon.

datos_sf_geo <- st_as_sf(datos_limpios, 
                         coords = c("Longitude", "Latitude"), 
                         crs = 4326)

# 3. Transformamos a un CRS Proyectado (UTM Zona 18N, EPSG:32618), buscando lograr medir distancias en metros, ya que es un modelo plano


datos_sf_utm <- st_transform(datos_sf_geo, crs = 32618)

# --- Análisis Exploratorio Inicial (EDA) ---

# Visualicemos los puntos de muestreo y su temperatura

ggplot(data = datos_sf_utm) +
  geom_sf(aes(color = Temperature), size = 2) +
  scale_color_viridis(option = "C", name = "Temperatura (°C)") +
  labs(title = "Distribución Espacial de la Temperatura",
       subtitle = "Muestreo de Árboles (Proyección UTM 18N)",
       x = "Coordenada Este (metros)",
       y = "Coordenada Norte (metros)") +
  theme_minimal()

# Cálculo del Semivariograma

Modelaremos la autocorrelación espacial, primero, calculamos el variograma empirico, para ello realizaremos el cálculo de semivariograma empírico usando la variable temperatura

# --- Cálculo del Semivariograma Empírico ---

# la temperatura como una función de su posición (Kriging Ordinario).

vario_empirico <- variogram(Temperature ~ 1, data = datos_sf_utm)

# Visualizamos el variograma empírico

plot(vario_empirico, 
     main = "Semivariograma Empírico de Temperatura",
     xlab = "Distancia (h) en metros",
     ylab = "Semivarianza")

Una vez realizado el Semivariograma , procederemos a realizar el ajuste de modelos teóricos y seleccionaremos el mejor ajuste basándonos en el menor Error de Suma de cuadrados SSE

# --- Ajuste de Modelos Teóricos ---

# Probaremos los tres modelos clásicos.

# 1. Modelo Esférico (Sph)

vgm_esferico <- fit.variogram(vario_empirico, model = vgm("Sph"))
plot(vario_empirico, model = vgm_esferico, main = "Ajuste Modelo Esférico (Sph)")

print("Modelo Esférico:")
## [1] "Modelo Esférico:"
print(vgm_esferico)
##   model    psill    range
## 1   Nug 0.753057  0.00000
## 2   Sph 2.149740 30.80213
# 2. Modelo Exponencial (Exp)

vgm_exponencial <- fit.variogram(vario_empirico, model = vgm("Exp"))
plot(vario_empirico, model = vgm_exponencial, main = "Ajuste Modelo Exponencial (Exp)")

print("Modelo Exponencial:")
## [1] "Modelo Exponencial:"
print(vgm_exponencial)
##   model     psill    range
## 1   Nug 0.4365545  0.00000
## 2   Exp 2.6356490 13.57746
# 3. Modelo Gaussiano (Gau)

vgm_gaussiano <- fit.variogram(vario_empirico, model = vgm("Gau"))
plot(vario_empirico, model = vgm_gaussiano, main = "Ajuste Modelo Gaussiano (Gau)")

print("Modelo Gaussiano:")
## [1] "Modelo Gaussiano:"
print(vgm_gaussiano)
##   model    psill    range
## 1   Nug 1.238174  0.00000
## 2   Gau 1.657167 14.84052
# --- Selección del Mejor Modelo ---

# 'fit.variogram' añade un atributo 'SSE' (Suma de Cuadrados del Error).

# El modelo con el menor SSE es, estadísticamente, el de mejor ajuste.

sse_sph <- attr(vgm_esferico, "SSE")
sse_exp <- attr(vgm_exponencial, "SSE")
sse_gau <- attr(vgm_gaussiano, "SSE")

cat(paste("SSE Esférico:", sse_sph, "\n"))
## SSE Esférico: 1.05537139848762
cat(paste("SSE Exponencial:", sse_exp, "\n"))
## SSE Exponencial: 1.11254445202679
cat(paste("SSE Gaussiano:", sse_gau, "\n"))
## SSE Gaussiano: 1.68073773507839
# Seleccionamos el modelo final

modelo_final <- vgm_esferico

La suma de cuadrados del error para los 3 modelos, el esferico, exponencial y gaussiano, donde podemos observar el menor error en el esférico

Predicción Espacial - Kriging Ordinario

Procederemos a interpolar los valores en una cuadricula (grid), con una celda 2,2 metros buscando cubrir toda la parcela, para ello usaremos el modelo de variograma seleccionado para predecir la temperatura

# --- Creación de la Cuadrícula (Grid) de Predicción ---


grid_prediccion <- st_make_grid(datos_sf_utm, 
                                cellsize = c(2, 2), 
                                what = "centers")

# Convertimos el grid a un objeto 'sf'

grid_sf <- st_sf(geometry = grid_prediccion)

# --- Ejecución del Kriging Ordinario ---

# Usamos gstat::krige para predecir 'Temperature ~ 1'

# en las 'newdata' (el grid), usando nuestro 'modelo_final'.

kriging_resultado <- krige(
  formula = Temperature ~ 1,
  locations = datos_sf_utm,
  newdata = grid_sf,
  model = modelo_final
)
## [using ordinary kriging]
# El resultado 'kriging_resultado' tiene dos columnas importantes:

# 1. 'var1.pred': La predicción de temperatura.
# 2. 'var1.var': La varianza del error de predicción.

#Generación de la Imagen de Proyección (Mapa)

Visualizaremos los resultados del Kriging, para crear 2 tipos de Mapa:

# --- Conversión del resultado a data.frame para ggplot2 ---


kriging_df <- as.data.frame(kriging_resultado)
kriging_df <- cbind(kriging_df, st_coordinates(kriging_resultado$geometry))

# --- Mapa de Predicción de Temperatura (Imagen 1) ---

ggplot() +
  # Capa del raster de predicción
  
  geom_raster(data = kriging_df, aes(x = X, y = Y, fill = var1.pred)) +
  
  # Capa de los puntos de muestreo originales
  
  geom_sf(data = datos_sf_utm, color = "black", size = 1.5, shape = 21, fill = "white") +
  
  # Paleta de color y etiquetas
  
  scale_fill_viridis(option = "C", name = "Temp. Predicha (°C)") +
  labs(title = "Mapa de Predicción Kriging de Temperatura",
       subtitle = "Estimación de temperatura en el lote de aguacates",
       x = "Coordenada Este (metros)",
       y = "Coordenada Norte (metros)") +
  theme_minimal() +
  coord_sf(datum = st_crs(datos_sf_utm)) # Asegura la escala correcta

# --- Mapa de Varianza del Error ---

ggplot() +
  # Capa del raster de varianza
  
  geom_raster(data = kriging_df, aes(x = X, y = Y, fill = var1.var)) +
  
  # Capa de los puntos de muestreo originales
  
  geom_sf(data = datos_sf_utm, color = "black", size = 1.5, shape = 21, fill = "white") +
  
  # Paleta de color y etiquetas
  
  scale_fill_viridis(option = "B", name = "Varianza") +
  labs(title = "Mapa de Varianza del Error (Kriging)",
       subtitle = "Incertidumbre de la predicción (zonas rojas = más incertidumbre)",
       x = "Coordenada Este (metros)",
       y = "Coordenada Norte (metros)") +
  theme_minimal() +
  coord_sf(datum = st_crs(datos_sf_utm))

Dentro de la interpretación del mapa de predicción de kriging de Temperatura, nos muestra zonas frias entre 23-25° que atraviesan el lote diagonalmente , esto sugiere influencia de acumulación de humedad, sobras parciales, puede llegar a ser por los árboles altos, barreras vivas, entre otras, asi mismo podemos ver zonas de mayor temperatura como es el caso del NorEste, esto puede llegar a ser asociado a zonas de mayor exposición solar, menor sombra vegetal, suelos con menor humedad a revisar.