Preparación de los datos

## Loading required package: readxl
## Loading required package: tidyverse
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.3     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: geoR
## 
## The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
## which was just loaded, will retire in October 2023.
## Please refer to R-spatial evolution reports for details, especially
## https://r-spatial.org/r/2023/05/15/evolution4.html.
## It may be desirable to make the sf package available;
## package maintainers should consider adding sf to Suggests:.
## The sp package is now running under evolution status 2
##      (status 2 uses the sf package in place of rgdal)
## 
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-4 (built on 2024-02-14) is now loaded
## --------------------------------------------------------------
## 
## 
## Loading required package: MASS
## 
## 
## Attaching package: 'MASS'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select
## 
## 
## Loading required package: knitr
## 
## Loading required package: kableExtra
## 
## 
## Attaching package: 'kableExtra'
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## 
## 
## Loading required package: leaflet
## 
## Loading required package: raster
## 
## Loading required package: sp
## 
## 
## Attaching package: 'raster'
## 
## 
## The following object is masked from 'package:MASS':
## 
##     select
## 
## 
## The following object is masked from 'package:dplyr':
## 
##     select

En primer lugar, vamos a importar la base de datos y dejar los registros que nos interesan, con la fecha del 01/10/2020 a las 10:11:12 a.m.

data <- read_excel("Datos_Completos_Aguacate.xlsx")
data <- data[data$FORMATTED_DATE_TIME == "01/10/2020  10:11:12 a, m,",]

La base de datos contiene 534 registros y 21 variables. Estas son id del árbol, latitud y longitud. Se registra la fecha y hora del muestreo y variables climáticas como la temperatura del bulbo húmedo, la presión atmosférica, la humedad relativa, y componentes del viento como el viento cruzado, el viento en contra, la dirección y la velocidad del viento. También se incluyen datos de temperatura, presión barométrica, y otros índices como el índice de estrés térmico, el punto de rocío, la altitud, densidad, y la sensación térmica del viento. Además, se registra el estado fenológico del árbol y el número o porcentaje de frutos afectados.

Verificamos que no existan valores faltantes en la base de datos.

data %>%summarise_all(funs(sum(is.na(.))))
## # A tibble: 1 × 21
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME Psychro_Wet_Bulb_Temperature
##      <int>    <int>     <int>               <int>                        <int>
## 1        0        0         0                   0                            0
## # ℹ 16 more variables: Station_Pressure <int>, Relative_Humidity <int>,
## #   Crosswind <int>, Temperature <int>, Barometric_Pressure <int>,
## #   Headwind <int>, Direction_True <int>, Direction_Mag <int>,
## #   Wind_Speed <int>, Heat_Stress_Index <int>, Altitude <int>, Dew_Point <int>,
## #   Density_Altitude <int>, Wind_Chill <int>,
## #   Estado_Fenologico_Predominante <int>, Frutos_Afectados <int>

Análisis Exploratorio de datos

#calcular estadisticas descriptivas
resumen <- t(summary(data[, 5:ncol(data)]))

kable(resumen,
      caption = "Tabla 1. Estadísticas Descriptivas",
      booktabs = TRUE,
      linesep = c("\\hline", "\\hline"),
      align = "c") %>%
  kable_styling(full_width = F)
Tabla 1. Estadísticas Descriptivas
Psychro_Wet_Bulb_Temperature Min. :19.00 1st Qu.:20.80 Median :21.50 Mean :21.63 3rd Qu.:22.40 Max. :24.90
Station_Pressure Min. :822.4 1st Qu.:824.8 Median :825.7 Mean :825.5 3rd Qu.:826.5 Max. :827.4
Relative_Humidity Min. :59.50 1st Qu.:67.33 Median :71.25 Mean :71.17 3rd Qu.:75.30 Max. :91.60
Crosswind Min. :0.0000 1st Qu.:0.0000 Median :0.1000 Mean :0.1989 3rd Qu.:0.3750 Max. :1.5000
Temperature Min. :22.20 1st Qu.:24.50 Median :25.80 Mean :25.83 3rd Qu.:27.18 Max. :29.70
Barometric_Pressure Min. :822.4 1st Qu.:824.7 Median :825.7 Mean :825.5 3rd Qu.:826.5 Max. :827.4
Headwind Min. :-0.7000 1st Qu.: 0.0000 Median : 0.1000 Mean : 0.1906 3rd Qu.: 0.3750 Max. : 1.3000
Direction_True Min. : 0.0 1st Qu.: 66.0 Median :287.5 Mean :222.5 3rd Qu.:311.8 Max. :359.0
Direction_Mag Min. : 0.0 1st Qu.: 66.0 Median :287.5 Mean :223.3 3rd Qu.:312.0 Max. :359.0
Wind_Speed Min. :0.0000 1st Qu.:0.0000 Median :0.4000 Mean :0.3315 3rd Qu.:0.5000 Max. :1.8000
Heat_Stress_Index Min. :22.80 1st Qu.:25.20 Median :27.20 Mean :27.40 3rd Qu.:29.18 Max. :34.50
Altitude Min. :1675 1st Qu.:1684 Median :1692 Mean :1693 3rd Qu.:1700 Max. :1724
Dew_Point Min. :17.70 1st Qu.:19.40 Median :20.00 Mean :20.16 3rd Qu.:20.80 Max. :24.00
Density_Altitude Min. :2.409 1st Qu.:2.506 Median :2.551 Mean :2.556 3rd Qu.:2.605 Max. :2.723
Wind_Chill Min. :22.20 1st Qu.:24.50 Median :25.80 Mean :25.79 3rd Qu.:27.10 Max. :29.60
Estado_Fenologico_Predominante Min. :717.0 1st Qu.:717.0 Median :718.0 Mean :717.8 3rd Qu.:718.0 Max. :719.0
Frutos_Afectados Min. :0.000000 1st Qu.:0.000000 Median :0.000000 Mean :0.005618 3rd Qu.:0.000000 Max. :1.000000

De la Tabla 1 se destaca que la Temperatura psicrométrica tiene un rango entre 19 y 24.9, con una media de 21.63, indicando una temperatura bastante constante en el área; la presión atmosférica oscila entre 822.4 y 827.4, con media de 825.5; la humedad relativa varía entre 59.5% y 91.6%; el viento cruzado tiene un rango desde 0 hasta 1.5, con una media de 0.1989; la temperatura está en un rango de 22.2 a 29.7, con una media de 25.83; el índice de estrés por calor varía entre 22.8 y 34.5; la altitud se mantiene constante entre 1675 y 1724 metros sobre el nivel del mar. En cuanto a los frutos afectados, el valor medio es muy bajo (0.0056), indicando que en general los frutos no están muy afectados.

leaflet() %>% addTiles() %>% addCircleMarkers(lng = data$Longitude,lat = data$Latitude,radius = 0.1,color = "green")

Figura 1. Localización de árboles

Análisis Geoestadístico

Antes de comenzar con la construcción del variograma, es necesario calcular la matriz de distancias.

summary(dist(data[,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-03

Podemos seleccionar que vaya desde 0 hasta el tercer cuartil (0.0009718) y que la separación sea de 0.00009718 (obtener alrededor de 10 puntos).

En primer lugar, transformaremos los datos en variables georegionalizadas. Para ello, hemos creado una función que recibe el DataFrame y la variable de interés como parámetros. En ella se grafica el variograma y retorna la variable georegionalizada.

convertir_geodata <- function(df, var){
  
  #convertir a variable regionalizada
  data_geo = as.geodata(data, coords.col =3:2 ,data.col = var)
  
  #graficar variograma
  variograma=variog(data_geo,option = "bin",uvec=seq(0,0.0009718,9.718e-05))
  plot(variograma)
  
  variograma_mc=variog.mc.env(data_geo,obj=variograma)
  
  lines(variograma_mc)

  return (data_geo)
}

Temperatura

data_temp = convertir_geodata(data, 'Temperature')
## 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

Figura 2. Variograma de la Temperatura.

Humedad relativa

data_humedad = convertir_geodata(data,'Relative_Humidity')
## 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

Figura 3. Variograma de la humedad relativa.

Presión Barométrica

data_presion = convertir_geodata(data,'Barometric_Pressure')
## 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

Figura 4. Variograma de la Presión.

De la Figura 1, la Figura 2, la Figura 3 se muestra en general que a medida que la distancia aumenta entre los árboles, las condiciones climáticas varían más. No osbastanete, por parte de la temperatura y humedad relativa, no hay un crecimiento constante, puesto que en algunas distancias, la varianza disminuye. A partir de este momento trabajaremos con la presión, ya que se observa una autocorrelación espacial.

Modelo teorico

Ajustamos el modelo del semivariograma.

ini.vals = expand.grid(seq(8,12,l=10), seq(6e-04,8e-04,l=10))

variograma = variog(data_presion,option = "bin",uvec=seq(0,0.0009718,9.718e-05))
## variog: computing omnidirectional variogram
plot(variograma)

#modelo exponencial
presion_exp = variofit(variograma, 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 1008603.27778179
#modelo gausiano
presion_gauss =variofit(variograma, 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 789254.080595493
#modelo esférico
presion_sph =variofit(variograma, 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 "8"     "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 3475061.90029651
lines(presion_exp,col="blue")
lines(presion_gauss,col="red")
lines(presion_sph,col="purple")

Figura 5. Modelos teóricos para la presión.

De acuerdo con la Figura 5 y el valor SCE, el modelo gaussiano es el que mejor se ajusta a los puntos del semivariograma muestral.

Predicción Espacial Kriging

Para generar el mapa, sacamos el minimo y máximo de la longitud y latitud.

c(min(data[,3]),max(data[,2]))
## [1] -76.711799   2.393634

Ahora vamos a realizar la predicción espacial.

data_grid <- expand.grid(lon=seq(-76.7118, -76.71022, l=100),
                                 lat=seq(2.392101, 2.393634, l=100))

plot(data_grid)

points(data[,3:2],col="green")

Figura 6. Terreno de los árboles

data_presion_pred <- krige.conv(
  data_presion, loc=data_grid, krige = krige.control(nugget=0,trend.d="cte", 
      trend.l="cte",cov.pars=c(sigmasq=12.7080, phi=0.0007 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
image(data_presion_pred, xlab="Longitud", ylab="Latitud")
par(new=TRUE)
contour(data_presion_pred,  xlab="Longitud", ylab="Latitud", drawlabels=TRUE)

Figura 7. Predicción Espacial

par(mfrow=c(1,2))

image(data_presion_pred, 
      val=sqrt(data_presion_pred$krige.var), xlab="Longitud", ylab="Latitud")

par(new=TRUE)

contour(data_presion_pred, xlab="Longitud", 
        val=sqrt(data_presion_pred$krige.var), ylab="Latitud", drawlabels=TRUE)

Figura 8. Error de la Predicción Espacial

Formato Raster

raster_pred <- rasterFromXYZ(cbind(data_grid, data_presion_pred$predict))

plot(raster_pred,
     xlab="Longitud", ylab="Latitud")

Figura 9 Predicción Espacial en formato Raster

# Creando mapa raster del error a partir de los datos
raster_pred_error <- rasterFromXYZ(cbind(data_grid, val=sqrt(data_presion_pred$krige.var)))

plot(raster_pred_error,
     xlab="Longitud", ylab="Latitud")

Figura 10 Error de la Predicción Espacial en formato Raster

presion_valid=xvalid(geodata = data_presion, 
                           model = presion_gauss)
## xvalid: number of data locations       = 534
## xvalid: number of validation locations = 534
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 
## xvalid: end of cross-validation
MAE=mean(abs(presion_valid$error))
print(MAE)
## [1] 0.6513166

El modelo gaussiano presenta un MAE de 0.65.

Conclusión

Aunque la temperatura y la humedad relativa mostraron autocorrelación espacial a ciertas distancias, sus variogramas evidenciaron comportamientos menos consistentes, como disminuciones en la varianza a medida que la distancia aumentaba. Por otro lado, la presión barométrica presenta un patrón de correlación espacial más definido, lo cual puede indicar que la presión está más influenciada por factores espaciales, mientras que las otras variables podrían estar afectadas por otros factores externos. Finalmente, el modelo gaussiano fue el que presentó el mejor ajuste para describir la presión barométrica, con un error absoluto medio (MAE) de 0.65 en la validación cruzada.