# function depends on phytools (& dependencies) and maps (& dependencies) # written by Liam J. Revell 2013 phylo.to.map<-function(tree,coords,rotate=TRUE,...){ # get optional arguments if(hasArg(xlim)) xlim<-list(...)$xlim else xlim<-c(-180,180) if(hasArg(ylim)) ylim<-list(...)$ylim else ylim<-c(-90,90) if(hasArg(fsize)) fsize<-list(...)$fsize else fsize<-1.0 if(hasArg(split)) split<-list(...)$split else split<-c(0.42,0.58) if(hasArg(psize)) psize<-list(...)$psize else psize<-1.0 if(hasArg(plot)) plot<-list(...)$plot else plot<-TRUE # create a map map<-map("world",xlim=xlim,ylim=ylim,plot=FALSE,fill=TRUE,resolution=0) # recompute ylim to leave space for the tree ylim<-c(ylim[1],ylim[2]+split[1]/split[2]*(ylim[2]-ylim[1])) # open & size a new plot par(mar=rep(0,4)) plot.new() plot.window(xlim=xlim,ylim=ylim,asp=1) map(map,add=TRUE,fill=TRUE,col="gray95",mar=rep(0,4)) # plot a white rectangle dx<-abs(diff(xlim)) rect(xlim[1]-1.04*dx,ylim[2]-split[1]*(ylim[2]-ylim[1]),xlim[2]+1.04*dx,ylim[2],col="white",border="white") # if rotate if(rotate) tree<-minRotate(tree,coords[,2]) # rescale tree so it fits in the upper half of the plot # with enough space for labels sh<-max(fsize*strwidth(tree$tip.label))/(par()$usr[2]-par()$usr[1])*(par()$usr[4]-par()$usr[3]) tree$edge.length<-tree$edge.length/max(nodeHeights(tree))*(split[1]*(ylim[2]-ylim[1])-sh) n<-length(tree$tip.label) # reorder cladewise to assign tip positions cw<-reorder(tree,"cladewise") x<-vector(length=n+cw$Nnode) x[cw$edge[cw$edge[,2]<=n,2]]<-0:(n-1)/(n-1)*(xlim[2]-xlim[1])+xlim[1] # reorder pruningwise for post-order traversal pw<-reorder(tree,"pruningwise") nn<-unique(pw$edge[,1]) # compute horizontal position of each edge for(i in 1:length(nn)){ xx<-x[pw$edge[which(pw$edge[,1]==nn[i]),2]] x[nn[i]]<-mean(range(xx)) } # compute start & end points of each edge Y<-ylim[2]-nodeHeights(cw) # plot vertical edges for(i in 1:nrow(Y)) lines(rep(x[cw$edge[i,2]],2),Y[i,],lwd=2,lend=2) # plot horizontal relationships for(i in 1:tree$Nnode+n) lines(range(x[cw$edge[which(cw$edge[,1]==i),2]]),Y[which(cw$edge[,1]==i),1],lwd=2,lend=2) # plot coordinates & linking lines coords<-coords[tree$tip.label,2:1] points(coords,pch=16,cex=psize,col="red") for(i in 1:n) lines(c(x[i],coords[i,1]),c(Y[which(cw$edge[,2]==i),2],coords[i,2]),col="red",lty="dashed") # plot tip labels for(i in 1:n) text(x[i],Y[which(cw$edge[,2]==i),2],tree$tip.label[i],pos=4,offset=0.1,srt=-90,cex=fsize) } # rotates all nodes to try and match tip ordering to longitude # written by Liam J. Revell 2013 minRotate<-function(tree,x){ tree<-reorder(tree) nn<-1:tree$Nnode+length(tree$tip.label) x<-x[tree$tip.label] for(i in 1:tree$Nnode){ tt<-read.tree(text=write.tree(rotate(tree,nn[i]))) oo<-sum(abs(order(x[tree$tip.label])-1:length(tree$tip.label))) pp<-sum(abs(order(x[tt$tip.label])-1:length(tt$tip.label))) if(oo>pp) tree<-tt message(paste("objective:",min(oo,pp))) } return(tree) }