library(StatDA)
data(ohorizon)
data(chorizon)
data(kola.background)
# only samples from same locations in both layers:
c.ind=NULL
o.ind=NULL
c.id=chorizon[,1]
o.id=ohorizon[,1]
for (i in 1:1000){
if (sum(c.id==i)+sum(o.id==i)==2){
c.ind=c(c.ind,which(c.id==i))
o.ind=c(o.ind,which(o.id==i))
}
}
X=ohorizon[o.ind,"XCOO"]
Y=ohorizon[o.ind,"YCOO"]
# true representation of x and y axis of map for plot
xwid=diff(range(X))/12e4
ywid=diff(range(Y))/12e4
pdf("fig-7-7.pdf",width=2*xwid,height=3*ywid)
par(mfrow=c(3,2),mar=c(1.5,1.5,1.5,1.5))
### NEW MAP:
elO=(ohorizon[o.ind,"Cu"])
elC=(chorizon[c.ind,"Cu"])
el=elO/elC
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Ratio Cu_O/Cu_C",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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)
### NEW MAP:
elO=(ohorizon[o.ind,"Pb"])
elC=(chorizon[c.ind,"Pb"])
elO=elO[!is.na(elC)]
elC=elC[!is.na(elC)]
el=elO/elC
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Ratio Pb_O/Pb_C",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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)
###############################################################################
### NEW MAP:
elO=(ohorizon[o.ind,"Cu"])
elC=(chorizon[c.ind,"Cu"])
el=elO
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Cu_O [mg/kg]",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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)
###############################################################################
### NEW MAP:
elO=(ohorizon[o.ind,"Pb"])
elC=na.omit(chorizon[c.ind,"Pb"])
elO=elO[!is.na(elC)]
elC=elC[!is.na(elC)]
el=elO
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Pb_O [mg/kg]",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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)
###############################################################################
### NEW MAP:
elO=(ohorizon[o.ind,"Cu"])
elC=(chorizon[c.ind,"Cu"])
el=elC
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Cu_C [mg/kg]",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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)
###############################################################################
### NEW MAP:
elO=(ohorizon[o.ind,"Pb"])
elC=na.omit(chorizon[c.ind,"Pb"])
elO=elO[!is.na(elC)]
elC=elC[!is.na(elC)]
el=elC
# generate plot with background
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")
plotbg(map.col=c("gray","gray","gray","gray"),add.plot=T)
SymbLegend(X,Y,el,type="percentile",qutiles<-c(0,0.05,0.25,0.75,0.95,1),symbtype="EDA",symbmagn=0.8,
leg.position="topright",leg.title="Pb_C [mg/kg]",leg.title.cex=0.95,leg.round=2,leg.wid=4,leg.just="right")
# Text
text(min(X)+diff(range(X))*4/7,max(Y),paste(qutiles*100,collapse=","),cex=1.0)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/25,"Percentiles",cex=1.0)
# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.2,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()