5.7.1 Construction of the (Semi)Variogram

# Fig. 5.7.: Basic information for variogram computation library(StatDA) # variogram example data(chorizon) data(kola.background) # Krige results: pdf("fig-5-7.pdf",width=8, height=8) par(mfrow=c(2,2),mar=c(4,4,2,2)) # As in C-Horizon X=chorizon[,"XCOO"]/1000 Y=chorizon[,"YCOO"]/1000 el=chorizon[,"As"] # binned variogram vario.b <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300) # variogram cloud vario.c <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, op="cloud") #binned variogram and stores the cloud vario.bc <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, bin.cloud=TRUE) # smoothed variogram vario.s <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, op="sm", band=10) # plotting the variograms: plotvario(vario.c,xlab="Distance [km]", ylab="Semivariogram",cex.lab=1.2,max.dist=300,pch=3,cex=0.5) plotvario(vario.bc, bin.cloud=TRUE, xlab="Distance classes",ylab="Semivariogram",cex.lab=1.2) plotvario(vario.s, type="l", xlab="Distance [km]",ylab="Semivariogram",cex.lab=1.2,max.dist=300) plotvario(vario.b,xlab="Distance [km]",ylab="Semivariogram",cex.lab=1.2,max.dist=300,pch=1,cex=1) dev.off()

# Fig. 5.8.: Variogram in different directions and model library(StatDA) data(chorizon) data(kola.background) # Krige results: pdf("fig-5-8.pdf",width=9, height=4.5) par(mfrow=c(1,2),mar=c(4,4,2,2)) # As in C-horizon X=chorizon[,"XCOO"]/1000 Y=chorizon[,"YCOO"]/1000 el=chorizon[,"As"] # binned variogram vario.b <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300) # variogram cloud vario.c <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, op="cloud") #binned variogram and stores the cloud vario.bc <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, bin.cloud=TRUE) # smoothed variogram vario.s <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, op="sm", band=10) # computing a directional variogram plot(0,0,xlab="Distance [km]",ylab="Semivariogram",xlim=c(0,300),ylim=c(0,1.6), cex.lab=1.2,type="n") title("As in C-horizon",cex.main=1.2) lines(vario.b,pch=1,type="p") vario.0 <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, dir=0, tol=pi/8) vario.90 <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, dir=pi/4, tol=pi/8) vario.45 <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, dir=pi/8, tol=pi/8) vario.120 <- variog(coords=cbind(X,Y), data=el, lambda=0, max.dist=300, dir=3*pi/8, tol=pi/8) lines(vario.0,lty=1,pch=2) lines(vario.90,lty=2,pch=3) lines(vario.45,lty=3,pch=4) lines(vario.120,lty=4,pch=5) legend("bottomright", legend=c("omnidirectional","N-S","E-W","NW-SE","NE-SW"), pch=c(1,2,3,4,5)) # Fit variogram for data with non-standardized distance: #res.eyefit=eyefit(vario.b) #dump("res.eyefit","../Data/res.eyefit.As_C") data(res.eyefit.As_C) v5=variofit(vario.b,res.eyefit.As_C,cov.model="spherical",max.dist=300) plot(0,0,xlab="Distance [km]",ylab="Semivariogram",xlim=c(0,300),ylim=c(0,1.6), cex.lab=1.2,type="n") title("As in C-horizon",cex.main=1.2) lines(vario.b,pch=1,type="p") lines(v5,col=1,lwd=2) #Range r=v5$cov.pars[2] # Nugget n=v5$nugget # Sill s=v5$cov.pars[1] arrows(0,0,0,n,length=0.08,code=3,col=gray(0.6)) text(2,n/2,paste("Nugget variance =",round(n,2)),cex=0.9,pos=4) abline(h=n,col=gray(0.6),lty=2) arrows(300,n,300,n+s,length=0.08,code=3,col=gray(0.6)) text(298,n+s/2,paste("Sill =",round(s,2)),cex=0.9,pos=2) arrows(0,1.3,r,1.3,length=0.08,code=3,col=gray(0.6)) text(r/2,1.3,paste("Range =",round(r,0)),cex=0.9,pos=3) dev.off()