The EEB & Flow's annual holiday caRd, 2015

caRd2015.R

##2015 Holiday caRd from the EEB & Flow, by C. Tucker

##If you don't have these libraries, please install before running...
#I didn't test this on different OS/versions of R this year, so if you run into problems refer to the blog post instead :-)

##Req'd libraries
require(ape)
require(rworldmap)
require(FD)


# function modified from on phytools (& dependencies) and rworldmap (& dependencies) (written by Liam J. Revell 2013)
phylo.to.map2<-function(tree,coords){
	# open & size a new plot
	par(mai=c(0.1,0.1,0.1,0.1))
	map<-getMap(resolution="low")
	plot(map,ylim=c(-30,100))
	# 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)+10),(Y[i,]+10),lwd=2,lend=2, col="forestgreen")
	# plot horizontal relationships
	for(i in 1:tree$Nnode+n)
		lines((range(x[cw$edge[which(cw$edge[,1]==i),2]])+10),(Y[which(cw$edge[,1]==i),1]+10),lwd=2,lend=2, col="forestgreen")
	# plot tip labels
	lab_short <- c("Coca \ncola Santa","Mall \nSanta","Salvation \nArmy.Santa","Kris.Kringle", "Santa \nClaus","Pere.Noel","Sinterklaas","Saint \nNicholas", "Father \nChristmas", "Ded \nMoroz") 
	for(i in 1:n)
		text(x[i],Y[which(cw$edge[,2]==i),2],lab_short[i],srt=0, col="red", cex=0.5)
	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")
}

#also from phytools
nodeHeights <- function (tree) 
{
    if (attr(tree, "order") != "cladewise" || is.null(attr(tree, 
        "order"))) 
        t <- reorder(tree)
    else t <- tree
    root <- length(t$tip.label) + 1
    X <- matrix(NA, nrow(t$edge), 2)
    for (i in 1:nrow(t$edge)) {
        if (t$edge[i, 1] == root) {
            X[i, 1] <- 0
            X[i, 2] <- t$edge.length[i]
        }
        else {
            X[i, 1] <- X[match(t$edge[i, 1], t$edge[, 2]), 2]
            X[i, 2] <- X[i, 1] + t$edge.length[i]
        }
    }
    if (attr(tree, "order") != "cladewise" || is.null(attr(tree, 
        "order"))) 
        o <- apply(matrix(tree$edge[, 2]), 1, function(x, y) which(x == 
            y), y = t$edge[, 2])
    else o <- 1:nrow(t$edge)
    return(X[o, ])
}
########################

##Make tree
edge <- matrix(c(11,12,12,13,13,14,14,9,14,15,15,6,15,7,13,8, 12, 10,11,16,16,17,17,18,18,1,18,2,17,3,16,19,19,4,19,5), nrow=18, ncol=2, byrow=TRUE)
edge.length <- c(0.1, 0.4, 0.1, 0.5, 0.3, 0.2, 0.2, 0.6, 1.0, 0.2, 0.3, 0.3, 0.3, 0.3, 0.6, 0.1, 0.8, 0.8)
tip.label <- c("Coca.cola.Santa", "Mall.Santa", "Salvation.Army.Santa", "Kriss.Kringle", "Santa.Claus",    "Pere.Noel", "Father.Christmas", "Sinterklaas", "Saint.Nicholas", "Ded.Moroz")
Nnode <- 9
xmas.tree <- list(edge, edge.length, Nnode, tip.label)
class(xmas.tree) <- "phylo"
names(xmas.tree) <- c("edge", "edge.length", "Nnode", "tip.label")

####Start

#taxa exploration
#functional ecology
heft <- c("fat", "fat", "fat", "fat", "fat", "thin", "fat","thin", "thin", "thin")

transport <- c("reindeer", "reindeer", "reindeer", "reindeer", "reindeer",  "donkey", "foot", "horse", "foot", "foot")
first.appearance <- c(1900, 1900, 1890, 1800, 1700, 1400, 1400, 400, 300, 1937)

santatraits <- cbind.data.frame(heft, transport, first.appearance)
rownames(santatraits) <- xmas.tree$tip.label

###Card plotting
#Make device with set size, works on multiple OS
dev.new <- function(width = 5, height = 5){
	platform <- sessionInfo()$platform 
if (grepl("linux",platform)) { 
	x11(width=width, height=height) 
	}else{ 
if (grepl("pc",platform)) { 
	windows(width=width, height=height) 
	}else{
if (grepl("w32",platform)) { 
	windows(width=width, height=height) 
	}else{		
if (grepl("apple", platform)){ 
	quartz(width=width, height=height)} }}}}

dev.new(10, 6)
par(mfrow=c(1,2), bg="darkseagreen1", mai=c(0.1,0.1,0.1,0.1), oma=c(0.1,0.1,0.1,0.1))

plot(xmas.tree, type = "c", FALSE, edge.color="darkgreen",  edge.lty=1, edge.width=18, label.offset = 1, direction="downward", font=3, tip.color="darkred")

#reconstruct ancestral state
cc <- ace(transport, xmas.tree, type="discrete")
co2 <- c("yellow", "gold", "darkorange", "red")
nodelabels(pie = cc$lik.anc, piecol = co2, cex = c(1.5, rep(1, 8)))

#Plot traits (fatness) against Santa
co1 <- c("blue", "purple")
tiplabels(pch = 12, col = co1[as.factor(heft)], cex = 3.5, adj=c(0.5, 0), lwd=2)

#Let's see the transportation mode trait too:
tiplabels(pch = 8, col = co2[as.factor(transport)], cex = 2, adj=c(0.5, 0), lwd=2)
#where do santa taxa live?

locales <- cbind(c(40.7, 40.7, 40.7, 70.93, 70.93, 47, 51.4, 52.2, 39.8, 55.7), c( -73.93, -73.93, -73.93, -95.4, -95.4,  2.7, -0.85, 5.7, 33.3, 37.6))
rownames(locales) <- xmas.tree$tip.label

phylo.to.map2(xmas.tree, locales)
text(x=-20, y=-120, "Happy Holidays from   \nthe EEB & Flow", col="darkred", cex=2, font=4)