## ----Loading_library,include=FALSE,echo=FALSE---------------------------------
knitr::opts_chunk$set(echo = TRUE, fig.align="center")
library(RGeostats)
constant.define("asp",1)


## -----------------------------------------------------------------------------
nx = 50
x0 = 0
dx = 1
dbTarget = db.create(nx=rep(nx+1,2),x0=rep(x0,2),dx=rep(dx,2))
dbTarget = db.locate(dbTarget,2:3,"x")

nObs = 5
set.seed(99)
coordObs = cbind(sample(seq(nx),nObs,replace=F),sample(seq(nx),nObs,replace=F))
valObs = rnorm(nObs)
dbObs = db.create(cbind(coordObs,valObs))
dbObs = db.locate(dbObs,2:3,"x")
dbObs = db.locate(dbObs,4,"z")
plot(dbObs,pos.legend=1,title="Observations")


## -----------------------------------------------------------------------------
exponentialCov <- function(h,range,sill){
  sill*exp(-h/range)
}
a = 10
C = 1


## -----------------------------------------------------------------------------
distObs = sqrt(outer(coordObs[,1],coordObs[,1],"-")^2 + outer(coordObs[,2],coordObs[,2],"-")^2)
cxx = exponentialCov(distObs,range=a,sill=C)

distTar = sqrt(outer(coordObs[,1],dbTarget[,"x1"],"-")^2 + outer(coordObs[,2],dbTarget[,"x2"],"-")^2)
cx0 = exponentialCov(distTar,range=a,sill=C)

c00 = exponentialCov(rep(0,dbTarget$nech),range=a,sill=C)


## -----------------------------------------------------------------------------
def_cov <- function(cxx,cx0,c00) {
  function(dist,db1=NA,iech1=NA,db2=NA,iech2=NA,incr=NA,x1=NA,x2=NA,...){
    result = 0
    if (db1 == 1 && db2 == 1)
        result = cxx[iech1,iech2] 
    else if (db1 == 1 && db2 == 2)
        result = cx0[iech1,iech2]
    else if (db1 == 2 && db2 == 1) 
        result = cx0[iech2,iech1]
    else if (db1 == 2 && db2 == 2)
        result = c00[iech1]
    else
        cat("db1=",db1," db2=",db2,"\n")
    result
  }
}

my_cov_func <- def_cov(cxx=cxx,cx0=cx0,c00=c00)


## -----------------------------------------------------------------------------
unique.neigh = neigh.create(type=0)
dbTarget = kriging(dbin=dbObs,dbout=dbTarget,model=my_cov_func,neigh=unique.neigh,uc="1")
plot(dbTarget,pos.legend=1,title="Kriging estimate")
plot(dbTarget,name="*.stdev",pos.legend=1,title="Standard deviation of kriging error")


## -----------------------------------------------------------------------------
model = model.create(vartype="exponential",scale=a,sill=C)
dbTarget = kriging(dbin=dbObs,dbout=dbTarget,model=model,neigh=unique.neigh,uc="1",radix="OK")


## -----------------------------------------------------------------------------
dummy = correlation(dbTarget,name1="Kriging.*.estim",name2="OK.*.estim",
                    flag.iso=T,flag.diag=T,
                    title="Estimations")
dummy = correlation(dbTarget,name1="Kriging.*.stdev",name2="OK.*.stdev",
                    flag.iso=T,flag.diag=T,
                    title="Standard deviations")