Projecting a phylogenetic tree onto a map with multiple geographic points per taxon

In response to a recent user request I recently added the capability to map more than one geographic location per taxon to the phytools plotting function phylo.to.map. To see the code click here. A phytools version containing this update can already be installed via GitHub using devtools.

In the following demo I'll use some fictional data to illustrate how it can be used - first on a world map projecting from the top, then on a single-country map projecting from the left:

library(phytools)
packageVersion("phytools")

## [1] '0.6.71'

library(mapdata)
## here are my data
World.tree

## 
## Phylogenetic tree with 26 tips and 25 internal nodes.
## 
## Tip labels:
##  B.ijyqhp, J.wjbgkp, I.zgm, U.otngy, E.bgzoej, L.cyl, ...
## 
## Rooted; includes branch lengths.

head(World,20)

##                      lat       long
## B.ijyqhp    -53.75945233   89.82698
## B.ijyqhp    -49.27526159   92.08920
## B.iownj      -1.78488475 -162.17161
## B.iownj      -5.28332361 -160.89541
## B.iownj      -1.15704991 -169.82635
## B.xiwt      -18.32161500  -39.76940
## B.xiwt      -21.81437265  -42.49659
## B.xiwt      -16.38054009  -42.98694
## E.bgzoej    -33.49914588  167.24997
## E.bgzoej    -34.83267341  168.81525
## E.bgzoej    -33.50117255  164.69732
## F.lbfqkh      6.04544745  -99.61004
## F.lbfqkh     -0.06604802  -95.04093
## F.lbfqkh      3.37645461 -101.15848
## G.jltpcivk  -16.65924503  -49.84077
## G.jltpcivk  -14.68628460  -48.19405
## G.thjdsean   10.95726052  -73.92540
## G.thjdsean   10.30423882  -72.88593
## G.thjdsean    9.79876619  -70.26079
## G.xqpudsjri   4.33975876   27.89086

obj<-phylo.to.map(World.tree,World,plot=FALSE)

## objective: 82
## objective: 82

## objective: 72
## objective: 72
## objective: 72
## objective: 72

## objective: 70
## objective: 70
## objective: 70
## objective: 70

## objective: 68

## objective: 46
## objective: 46
## objective: 46

## objective: 40
## objective: 40
## objective: 40

## objective: 38
## objective: 38
## objective: 38
## objective: 38
## objective: 38
## objective: 38

## objective: 36

## objective: 34

plot(obj,ftype="i")

If we like, it is very easy to plot the points & linking lines to our species to be shown in different colors. Here's a simple example:

cols<-setNames(sample(rainbow(n=Ntip(World.tree))),
    World.tree$tip.label)
plot(obj,colors=cols,ftype="i",fsize=0.6,cex.points=c(0.7,1.2))

This also should work for sideways plotted trees. For the following it's necessary to install both mapdata (for the higher resolution country map) and viridis (for the color palette that I used).

library(mapdata)
library(viridis)
Argentina.tree

## 
## Phylogenetic tree with 12 tips and 11 internal nodes.
## 
## Tip labels:
##  J.zwilcgka, Z.higpu, X.usrfnti, H.noswc, K.jgquihs, P.xdwzvfj, ...
## 
## Rooted; includes branch lengths.

Argentina

##                   lat      long
## J.zwilcgka  -49.22635 -70.70463
## Z.higpu     -48.56640 -68.50478
## X.usrfnti   -47.35648 -70.04468
## H.noswc     -44.05671 -70.53964
## K.jgquihs   -42.90179 -67.73484
## P.xdwzvfj   -40.31697 -65.47999
## G.xydpt     -36.19226 -67.29487
## N.yinzpjo   -33.44245 -68.22980
## H.hpsycu    -31.46258 -64.05009
## Q.fstnuprx  -29.70271 -60.42035
## B.pnhma     -27.17288 -62.01523
## L.lfbpxwutn -24.25808 -64.98503
## J.zwilcgka  -49.92165 -70.27778
## Z.higpu     -47.80319 -69.42856
## X.usrfnti   -48.19909 -71.25053
## X.usrfnti   -48.46782 -70.90932
## H.noswc     -44.78242 -70.02675
## H.noswc     -43.45621 -71.00061
## K.jgquihs   -44.22061 -65.83883
## K.jgquihs   -42.57006 -67.32508
## G.xydpt     -35.57724 -67.16105
## N.yinzpjo   -33.73010 -68.13654
## H.hpsycu    -30.38485 -61.98099
## Q.fstnuprx  -28.14004 -60.40019
## Q.fstnuprx  -29.56328 -60.08278
## B.pnhma     -26.15989 -62.03782
## B.pnhma     -25.90764 -62.64044
## L.lfbpxwutn -24.41996 -65.55994

obj<-phylo.to.map(Argentina.tree,Argentina,database="worldHires",
    regions="Argentina",plot=FALSE)

## objective: 22
## objective: 22
## objective: 22

## objective: 20

## objective: 18
## objective: 18

## objective: 16
## objective: 16
## objective: 16

## objective: 14

## objective: 12

cols<-setNames(sample(viridis(n=Ntip(Argentina.tree))),
    Argentina.tree$tip.label)
plot(obj,direction="rightwards",colors=cols,ftype="off",cex.points=c(0,1.5),
    pts=FALSE)

This is brand new & involved me working on code that I haven't really touched in a few years so please let me know if you come across any bugs!

For those interested in such things, here's how the data for this example were generated:

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)
for(i in 1:Ntip(World.tree)){
    ni<-sample(0:2,1)
    for(j in 1:ni){ 
        World<-rbind(World,c(World[i,1]+rnorm(n=1,sd=4),
            World[i,2]+rnorm(n=1,sd=4)))
        rownames(World)[nrow(World)]<-rownames(World)[i]
    }
}
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
for(i in 1:Ntip(Argentina.tree)){
    ni<-sample(0:2,1)
    if(ni>0){
        for(j in 1:ni){
            Argentina<-rbind(Argentina,c(Argentina[i,1]+rnorm(n=1),
                Argentina[i,2]+rnorm(n=1)))
            rownames(Argentina)[nrow(Argentina)]<-rownames(Argentina)[i]
        }
    }
}