-
-
Save cmtucker/8e5677bdd5c409d70738 to your computer and use it in GitHub Desktop.
The EEB & Flow's annual holiday caRd, 2015
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
##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) |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment