Saturday, June 4, 2016

Some updates to phylo.to.map plotting method

I just finished adding some features to the phytools function phylo.to.map for projecting a phylogeny onto a geographic map, and, in particular, to the S3 plotting method.

Mostly, I have added features that I think better automate the spacing of tree, tree labels, links, and the map; but I also contributed some attributes that are aesthetic decisions regarding how the tree and plotted maps look. For instance, I changed the default point type from pch=19 to pch=21, and, by default, I added points to the tips of the tree - although these can be removed using the argument pts.

Let me know if you encounter any new bugs that might be due to this update!

Here is a demo using totally simulated, but somewhat realistic looking data:

library(phytools)
library(mapdata)

## load source code of function from GitHub
## also can be obtained by updating phytools from GitHub using devtools
source("https://raw.githubusercontent.com/liamrevell/phytools/master/R/phylo.to.map.R")

## first a world tree & dataset using default settings:
obj<-phylo.to.map(World.tree,World,plot=FALSE)
## objective: 130
## objective: 122
## objective: 122
## objective: 122
## objective: 118
## objective: 118
## objective: 118
## objective: 116
## objective: 116
## objective: 116
## objective: 104
## objective: 104
## objective: 100
## objective: 100
## objective: 100
## objective: 96
## objective: 96
## objective: 96
## objective: 96
## objective: 92
## objective: 88
## objective: 88
## objective: 88
## objective: 88
## objective: 88
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 26 tips and 25 internal nodes.
## 
## (2) A geographic map with range:
##      -85.19N, 83.6N
##      -180W, 180W.
## 
## (3) A table containing 26 geographic coordinates.
plot(obj,ftype="i",psize=1.5)

plot of chunk unnamed-chunk-1

## next, a simulated dataset for Argentina & the package mapdata
obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
    regions="Argentina",plot=FALSE)
## objective: 30
## objective: 22
## objective: 16
## objective: 16
## objective: 14
## objective: 14
## objective: 14
## objective: 14
## objective: 14
## objective: 14
## objective: 12
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 12 tips and 11 internal nodes.
## 
## (2) A geographic map with range:
##      -55.05N, -21.79N
##      -73.58W, -53.65W.
## 
## (3) A table containing 12 geographic coordinates.
plot(obj,direction="rightwards",ftype="i",colors=c("blue","white"),
    pts=FALSE,psize=1.2)

plot of chunk unnamed-chunk-2

## finally, another simulated dataset, this time for the island
## of Puerto Rico
obj<-phylo.to.map(PuertoRico.tree,PuertoRico,database="worldHires",
    regions="Puerto Rico",plot=FALSE,rotate=FALSE)
obj
## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 20 tips and 19 internal nodes.
## 
## (2) A geographic map with range:
##      17.92N, 18.52N
##      -67.94W, -65.27W.
## 
## (3) A table containing 20 geographic coordinates.
plot(obj,ftype="off",colors="blue",pch=24,psize=1.5)

plot of chunk unnamed-chunk-3

The data were simulated more or less as follows:

World.tree<-pbtree(n=26,scale=100)
World.tree$tip.label<-replicate(Ntip(World.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
lat<-fastBM(World.tree,sig2=10,bounds=c(-90,90))
long<-fastBM(World.tree,sig2=80,bounds=c(-180,180))
World<-cbind(lat,long)

Argentina<-locator(n=12) ## with world map plotted
Argentina<-cbind(Argentina$y,Argentina$x)
colnames(Argentina)<-c("lat","long")
Argentina.tree<-pbtree(n=nrow(Argentina),scale=100)
Argentina.tree$tip.label<-replicate(Ntip(Argentina.tree),
    paste(sample(LETTERS,1),".",
    paste(sample(letters,round(runif(n=1,min=3,max=10))),
    collapse=""),
    sep=""))
rownames(Argentina)<-Argentina.tree$tip.label

PuertoRico<-locator(n=20) ## with world map plotted
PuertoRico<-cbind(PuertoRico$y,PuertoRico$x)
colnames(PuertoRico)<-c("lat","long")
PuertoRico.tree<-pbtree(n=20,scale=1)
rownames(PuertoRico)<-PuertoRico.tree$tip.label

3 comments:

  1. Thanks for the best blog.it was very useful for me.keep sharing such ideas in the future as well.this is actually what i was looking for,and i am glad to came here. Great to know some new information. Keep writing and sharing more informative posts :)
    essay writing service reviews

    ReplyDelete
  2. Hi Liam, the package is very cool and I could also modify your code with the information provided in your blog here. Thank you. However, I'm struggeling in implementing projections in the map function call, e.g., with mapproj. Have you thought about how to use phylo.to.map() and plot.phylo.to.map() with map projections in mapproj or other packages used for drawing maps.

    ReplyDelete

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.