Se cargan las librerías y los datos

Previamente se filtró la base de datos para las observaciones de 2020, esto con el fin de acotar el análisis a un solo lote cerca de la ciudad de Popayán en la vereda La Chorrera.

leaflet() %>% addTiles() %>% addCircleMarkers(lng = datos$Longitude,lat = datos$Latitude,radius = 0.5,color = "yellow")

Variable regionalizada

Convertimos nuestra base en una variable regionalizada con las columnas de longitud, latitud y los valores de temperatura. Esta oscila entre 22 y 30 grados, con una media en 26. En colores podemos observar como los árboles se agrupan entre la misma temperatura.

temp=as.geodata(datos,coords.col = 3:2,data.col = 9)
plot(temp)

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

Variograma

Aunque el semiovariograma que observamos nos puede indicar que la variable tiene el comportamiento esperado, ajustaremos diversos modelos para decidir cual de todos es el que se ajusta mejor a nuestros datos y a lo que se espera de un semiovariograma.

variograma=variog(temp,option = "bin",uvec=seq(0,0.000917,9.178e-05))
## variog: computing omnidirectional variogram
plot(variograma)
variograma_temp=variog.mc.env(temp,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 envelops
lines(variograma_temp)

show.vgms(par.strip.text = list(cex = 0.75))

Modelos ajustados de semivarianza

Ajustando los modelos elegidos, podemos observar como el modelo Gaussiano y el Exponencial son los que mejor se ajustan a nuestros datos, sin embargo el mejor ajuste se alcanza con el modelo exponencial.

ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0001,0.0008,l=10))
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: 240018.669934895
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: 231440.514329651
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: 226841.340327613
plot(variograma)

lines(exp ,col="green")
lines(gaus ,col="red")
lines(spe ,col="yellow")

La intuición visual se confirma con el resumen de los modelos, siendo el exponencial el que presenta menor suma de cuadrados de los errores.

exp
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7849  2.2672  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0003390645
## 
## variofit: minimised weighted sum of squares = 3907.716
gaus
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7594  2.2508  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
## 
## variofit: minimised weighted sum of squares = 5241.893
spe
## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: spherical
## fixed value for tausq =  0 
## parameter estimates:
## sigmasq     phi 
##  2.9538  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 12203.46

Malla

Con los mínimos de las coordenadas de latitud y longitud creamos una malla para sobreponer la predicción que lograremos con el variograma.

c(min(datos[,3]),
  max(datos[,3]),
  min(datos[,2]),
  max(datos[,2]))
## [1] -76.711799 -76.710215   2.392101   2.393634
geodatos_grid=expand.grid( lon=seq(-76.711799,-76.710215,l=100),lat=seq(2.392101 ,2.393634 ,l=200)) 
plot(geodatos_grid)
points(temp$coords,col="purple")

Predicción

Se procede a la construcción del semivariograma teórico con sigmasq=2.2672, phi=0.0001, correspondientes el modelo exponencial antes definido.

prediccion=krige.conv(temp, loc=geodatos_grid,
                       krige= krige.control(nugget=0,trend.d="cte", 
                                            trend.l="cte",cov.pars=c(sigmasq=2.2672, phi=0.0001 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,2))
contour(prediccion,main="kriging Predict", drawlabels=TRUE)
image(prediccion, main="kriging Predict", xlab="East", ylab="North")

Mejoras Gráficas

Plot simple con colores

Podemos utilizar una gama de colores para visualizar nuestros resultados.

pred=cbind(geodatos_grid, prediccion$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid, prediccion$predict))
colores <- colorRampPalette(c("blue", "green", "yellow", "red"))
plot(temp_predict, col.regions= colores, main = "Predicción de Temperatura con kriging")
## Warning in plot.window(...): "col.regions" is not a graphical parameter
## Warning in plot.xy(xy, type, ...): "col.regions" is not a graphical parameter
## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter

## Warning in axis(side = side, at = at, labels = labels, ...): "col.regions" is
## not a graphical parameter
## Warning in box(...): "col.regions" is not a graphical parameter
## Warning in title(...): "col.regions" is not a graphical parameter

Levelplot

Con el Level Plot podemos generar un gráfico con la gráfica de distribución de los datos adjunta. Esto nos puede ayudar a interpretar que entre menor sea la coordenada X y Y, mayor es la temperatura.

levelplot(temp_predict,par.settings =BuRdTheme, col.regions= colores, main = "Mapa de calor del kriging",
          xlab = "Coordenada X",
          ylab = "Coordenada Y")

Contornos

A su vez podemos sobreponer los valores de la predicción y los niveles en un Mapa de contornos que pueden ser un predictor de la altitud.

r <- rasterFromXYZ(pred[, c( "lon", "lat", "prediccion$predict")])

contourplot(r, 
            main = list(label = "Mapa de Contornos del Kriging", cex = 1.5, font = 2),
            xlab = list(label = "Coordenada X", cex = 1.2),
            ylab = list(label = "Coordenada Y", cex = 1.2),
            col.regions= colores,
            labels = list(cex = 1.0),
            region = TRUE)

Mapa interactivo

Con el mapa interactivo podemos observar la estimación puntual para cada punto y con esto decidir con mayor precisión el lugar donde debemos plantar un árbol tomando en cuenta la temperatura esperada.

p <- ggplot(pred, aes(x = lon, y = lat, color = pred$`prediccion$predict`)) +
  geom_point() +
  scale_color_viridis() +
  labs(title = "Mapa de Puntos interactivo del Kriging",
       x = "Coordenada X",
       y = "Coordenada Y",
       color = colores) +
  theme_minimal()


p_interactive <- ggplotly(p)
## Warning: Use of `` pred$`prediccion$predict` `` is discouraged.
## ℹ Use `prediccion$predict` instead.
p_interactive

Mapa 3D

Con el mapa interactivo 3D se puede observar en 3 coordenadas los valores estimados.

z <- matrix(pred$`prediccion$predict`, nrow = length(unique(pred$lon)), ncol = length(unique(pred$lat)))
x <- unique(pred$lon)
y <- unique(pred$lat)


plot_ly(z = ~z, x = ~x, y = ~y, type = "surface", colors = viridis(100)) %>%
  layout(title = "Diagrama de Superficie 3D del Kriging",
         scene = list(
           xaxis = list(title = "Coordenada X"),
           yaxis = list(title = "Coordenada Y"),
           zaxis = list(title = "Valor Estimado")
         ))

Conclusiones

  1. La geomorfología del lote nos hace notar que en el suroeste y en el centro “elevado” se presentan las temperaturas más altas.
  2. Dado que la temperatura ideal del aguacate esta entre 17 a 24 grados. Las zonas centrales elevadas y la zona suroeste no es óptima para su rendimiento.
  3. La mayor parte del terreno es apto para el cultivo de aguacate, siendo las zonas calientes ideales para crear infraestructura de sostenimiento del cultivo.