5.8 Colour Maps

# Fig. 5.12.: Color smoothed surface maps
library(StatDA)
data(chorizon)
X=chorizon[,"XCOO"]
Y=chorizon[,"YCOO"]
el=log10(chorizon[,"As"])
data(kola.background)
data(bordersKola)

# true representation of x and y axis of map for plot
xwid=diff(range(X))/12e4
ywid=diff(range(Y))/12e4


pdf("fig-5-12.pdf",width=2*xwid*0.8,height=3*ywid*0.8)
par(mfrow=c(3,2),mar=c(1.5,1.5,1.5,1.5))

######################################################################################
# Rainbow map based on Percentiles:

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="percentile",whichcol="rainbow",
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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################

# Inverted rainbow map based on Percentiles:

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="percentile",whichcol="rainbow.inv",
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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################

# Terrain map based on Percentiles:

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="percentile",whichcol="terrain",
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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################

# Topo map based on Percentiles:

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="percentile",whichcol="topo",
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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################

# Rainbow map based on continuous scale (Finnish):

## Finnish color scale:
#nc=98 # number of classes
#Q10=quantile(whichdraw[!is.na(whichdraw)],0.10)
#const=floor(Q10)
#Q99=quantile(whichdraw[!is.na(whichdraw)],0.99)
#im.br=10^( log10(Q99-const)*(1:nc-1)/(nc-2) + log10(Q10-const)*(nc-1-1:nc)/(nc-2) ) +const
#im.br=c(min(whichdraw[!is.na(whichdraw)]),im.br,max(whichdraw[!is.na(whichdraw)]))
#im.col=rev(rainbow(length(im.br)-1,start=0,end=0.7))

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="contin",whichcol="rainbow",
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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################

# Rainbow map based on 75% green

# generate plot
plot(X,Y,frame.plot=FALSE,xaxt="n",yaxt="n",xlab="",ylab="",type="n")

SmoothLegend(X,Y,el,resol=200,type="percentile",whichcol="rainbow.trunc",
qutiles=c(0,0.75,0.80,0.85,0.90,0.95,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.90, leg.numb.cex=0.90, leg.round=2,leg.wid=4,
leg.numb.xshift=0.7e5,leg.perc.xshift=0.4e5,leg.perc.yshift=0.2e5,tit.xshift=0.35e5)

# plot background
plotbg(map.col=c("gray","gray","gray","gray"),map.lwd=c(1,1,1,1),add.plot=T)

# Text
text(min(X)+diff(range(X))*4/7,max(Y),"As",cex=1)
text(min(X)+diff(range(X))*4/7,max(Y)-diff(range(Y))/28,"in C-horizon",cex=0.9)

# scalebar
scalebar(761309,7373050,861309,7363050,shifttext=-0.1,shiftkm=37e3,sizetext=0.75)
# North arrow
Northarrow(362602,7818750,362602,7878750,362602,7831750,Alength=0.15,Aangle=15,Alwd=1.3,Tcex=1.6)

######################################################################################
dev.off()