Unidad 3 Temperatura Aguacate
Diana Socha
2024-06-14
Utilizando la geoestadística, es necesario llevar a cabo un análisis exhaustivo del clima y de las condiciones geográficas óptimas para la siembra y el cultivo de aguacates.
Primero que todo cargaremos los datos completos del aguacate desde la fuente de datos correspondiente.
Aguacates <- read_excel("Datos_Completos_Aguacate.xlsx")
summary(Aguacates)## id_arbol Latitude Longitude FORMATTED_DATE_TIME
## Length:20271 Min. :2.316 Min. :-76.71 Length:20271
## Class :character 1st Qu.:2.371 1st Qu.:-76.71 Class :character
## Mode :character Median :2.373 Median :-76.61 Mode :character
## Mean :2.367 Mean :-76.66
## 3rd Qu.:2.392 3rd Qu.:-76.61
## Max. :2.394 Max. :-76.61
##
## Psychro_Wet_Bulb_Temperature Station_Pressure Relative_Humidity
## Min. :11.10 Min. :798.7 Min. : 21.60
## 1st Qu.:16.70 1st Qu.:804.5 1st Qu.: 50.70
## Median :18.10 Median :807.4 Median : 69.30
## Mean :18.09 Mean :811.9 Mean : 67.78
## 3rd Qu.:19.50 3rd Qu.:821.8 3rd Qu.: 83.20
## Max. :27.90 Max. :830.1 Max. :100.00
##
## Crosswind Temperature Barometric_Pressure Headwind
## Min. :0.0000 Min. :12.90 Min. :798.7 Min. :-5.30000
## 1st Qu.:0.0000 1st Qu.:20.10 1st Qu.:804.5 1st Qu.: 0.00000
## Median :0.0000 Median :22.70 Median :807.4 Median : 0.00000
## Mean :0.2124 Mean :22.91 Mean :811.8 Mean : 0.06135
## 3rd Qu.:0.4000 3rd Qu.:25.70 3rd Qu.:821.7 3rd Qu.: 0.20000
## Max. :6.6000 Max. :37.80 Max. :830.0 Max. : 4.60000
##
## Direction_True Direction_Mag Wind_Speed Heat_Stress_Index
## Min. : 0.0 Min. : 0.0 Min. :0.0000 Min. :13.30
## 1st Qu.:124.0 1st Qu.:125.0 1st Qu.:0.0000 1st Qu.:20.30
## Median :237.0 Median :238.0 Median :0.3000 Median :22.80
## Mean :208.1 Mean :208.1 Mean :0.3229 Mean :23.14
## 3rd Qu.:294.0 3rd Qu.:294.0 3rd Qu.:0.5000 3rd Qu.:25.60
## Max. :359.0 Max. :359.0 Max. :7.2000 Max. :48.90
##
## Altitude Dew_Point Density_Altitude Wind_Chill
## Min. :1648 Min. : 5.40 Min. :2.105 Min. :12.80
## 1st Qu.:1730 1st Qu.:14.10 1st Qu.:2.484 1st Qu.:20.10
## Median :1872 Median :16.50 Median :2.596 Median :22.70
## Mean :1828 Mean :15.95 Mean :2.601 Mean :22.86
## 3rd Qu.:1901 3rd Qu.:18.10 3rd Qu.:2.707 3rd Qu.:25.60
## Max. :1960 Max. :27.50 Max. :3.152 Max. :37.80
##
## Estado_Fenologico_Predominante Frutos_Afectados
## Min. : 0 Min. : 0.0000
## 1st Qu.:715 1st Qu.: 0.0000
## Median :716 Median : 0.0000
## Mean :715 Mean : 0.1408
## 3rd Qu.:717 3rd Qu.: 0.0000
## Max. :719 Max. :60.0000
## NA's :22top5_data <- head(Aguacates, 5)
kable(top5_data)| id_arbol | Latitude | Longitude | FORMATTED_DATE_TIME | Psychro_Wet_Bulb_Temperature | Station_Pressure | Relative_Humidity | Crosswind | Temperature | Barometric_Pressure | Headwind | Direction_True | Direction_Mag | Wind_Speed | Heat_Stress_Index | Altitude | Dew_Point | Density_Altitude | Wind_Chill | Estado_Fenologico_Predominante | Frutos_Afectados |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 2.381290 | -76.61264 | 21/08/2019 9:22:57 a, m, | 14.8 | 805.1 | 33.6 | 0.2 | 25.7 | 805.0 | 0.7 | 166 | 165 | 0.8 | 24.1 | 1896 | 8.6 | 2.743 | 25.7 | 715 | 0 |
| 2 | 2.381320 | -76.61259 | 21/08/2019 9:27:13 a, m, | 11.6 | 805.2 | 36.8 | 3.6 | 20.8 | 805.2 | 3.5 | 314 | 313 | 5.1 | 19.5 | 1895 | 5.5 | 2.570 | 20.8 | 715 | 3 |
| 3 | 2.381353 | -76.61254 | 21/08/2019 9:36:36 a, m, | 12.9 | 805.8 | 31.5 | 0.4 | 23.7 | 805.7 | 0.7 | 332 | 331 | 0.8 | 22.0 | 1889 | 5.8 | 2.661 | 23.7 | 715 | 0 |
| 4 | 2.381374 | -76.61261 | 21/08/2019 9:38:02 a, m, | 14.1 | 805.7 | 33.2 | 0.6 | 25.0 | 805.7 | 0.7 | 139 | 139 | 0.9 | 23.2 | 1890 | 7.7 | 2.707 | 24.9 | 715 | 0 |
| 5 | 2.381335 | -76.61266 | 21/08/2019 9:39:38 a, m, | 14.3 | 805.2 | 34.3 | 0.4 | 25.0 | 805.2 | 0.4 | 129 | 128 | 0.6 | 23.3 | 1894 | 8.1 | 2.714 | 24.9 | 715 | 1 |
Tras cargar la información de la fuente de datos, es aconsejable mapear la ubicación de las plantaciones de aguacates.
leaflet() %>% addTiles() %>% addCircleMarkers(lng = Aguacates$Longitude,lat = Aguacates$Latitude,radius = 0.2,color = "#12A17A")Geodata
Ahora vamos a regionalizar la variable
geod_aguacate=as.geodata(Aguacates,coords.col = 3:2,data.col = 9)## as.geodata: 18586 replicated data locations found.
## Consider using jitterDupCoords() for jittering replicated locations.
## WARNING: there are data at coincident or very closed locations, some of the geoR's functions may not work.
## Use function dup.coords() to locate duplicated coordinates.
## Consider using jitterDupCoords() for jittering replicated locationsplot(geod_aguacate)
Semivariograma
summary(dist(geod_aguacate$coords))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000000 0.001756 0.077784 0.063476 0.106252 0.109962Variograma
variograma=variog(geod_aguacate,option = "bin",uvec=seq(0,0.001,0.00002))## variog: computing omnidirectional variogram
## variog: co-locatted data found, adding one bin at the origindatos.env=variog.mc.env(geod_aguacate,obj=variograma)## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelopsSe realiza la gráfica del Variograma
plot(variograma, main="Grafico variograma")
lines(datos.env)
Modelo ajustado
Para realizar el ajuste del modelo de semivarianza se tiene en cuenta el uso de modelos como el exponencial, esferico, y gaussiano.
primero el modelo exponencial
ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,0.0008,l=10))
model_mco_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 "1.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 5351625789.89722Segundo el modelo esférico.
model_mco_spe=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 "1.5" "0" "0" "0.5"
## status "est" "est" "fix" "fix"
## loss value: 5297606635.41097Tercero el modelo Gaussiano.
model_mco_gaus=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 "1.5" "0" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 5325691774.74162Gráfica de los modelos
plot(variograma)
lines(model_mco_exp,col="lightblue")
lines(model_mco_gaus,col="pink")
lines(model_mco_spe,col="lightyellow")
Resultados
model_mco_exp## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 12.5756 0.9163 0.0013
## Practical Range with cor=0.05 for asymptotic range: 0.004020816
##
## variofit: minimised weighted sum of squares = 1234962Por tanto se tiene que el ajuste del modelo de semivarianza se realizará mediante el modelo exponencial debido a que es el que muestra la menor suma de cuadrados.
model_mco_gaus## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
## tausq sigmasq phi
## 12.8480 0.0011 0.0013
## Practical Range with cor=0.05 for asymptotic range: 0.002165812
##
## variofit: minimised weighted sum of squares = 1484166Grafico Variograma mediante la función exponencial
model_mco_spe## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq = 0
## parameter estimates:
## sigmasq phi
## 12.8448 0.0000
## Practical Range with cor=0.05 for asymptotic range: 0
##
## variofit: minimised weighted sum of squares = 23884586Predicción espacial
En este proceso se estimara los valores de la variable temperatura en ubicaciones no muestreadas dentro de un área geográfica, basado en datos de ubicaciones conocidas.
Primero se creara una vecindad de puntos máximos y mínimos de la longitud y la latitud
c(min(Aguacates[,3]),
max(Aguacates[,3]),
min(Aguacates[,2]),
max(Aguacates[,2]))## [1] -76.711799 -76.606710 2.316405 2.393634y se creará la malla.
geodatos_grid=expand.grid( lon=seq(-76.606710,-76.711799,l=100),lat=seq(2.316405 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geod_aguacate$coords,col="pink")
Predicción del modelo
x.range <- range(geod_aguacate$coords[, 1])
y.range <- range(geod_aguacate$coords[, 2])
grid_resolution <- 1
geodatos_grid <- expand.grid(
x = seq(x.range[1], x.range[2], by = grid_resolution),
y = seq(y.range[1], y.range[2], by = grid_resolution)
)
cat("Tamaño del grid:", nrow(geodatos_grid), "puntos\n")## Tamaño del grid: 1 puntoshead(geodatos_grid)## x y
## 1 -76.7118 2.316405