Creating a layered plotTree.wBars style fan-tree visualization for multiple quantitative traits

In response to a phytools user inquiry, I've been playing around with creating a multilayer plotTree.wBars style graph, showing the values of various continuous traits at the tips of a fan-style tree using bars.

One solution that I devised is to do this using thickened lines.

My solution is a bit wild & hacky, because to get the specific place & orientation of each line, I just continously transparently overplotted my tree, each time first lengthening all the tip edges, and then extracting the x, y coordinates of the new tip!

Here is the result, using some data for various morphological traits in phrynosomatid lizards from Harmon et al. (2010) that also feature in my upcoming book.

First load packages.

library(phytools)
library(geiger)

Then read tree and subsample data to match.

phryn.tree<-read.tree(file="http://www.phytools.org/Rbook/5/Phrynosomatidae.phy")
phryn.tree$edge.length<-phryn.tree$edge.length*100
iguania.data<-read.csv(file="http://www.phytools.org/Rbook/5/Iguania.csv",
    row.names=1)
chk<-name.check(phryn.tree,iguania.data)
summary(chk)

## 212 taxa are present in the data but not the tree:
##     amcris,
##     anacut,
##     anaen,
##     anahl,
##     analin,
##     analli,
##     ....
## 
## To see complete list of mis-matched taxa, print object.

phryn.data<-iguania.data[phryn.tree$tip.label,]

Now rescale data.

phryn.norm<-phryn.data/max(phryn.data,na.rm=TRUE)

Finally, create the plot.

h<-max(nodeHeights(phryn.tree))
m<-ncol(phryn.data)
cols<-RColorBrewer::brewer.pal(n=4,"Spectral")
xlim<-ylim<-c(-h,h)+c(-1,1)*0.1*m*h+0.02*c(-h,h)
plotTree(phryn.tree,type="fan",xlim=xlim,ylim=ylim,
    lwd=1,ftype="off")
for(i in 1:m){
    tt<-phryn.tree
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
        tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+
        0.02*h+(i-1)*0.1*h
    plotTree(tt,color="transparent",type="fan",xlim=xlim,
        ylim=ylim,lwd=1,ftype="off",add=TRUE)
    pp1<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]<-
        tt$edge.length[which(tt$edge[,2]<=Ntip(tt))]+
        0.09*h
    plotrix::draw.circle(0,0,radius=h+0.02*h+(i-1)*0.1*h,
        border="grey")
    plotTree(tt,color="transparent",type="fan",xlim=xlim,
        ylim=ylim,lwd=1,ftype="off",add=TRUE)
    pp2<-get("last_plot.phylo",envir=.PlotPhyloEnv)
    par(lend=1)
    for(j in 1:Ntip(phryn.tree)){
        ii<-which(rownames(phryn.norm)==phryn.tree$tip.label[j])
        dx<-(pp2$xx[j]-pp1$xx[j])*phryn.norm[ii,i]
        dy<-(pp2$yy[j]-pp1$yy[j])*phryn.norm[ii,i]
        #lines(pp1$xx[j]+c(0,1.05*dx),pp1$yy[j]+c(0,1.05*dy),lwd=10,
        #   col=par()$fg)
        lines(pp1$xx[j]+c(0,dx),pp1$yy[j]+c(0,dy),lwd=10,
            col=cols[i])
    }
}
xx<-rep(0.85*par()$usr[4],m)
yy<-0.95*par()$usr[4]-1.5*0:(m-1)*strheight("W")
text(xx,yy,colnames(phryn.norm),pos=4,cex=0.8)
scale.bar<-apply(phryn.norm,2,max,na.rm=TRUE)*0.09*h
xx2<-xx-scale.bar
segments(x0=xx,y0=yy,x1=xx2,y1=yy,lwd=10,col=cols)
text(rep(min(xx2),m),yy,paste(formatC(apply(phryn.data,2,max,
    na.rm=TRUE),digits=1,format="f"),"cm",sep=" "),cex=0.8,
    pos=2)