library(StatDA)
# need data:
data(chorizon)
data(ohorizon)
data(kola.background)
data(bordersKola)
X=ohorizon[,"XCOO"]
Y=ohorizon[,"YCOO"]
xwid=diff(range(X))/12e4
ywid=diff(range(Y))/12e4
pdf("fig-18-4.pdf",width=2*xwid,height=1*ywid)
par(mfrow=c(1,2),mar=c(1.5,1.5,1.5,1.5))
# Ca in O-Horizon
X=ohorizon[,"XCOO"]
Y=ohorizon[,"YCOO"]
# Fit variogram for data with non-standardized distance:
el="Ca"
###vario.b=varioCalc(X,Y,el=ohorizon[,el],max.dist=200000,title=paste(el,"in O-horizon"))
###res.eyefit=eyefit(vario.b)
###dump("res.eyefit","../Data/res.eyefit.Ca_O")
data(res.eyefit.Ca_O)
vario.b=varioCalc(X,Y,el=ohorizon[,el],max.dist=200000,title=paste(el,"in O-horizon"),km=F,
plot=FALSE)
v5=variofit(vario.b,res.eyefit.Ca_O,cov.model="matern",max.dist=200000,fix.nugget=T,
nugget=0.12,minimisation.function="nlm")
el=ohorizon[,"Ca"]
# Kriging
# generate plot with packground
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
KrigeLegend(X,Y,el,resol=50,vario=v5,type="percentile",whichcol="gray",
qutiles=c(0,0.05,0.25,0.50,0.75,0.90,0.98,1),borders="bordersKola",
leg.xpos.min=7.8e5,leg.xpos.max=8.0e5,leg.ypos.min=77.6e5,leg.ypos.max=78.7e5,
leg.title="mg/kg", leg.title.cex=0.7, leg.numb.cex=0.7, leg.round=2,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)
text(min(X)+diff(range(X))*4/7,max(Y),"Ca",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in O-horizon",cex=0.8)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)
# Ca in C-Horizon
X=chorizon[,"XCOO"]
Y=chorizon[,"YCOO"]
# Fit variogram for data with non-standardized distance:
el="Ca"
###vario.b=varioCalc(X,Y,el=chorizon[,el],max.dist=200000,title=paste(el,"in C-horizon"))
###res.eyefit=eyefit(vario.b)
###dump("res.eyefit","../Data/res.eyefit.Ca_C")
data(res.eyefit.Ca_C)
vario.b=varioCalc(X,Y,el=chorizon[,el],max.dist=200000,title=paste(el,"in C-horizon"),km=F,plot=FALSE)
v5=variofit(vario.b,res.eyefit.Ca_C,cov.model="spherical",max.dist=200000)
el=chorizon[,"Ca"]
# Kriging
# generate plot with packground
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
KrigeLegend(X,Y,el,resol=50,vario=v5,type="percentile",whichcol="gray",
qutiles=c(0,0.05,0.25,0.50,0.75,0.90,0.98,1),borders="bordersKola",
leg.xpos.min=7.8e5,leg.xpos.max=8.0e5,leg.ypos.min=77.6e5,leg.ypos.max=78.7e5,
leg.title="mg/kg", leg.title.cex=0.7, leg.numb.cex=0.7, leg.round=2,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)
text(min(X)+diff(range(X))*4/7,max(Y),"Ca",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/33,"in C-horizon",cex=0.8)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.5,shiftkm=37e3,sizetext=0.8)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7838750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)
dev.off()
library(StatDA)
data(chorizon)
data(ohorizon)
data(kola.background)
# Plot variogram
pdf("fig-18-5.pdf",width=9,height=4.5)
par(mfrow=c(1,2),mar=c(4,4,2,2))
# Ca in O-Horizon
X=ohorizon[,"XCOO"]
Y=ohorizon[,"YCOO"]
# Fit variogram for data with non-standardized distance:
el="Ca"
###vario.b=varioCalc(X,Y,el=ohorizon[,el],max.dist=200000,title=paste(el,"in O-horizon"))
###res.eyefit=eyefit(vario.b)
###dump("res.eyefit","../Data/res.eyefit.Ca_O")
data(res.eyefit.Ca_O)
vario.b=varioCalc(X,Y,el=ohorizon[,el],max.dist=200000,title=paste(el,"in O-horizon"),km=F,
plot=T)
v5=variofit(vario.b,res.eyefit.Ca_O,cov.model="matern",max.dist=200000,fix.nugget=T,
nugget=0.12,minimisation.function="nlm")
lines(v5,col=1,lwd=2)
# Ca in C-Horizon
X=chorizon[,"XCOO"]
Y=chorizon[,"YCOO"]
# Fit variogram for data with non-standardized distance:
el="Ca"
###vario.b=varioCalc(X,Y,el=chorizon[,el],max.dist=200000,title=paste(el,"in C-horizon"))
###res.eyefit=eyefit(vario.b)
###dump("res.eyefit","../Data/res.eyefit.Ca_C")
data(res.eyefit.Ca_C)
vario.b=varioCalc(X,Y,el=chorizon[,el],max.dist=200000,title=paste(el,"in C-horizon"),km=F,plot=T)
v5=variofit(vario.b,res.eyefit.Ca_C,cov.model="spherical",max.dist=200000)
lines(v5,col=1,lwd=2)
dev.off()