# function depends on phytools (& dependencies) and maps (& dependencies) # written by Liam J. Revell 2013 phylo.to.map<-function(tree,coords,rotate=TRUE){ # open & size a new plot par(mai=c(0.1,0.1,0.1,0.1)) # map<-getMap(resolution="low") map("world",ylim=c(-90,180)) # 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(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))*(90-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)*360-180 # 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<-180-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 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) coords<-coords[tree$tip.label,2:1] points(coords,pch=16,cex=1,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") } # 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) }