Análisis geoespacial

2. Geoestadística

El objetivo del análisis estadística evaluar si existe una correlación espacial en las mediciones de temperatura, usando técnicas geoestadísticas en R. Inicialmente se crea el objeto geodata el cual es fundamental para el análisis.

2.1. Crear geodata

El objeto se crea a partir de la base cargada anteriormente indicando las columnas de longitud y latitud de las mediciones y la columna de la variable de interés (Temperatura) que en este caso es la 5.

geo_temp=as.geodata(geo_datos,coords.col = 3:2,data.col = 5)
plot(geo_temp)

En este reporte inicial se ve que las temperaturas se encuentran entre los 19 y 25 °C y que hay ciertos lugares donde se concentran en su mayoria temperaturas altas o bajas.

2.2. Variograma

Primero se hace el análisis de las distancias entre los puntos, para evaluar el rango que tendrá la secuencia en la creación del variograma. Los rangos de distancia estarán entre 0 y 0.001, separados por 0.00002.

options(scipen=999)
summary(dist(geo_datos[,3:2]))

##       Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
## 0.00001712 0.00040512 0.00064078 0.00068267 0.00091776 0.00195913

variograma=variog(geo_temp,option = "bin",uvec=seq(0,0.001,0.00002))

## variog: computing omnidirectional variogram

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

plot(variograma)
lines(variograma_mc)

El semivariograma indica que posiblemente existe una correlación espacial ya que para distancias más pequeñas la semivarianza es más baja. Para las distancias más pequeñas la semivarianza inicia en 0.7 aproximadamente y se estabiliza en 1.3 a una distancia de 0.00025 aproximandamente.

2.3. Ajuste del modelo teórico

Se ajustan tres modelos teóricos: - Exponencial - Gaussiano - Esférico

ini.vals = expand.grid(seq(1.2,1.5,l=10), seq(0.0003,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: 6624.4190268492

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.47"  "0"   "0"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 11320.824524992

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.33"  "0"   "0"   "0.5"
## status        "est"   "est" "fix" "fix"
## loss value: 1216.01931350849

plot(variograma)

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

model_mco_exp

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

Según el gráfico y la suma de cuadrados el modelo que mejor se ajuste es el exponencial. Se estima un valor de sigmasq = 1.3138 y phi = 0.0001.

2.4. Interpolación (Predicción espacial)

Para la interpolación se crea la malla con la extensión de las mediciones, la cual sirve como entrada del kriggin con los datos y los parámetros estimados en el modelo exponencial.

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

## [1] -76.711799 -76.710215   2.392101   2.393634

geodatos_grid=expand.grid( lon=seq(-76.71022,-76.7118,l=100),lat=seq(2.392101 ,2.393634 ,l=100))
plot(geodatos_grid)
points(geo_datos[,3:2],col="red")

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

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

El resultado de la predicción espacial es el siguiente.

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

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

Y la desviación estandar espacial:

par(mfrow=c(1,2))

contour(geodatos_ko,main="kriging StDv Predict",val=sqrt(geodatos_ko$krige.var),  drawlabels=TRUE)
image(geodatos_ko, main="kriging StDv Predicted",val=sqrt(geodatos_ko$krige.var), xlab="East", ylab="North")

2.5. Transformar imagen a raster

Finalmente se transforma el resultado anterior en un objeto raster.

pred=cbind(geodatos_grid,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_grid,geodatos_ko$predict))
levelplot(temp_predict,par.settings =BuRdTheme)

Con este resultado ya se tiene un conjunto de predicciones de temperatura en lugares donde esta no se midió.