Vario.calc with measurement error

How to use RGeostats package step-by-step.

Vario.calc with measurement error

Another tutorial for vario.calc:

First example

First, we simulate the data without measurement error and we compute the experimental variogram as a reference.
Code: Select all
`set.seed(11111)db=db.create(x1=sort(runif(400)))model=model.create(melem.name(3),range=.2,ndim=1)db=simtub(dbout=db,model=model,nbtuba=1000)vario_true=vario.calc(db,lag=0.03,nlag=15)`

Then we add some measurement error on 90% of the points, half with variance 1 and half with variance 4:

Code: Select all
`error_sd=sample(c(rep(0,10),rep(1,90),rep(3,300)))error=rnorm(400,0,error_sd)db[,3]=db[,3]+errordb=db.add(db,error_sd^2,locnames="v1")`

We can compute the experimental variogram with only the points without measurement error:

Code: Select all
`db=db.add(db,(error_sd==0),locnames="sel")vario_only=vario.calc(db,lag=0.03,nlag=15)db=db.locate(db,5)`

We can display the data:

Code: Select all
`plot(db[,2],db[,3],col=1+(error_sd==1)+2*(error_sd==3),pch=16,cex=.8,xlab="x",ylab="z",ylim=c(-12,12))lines(db[,2],db[,3]-error)legend(0.4,-7,legend=c("True signal","Data without error","Data with error variance 1","Data with error variance 9"),lty=c(1,NA,NA,NA),pch=c(NA,16,16,16),col=c(1,1,2,3),box.lty=0)`

We then compute:
The usual experimental variogram
The experimental variogram with bias correction
The two iterative variograms

Code: Select all
`vario0=vario.calc(db,lag=0.03,nlag=15)vario1=vario.calc(db,lag=0.03,nlag=15,verr.mode=1)vario2=vario.calc(db,lag=0.03,nlag=15,verr.mode=2)vario3=vario.calc(db,lag=0.03,nlag=15,verr.mode=3)`

Then we can display and compare the results:

Code: Select all
`plot(vario0,col=2,varline=F,ylim=c(-.5,8),xlab="h",ylab=expression(gamma))plot(vario_only,col=6,lty=5,add=T,varline=F)plot(vario1,col=3,lty=2,add=T,varline=F)plot(vario2,col=4,lty=3,add=T,varline=F)plot(vario3,col=5,lty=4,add=T,varline=F)plot(vario_true,lwd=2,add=T,varline=F)legend(0.1,5,legend=c("Reference","With points with no measurement error",expression(gamma(h)),expression(gamma[a]^"+"~(h)),expression(gamma[b]^"+"~(h)),expression(gamma[c](h))),col=c(1,6,2,3,4,5),lwd=c(2,1,1,1,1,1,1),lty=c(1,5,1,2,3,4),box.col=0)`

And a zoom of the previous figure:

Code: Select all
`plot(vario_only,col=6,lty=5,varline=F,xlab="h",ylab=expression(gamma),ylim=c(0,2.5))plot(vario1,col=3,lty=2,add=T,varline=F)plot(vario2,col=4,lty=3,add=T,varline=F)plot(vario3,col=5,lty=4,add=T,varline=F)plot(vario_true,lwd=2,add=T,varline=F)legend(0.1,2.4,legend=c("Reference","With points with no measurement error",expression(gamma[a]^"+"~(h)),expression(gamma[b]^"+"~(h)),expression(gamma[c](h))),col=c(1,6,3,4,5),lwd=c(2,1,1,1,1,1),lty=c(1,5,2,3,4),box.col=0)`

Nicolas Desassis

Posts: 196
Joined: Thu Sep 20, 2012 4:22 pm