On 01/30/2012 11:09 PM, Roopashree Shrivastava wrote:
Dear all,
I am new to plotting Taylor Diagram using plotrix package within R, hence
this post. I have written a script which plots Taylor Diagram with one
reference and 7 model values. However the font size, line width and line
type are not clear when saving the diagram as a jpeg file. I tried the
functions lty, lwd and font but no apparent change. I am attaching the
script here. Any help would be greatly appreciated. The script is...
Hi Roopa,
If by "not clear" you mean that the figure (especially the text) is
generally a bit fuzzy, you probably want to use a device like PDF or
Postscript.
If you mean that you want larger text on the axes, I have added an
argument "cex.axis" to the function below that allows the user to get a
bit of text expansion. This will appear in the next version of plotrix.
About line width and type, there seems to be a convention on the line
types used in Taylor diagrams, but additional arguments could be added
to modify these if you really need them. If you want to change the font,
do it as in the example below.
I found that using only the positive correlations produced a much
clearer plot, so you may want to try this somewhat altered example:
ref<-c(0.00640091,0.00533091,0.00381636,0.00275519,0.00277649,
0.00280806,0.00267945,0.00237123,0.000970663,0.000986191,
0.00100226,0.00086391,0.000622819,0.000485319,
0.000362976,0.000246112,0.000165615,0.00008184,0.00004)
m1<-c(0.0124827,0.011662,0.0102956,0.0091183,0.00813907,
0.007192,0.00662517,0.00433745,0.00184044,0.000649477,
0.00024642,5.43E-05,0.000097696,0.000194817,
0.000182709,0.000134398,0.000106024,8.92E-05,6.28E-05)
m2<-c(0.0101348,0.00920886,0.0086196,0.00785134,0.00723838,
0.00675833,0.00579093,0.00540478,0.00226489,0.000809049,
0.00019625,3.95E-05,8.89E-05,0.000195028,
0.000185004,0.000131202,0.000109852,9.98E-05,6.80E-05)
m3<-c(0.0123251,0.0120384,0.00871793,0.00678519,0.00628331,
0.00532673,0.00486861,0.0048328,0.0038655,0.00143683,
0.00022057,8.61E-06,7.79E-05,0.000184976,
0.000185927,0.000133771,0.000104613,9.26E-05,6.38E-05)
m4<-c(0.0134251,0.0126776,0.012559,0.0121933,0.0099911,
0.00727952,0.00475407,0.00227909,0.00130748,0.000705607,
0.000304828,5.70E-05,0.000109972,0.000187504,
0.0002016,0.000133706,0.000109697,9.54E-05,6.35E-05)
m5<-c(0.0124275,0.0112242,0.00886243,0.00793019,0.0067846,
0.00603205,0.00566561,0.00530552,0.00318331,0.000961854,
0.000218234,3.66E-05,7.99E-05,0.000182724,
0.000196627,0.000136862,0.000104907,0.000094622,6.20E-05)
m6<-c(0.0142817,0.0134474,0.0129694,0.0113914,0.0102208,
0.00920309,0.00555206,0.00289796,0.00143831,0.000706277,
0.000277201,5.60E-05,0.000114714,0.000186412,
0.000198743,0.000134991,0.000108689,9.43E-05,6.16E-05)
m7<-c(0.0120621,0.0117936,0.00854782,0.00734006,0.00669576,
0.00629334,0.00595018,0.00564455,0.00396859,0.000991006,
0.000171742,9.68E-06,8.10E-05,0.000186982,
0.00018854,0.000136548,0.000104581,9.18E-05,6.20E-05)
par(family="serif")
taylor.diagram(ref,m1,ngamma=3,pch=19,pcex=2,cex.axis=1.2,
lty=1,lwd=10,font=5)
taylor.diagram(ref,m2,add=TRUE,pch=19,pcex=2,col="blue",
lty="solid",lwd=3)
taylor.diagram(ref,m3,add=TRUE,pch=19,pcex=2,col="orange",
lty="solid",lwd=3)
taylor.diagram(ref,m4,add=TRUE,pch=19,pcex=2,col="pink",
lty="solid",lwd=3)
taylor.diagram(ref,m5,add=TRUE,pch=19,pcex=2,col="purple",
lty="solid",lwd=3)
taylor.diagram(ref,m6,add=TRUE,pch=19,pcex=2,col="brown",
lty="solid",lwd=3)
taylor.diagram(ref,m7,add=TRUE,pch=19,pcex=2,col="cyan",
lty="solid",lwd=3)
lpos<-3.5*sd(ref)
legend(lpos,1.1*lpos,
legend=c("1111","1171","1211","2121","2221","4141","5251"),
pch=19,
col=c("red","blue","orange","pink","purple","brown","cyan"))
Jim
# Display a Taylor diagram
# Taylor K.E. (2001)
# Summarizing multiple aspects of model performance in a single diagram
# Journal of Geophysical Research, 106: 7183-7192.
# version 1.0
# progr. Olivier.Eterradossi, 12/2007
# 2007-01-12 - modifications and Anglicizing - Jim Lemon
# version 2.0
# progr. initiale OLE, 8/01/2007
# rev. OLE 3/09/2008 : remove samples with NA's from mean, sd and cor
calculations
# 2008-09-04 - integration and more anglicizing - Jim Lemon
# 2008-12-09 - added correlation radii, sd arcs to the pos.cor=FALSE routine
# and stopped the pos.cor=FALSE routine from calculating arcs for zero
radius
# Jim Lemon
# 2010-4-30 - added the gamma.col argument for pos.cor=TRUE plots - Jim
Lemon
# 2010-6-24 - added mar argument to pos.cor=TRUE plots - Jim Lemon
# 2012-1-31 - added the cex.axis argument - Jim Lemon
taylor.diagram<-function(ref,model,add=FALSE,col="red",pch=19,pos.cor=TRUE,
xlab="",ylab="",main="Taylor
Diagram",show.gamma=TRUE,ngamma=3,gamma.col=8,
sd.arcs=0,ref.sd=FALSE,grad.corr.lines=c(0.2,0.4,0.6,0.8,0.9),pcex=1,
cex.axis=1,normalize=FALSE,mar=c(5,4,6,6),...) {
grad.corr.full<-c(0,0.2,0.4,0.6,0.8,0.9,0.95,0.99,1)
R<-cor(ref,model,use="pairwise")
sd.r<-sd(ref)
sd.f<-sd(model)
if(normalize) {
sd.f<-sd.f/sd.r
sd.r<-1
}
maxsd<-1.5*max(sd.f,sd.r)
oldpar<-par("mar","xpd","xaxs","yaxs")
if(!add) {
# display the diagram
if(pos.cor) {
if(nchar(ylab) == 0) ylab="Standard deviation"
par(mar=mar)
plot(0,xlim=c(0,maxsd),ylim=c(0,maxsd),xaxs="i",yaxs="i",axes=FALSE,
main=main,xlab=xlab,ylab=ylab,type="n",cex=cex.axis,...)
if(grad.corr.lines[1]) {
for(gcl in grad.corr.lines)
lines(c(0,maxsd*gcl),c(0,maxsd*sqrt(1-gcl^2)),lty=3)
}
# add the axes
segments(c(0,0),c(0,0),c(0,maxsd),c(maxsd,0))
axis.ticks<-pretty(c(0,maxsd))
axis.ticks<-axis.ticks[axis.ticks<=maxsd]
axis(1,at=axis.ticks,cex.axis=cex.axis)
axis(2,at=axis.ticks,cex.axis=cex.axis)
if(sd.arcs[1]) {
if(length(sd.arcs) == 1) sd.arcs<-axis.ticks
for(sdarc in sd.arcs) {
xcurve<-cos(seq(0,pi/2,by=0.03))*sdarc
ycurve<-sin(seq(0,pi/2,by=0.03))*sdarc
lines(xcurve,ycurve,col="blue",lty=3)
}
}
if(show.gamma[1]) {
# if the user has passed a set of gamma values, use that
if(length(show.gamma) > 1) gamma<-show.gamma
# otherwise make up a set
else gamma<-pretty(c(0,maxsd),n=ngamma)[-1]
if(gamma[length(gamma)] > maxsd) gamma<-gamma[-length(gamma)]
labelpos<-seq(45,70,length.out=length(gamma))
# do the gamma curves
for(gindex in 1:length(gamma)) {
xcurve<-cos(seq(0,pi,by=0.03))*gamma[gindex]+sd.r
# find where to clip the curves
endcurve<-which(xcurve<0)
endcurve<-ifelse(length(endcurve),min(endcurve)-1,105)
ycurve<-sin(seq(0,pi,by=0.03))*gamma[gindex]
maxcurve<-xcurve*xcurve+ycurve*ycurve
startcurve<-which(maxcurve>maxsd*maxsd)
startcurve<-ifelse(length(startcurve),max(startcurve)+1,0)
lines(xcurve[startcurve:endcurve],ycurve[startcurve:endcurve],
col=gamma.col)
if(xcurve[labelpos[gindex]] > 0)
boxed.labels(xcurve[labelpos[gindex]],ycurve[labelpos[gindex]],
gamma[gindex],border=FALSE)
}
}
# the outer curve for correlation
xcurve<-cos(seq(0,pi/2,by=0.01))*maxsd
ycurve<-sin(seq(0,pi/2,by=0.01))*maxsd
lines(xcurve,ycurve)
bigtickangles<-acos(seq(0.1,0.9,by=0.1))
medtickangles<-acos(seq(0.05,0.95,by=0.1))
smltickangles<-acos(seq(0.91,0.99,by=0.01))
segments(cos(bigtickangles)*maxsd,sin(bigtickangles)*maxsd,
cos(bigtickangles)*0.97*maxsd,sin(bigtickangles)*0.97*maxsd)
par(xpd=TRUE)
if(ref.sd) {
# the inner curve for reference SD
xcurve<-cos(seq(0,pi/2,by=0.01))*sd.r
ycurve<-sin(seq(0,pi/2,by=0.01))*sd.r
lines(xcurve,ycurve)
}
points(sd.r,0,cex=pcex)
text(cos(c(bigtickangles,acos(c(0.95,0.99))))*1.05*maxsd,
sin(c(bigtickangles,acos(c(0.95,0.99))))*1.05*maxsd,
c(seq(0.1,0.9,by=0.1),0.95,0.99))
text(maxsd*0.8,maxsd*0.8,"Correlation",srt=315)
segments(cos(medtickangles)*maxsd,sin(medtickangles)*maxsd,
cos(medtickangles)*0.98*maxsd,sin(medtickangles)*0.98*maxsd)
segments(cos(smltickangles)*maxsd,sin(smltickangles)*maxsd,
cos(smltickangles)*0.99*maxsd,sin(smltickangles)*0.99*maxsd)
}
else {
x<- ref
y<- model
R<-cor(x,y,use="pairwise.complete.obs")
E<-mean(x,na.rm=TRUE)-mean(y,na.rm=TRUE) # overall bias
xprime<-x-mean(x,na.rm=TRUE)
yprime<-y-mean(y,na.rm=TRUE)
sumofsquares<-(xprime-yprime)^2
Eprime<-sqrt(sum(sumofsquares)/length(complete.cases(x))) # centered
pattern RMS
E2<-E^2+Eprime^2
if (add==FALSE) {
# pourtour du diagramme (display the diagram)
maxray<-1.5*max(sd.f,sd.r)
plot(c(-maxray,maxray),c(0,maxray),type="n",asp=1,bty="n",xaxt="n",yaxt="n",
xlab=xlab,ylab=ylab,main=main,cex=cex.axis)
discrete<-seq(180,0,by=-1)
listepoints<-NULL
for (i in discrete){
listepoints<-cbind(listepoints,maxray*cos(i*pi/180),maxray*sin(i*pi/180))
}
listepoints<-matrix(listepoints,2,length(listepoints)/2)
listepoints<-t(listepoints)
lines(listepoints[,1],listepoints[,2])
# axes x,y
lines(c(-maxray,maxray),c(0,0))
lines(c(0,0),c(0,maxray))
# lignes radiales jusqu'� R = +/- 0.8
for (i in grad.corr.lines){
lines(c(0,maxray*i),c(0,maxray*sqrt(1-i^2)),lty=3)
lines(c(0,-maxray*i),c(0,maxray*sqrt(1-i^2)),lty=3)
}
# texte radial
for (i in grad.corr.full){
text(1.05*maxray*i,1.05*maxray*sqrt(1-i^2),i,cex=0.6)
text(-1.05*maxray*i,1.05*maxray*sqrt(1-i^2),-i,cex=0.6)
}
# sd concentriques autour de la reference
seq.sd<-seq.int(0,2*maxray,by=(maxray/10))[-1]
for (i in seq.sd){
xcircle<-sd.r+(cos(discrete*pi/180)*i)
ycircle<-sin(discrete*pi/180)*i
for (j in 1:length(xcircle)){
if ((xcircle[j]^2+ycircle[j]^2)<(maxray^2)){
points(xcircle[j],ycircle[j], col="darkgreen",pch=".")
if(j==10)
text(xcircle[j],ycircle[j],signif(i,2),cex=0.5,col="darkgreen")
}
}
}
# sd concentriques autour de l'origine
seq.sd<-seq.int(0,maxray,length.out=5)
for (i in seq.sd){
xcircle<-(cos(discrete*pi/180)*i)
ycircle<-sin(discrete*pi/180)*i
if(i) lines(xcircle,ycircle,lty=3,col="blue")
text(min(xcircle),-0.03*maxray,signif(i,2),cex=0.5,col="blue")
text(max(xcircle),-0.03*maxray,signif(i,2),cex=0.5,col="blue")
}
text(0,-0.08*maxray,"Standard Deviation",cex=0.7,col="blue")
text(0,-0.12*maxray,"Centered RMS Difference",cex=0.7,col="darkgreen")
points(sd.r,0,pch=22,bg="darkgreen",cex=1.1)
text(0,1.1*maxray,"Correlation Coefficient",cex=0.7)
}
S<-(2*(1+R))/(sd.f+(1/sd.f))^2
# Taylor<-S
}
}
# display the points
points(sd.f*R,sd.f*sin(acos(R)),pch=pch,col=col,cex=pcex)
invisible(oldpar)
}
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
