Vario.calc with measurement error

How to use RGeostats package step-by-step.
Direct access to Tutorial index and Frequently Asked Questions

Vario.calc with measurement error

Postby Nicolas Desassis » Sat Feb 08, 2014 3:00 pm

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]+error
db=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)




Image

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)



Image


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)



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

Return to Tutorials and FAQ

Who is online

Users browsing this forum: No registered users and 0 guests

cron