require(geoR)
## Loading required package: geoR
## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.8-1 (built on 2020-02-08) is now loaded
## --------------------------------------------------------------
library(readxl)
library(readxl)
require(rasterVis)
## Loading required package: rasterVis
## Loading required package: raster
## Loading required package: sp
## Loading required package: lattice
## Loading required package: latticeExtra
require(raster)
library(readxl)
datos <- read_excel("C:/Users/linam/Desktop/Doctorado-20201031T024248Z-001/Doctorado/Semestre V/Analisis Espacial de Datos/base_cienaga.xls")
head(datos)
## # A tibble: 6 x 6
## Este Norte prof temp sali oxig
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 976952 1706444 1.75 26.4 29.0 6.34
## 2 970883. 1704880. 1.4 30.3 13.0 9.42
## 3 972704. 1704827. 1.1 30.4 19.2 7.88
## 4 974718 1704874 1.5 28.3 35.0 5.67
## 5 976538 1704874 0.5 27 34.0 5.49
## 6 955344. 1703160. 1.62 26 18.0 5.07
plot(datos[,1:2])
##names(datos)=c("Este","Norte","Temp")
points(datos[,1:2], col="red", pch=16)
geodatos=as.geodata(datos,coords.col=1:2,data.col=4)
plot(geodatos)
##Semivariograma experimental
summary(dist(geodatos$coords))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1373 6994 11204 11531 15606 31924
summary(dist(geodatos$coords[1:5,]))
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1624 1870 3282 3415 4368 6268
##geodatos$data=sample(geodata$data)
variograma=variog(geodatos,option = "bin",uvec=seq(0,20000,1500))
## variog: computing omnidirectional variogram
plot(variograma,pch=16)
##lines(variograma_aleatorio,col="red")
geodatos.env=variog.mc.env(geodatos,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, pch=16, main=names(datos)[4], envelope=geodatos.env)
names(datos)
## [1] "Este" "Norte" "prof" "temp" "sali" "oxig"
##Ajuste del Modelo de Semivarianza
ini.vals=expand.grid(seq(5,7,l=10), seq(10000,5000,l=10))
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 "6.33" "7777.78" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 353.860197210756
model_mco_exp=variofit(variograma,ini=ini.vals,cov.model="exponential",wei="npair", min="optim", nugget=0)
## 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 "7" "8888.89" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 2484.41516301323
model_mco_sphe=variofit(variograma,ini=ini.vals,cov.model="spheric",wei="npair", min="optim", nugget=0)
## 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 "5.22" "10000" "0" "0.5"
## status "est" "est" "est" "fix"
## loss value: 6345.73259045248
lines(model_mco_gaus, col="red")
lines(model_mco_exp, col="blue")
lines( model_mco_sphe,col="purple")
##Predicciòn Espacial Kriging
geodatos_gird=expand.grid(Este=seq(952023,979916,l=100),Norte=seq(1677586, 1708649,l=100))
plot(geodatos_gird)
geodatos_ko=krige.conv(geodatos,loc=geodatos_gird,krige=krige.control(nugget=0,trend.d = "cte", trend.l= "cte", cov.pars =c(sigmasq=6.2810, phi=7971.9857 )))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow=c(1,1))
image(geodatos_ko,main="kriging predict", xlab="East", ylab="North")
contour(geodatos_ko,main="Kriging Predict", add=TRUE, drawlabels=TRUE)
image(geodatos_ko,main="kriging StDv predict", val=sqrt(geodatos_ko$krige.var),xlab="East", ylab="North")
contour(geodatos_ko,main="kriging StDv predict", val=sqrt(geodatos_ko$krige.var),xlab="East", ylab="North")
require(rasterVis)
require(raster)
pred=cbind(geodatos_gird,geodatos_ko$predict)
temp_predict=rasterFromXYZ(cbind(geodatos_gird,geodatos_ko$predict))
plot(temp_predict)
temp_error=rasterFromXYZ(cbind(geodatos_gird,geodatos_ko$krige.var))
plot(temp_error)
