Unidad 3 - Análisis de información geográfica

El objetivo del presente Script es realizar una exploración de datos espaciales utilizando geoestadística, sobre la producción de aguacate en un conjunto de arboles determinado en la región del Cauca, y así identificar si existe una correlación espacial de la variable temperatura.

#Cargar la información de los aguacates
library(readxl)
Datos_aguacates <- read_excel("C:/Users/PACHO/Documents/Maestria/SIG/Datos_Completos_Aguacate.xlsx")
head(Datos_aguacates)

## # A tibble: 6 × 21
##   id_arbol Latitude Longitude FORMATTED_DATE_TIME        Psychro_Wet_Bulb_Temp…¹
##      <dbl>    <dbl>     <dbl> <chr>                                        <dbl>
## 1        1     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22  
## 2        2     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.4
## 3        3     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.8
## 4        4     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.8
## 5        5     2.39     -76.7 01/10/2020  10:11:12 a, m,                    22.6
## 6        6     2.39     -76.7 01/10/2020  10:11:12 a, m,                    21.5
## # ℹ abbreviated name: ¹Psychro_Wet_Bulb_Temperature
## # ℹ 16 more variables: Station_Pressure <dbl>, Relative_Humidity <dbl>,
## #   Crosswind <dbl>, Temperature <dbl>, Barometric_Pressure <dbl>,
## #   Headwind <dbl>, Direction_True <dbl>, Direction_Mag <dbl>,
## #   Wind_Speed <dbl>, Heat_Stress_Index <dbl>, Altitude <dbl>, Dew_Point <dbl>,
## #   Density_Altitude <dbl>, Wind_Chill <dbl>,
## #   Estado_Fenologico_Predominante <dbl>, Frutos_Afectados <dbl>

#Graficar la unicación de los arboles
require(leaflet)

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

Crear el objeto Geodata

#Se convierten los datos en una variable regionalizada, que son las coordenadas mas la variable
require(geoR)
geo_aguacate=as.geodata(Datos_aguacates,coords.col = 3:2,data.col = 9)
plot(geo_aguacate)

De las imagenes anteriores, se puede observar que la temperatura oscila entre los 22 y 30 grados, y también, permite observar como se encuentren ahrupados las arboles con temperaturas similares.

Crear el semivariograma

summary(dist(geo_aguacate$coords))

##      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=variog(geo_aguacate,option = "bin",uvec=seq(0,0.001,0.00002))

## variog: computing omnidirectional variogram

datos.env=variog.mc.env(geo_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 envelops

plot(variograma)
lines(datos.env)

Ajuste del modelo de semivarianza

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: 296807.231725087

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: 288205.86784414

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: 283743.808438713

plot(variograma)

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

#Resultado modelo exponencial
model_mco_exp

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7969  2.2860  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.000232806
## 
## variofit: minimised weighted sum of squares = 9158.212

#Resultado del modelo Gaussiano
model_mco_gaus

## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: gaussian
## parameter estimates:
##   tausq sigmasq     phi 
##  0.7902  2.2842  0.0001 
## Practical Range with cor=0.05 for asymptotic range: 0.0002330932
## 
## variofit: minimised weighted sum of squares = 8947.68

#Resultado del modelo esférico
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 
##  3.0186  0.0000 
## Practical Range with cor=0.05 for asymptotic range: 0
## 
## variofit: minimised weighted sum of squares = 17731.76

De acuerdo con los tres modelos observados, el que presenta menor suma de cuadrados de los errores es el modelo esférico con un valor de 17731.

Predicción espacial

#Se obtiene los limites de la malla a crear

c(min(Datos_aguacates[,3]),
max(Datos_aguacates[,3]),
min(Datos_aguacates[,2]),
max(Datos_aguacates[,2]))

## [1] -76.711799 -76.710215   2.392101   2.393634

#Se crea la malla
geodatos_grid=expand.grid( lon=seq(-76.710215,-76.711799,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_aguacate$coords,col="red")

#Predicción del modelo utilizando el modelo esférico

geodatos_ko=krige.conv(geo_aguacate, loc=geodatos_grid,
                       krige= krige.control(nugget=0,trend.d="cte", 
                                            trend.l="cte",cov.pars=c(sigmasq=3.0186, phi=0.0001 )))

## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood

par(mfrow=c(1,2))
contour(geodatos_ko,main="kriging Predict", drawlabels=TRUE)

image(geodatos_ko, main="kriging Predict", xlab="East", ylab="North")

Transformar a imagen raster