Updates to phytools function phylo.to.map for projecting a phylogeny onto a geographic map

The phytools function phylo.to.map is very popular but objectively kind of embarrassing.

I wrote this function many years ago & have hardly touched the code base since. The plotting method associated with the function is very inflexible & users often struggle with issues such as the map improperly overlapping the phylogeny and other things.

I just pushed a number of small improvements to the phytools GitHub page. Primarily, the function now (by default) limits the map area to the range of the tip species geographic coordinates & also clips map regions that don’t fall within the map area.

To obtain these updates, users should install phytools from GitHub.

library(phytools)
packageVersion("phytools")

## [1] '1.9.23'

To see how it works, we can run a couple of quick demos.

First, let’s use a simple dataset (mammal.geog) based on iNaturalist citizen science records. I’ll restrict the distribution of points to North America.

data(mammal.tree)
data(mammal.geog)

This is just to geographically subsample to North America.

na_mammals<-mammal.geog[
  intersect(which(mammal.geog[,"lat"]>0),
    which(mammal.geog[,"long"]<(-20))),]
na_tree<-keep.tip(mammal.tree,
  unique(rownames(na_mammals)))

Now let’s create our phylo.to.map object.

mammal.phymap<-phylo.to.map(na_tree,na_mammals,
  plot=FALSE)

## objective: 92
## objective: 88
## objective: 88
## objective: 88
## objective: 88
## objective: 86
## objective: 86
## objective: 84
## objective: 82
## objective: 76
## objective: 74
## objective: 74
## objective: 70
## objective: 68
## objective: 68
## objective: 68

If we plot this using the new default settings this is what we obtain.

plot(mammal.phymap,fsize=0.6,ftype="i")

plot(mammal.phymap,xlim=c(-175,-35),ylim=c(0,90),
  fsize=0.6,ftype="i",asp=0.8,lty="solid",
  map.bg="darkgreen",map.fill="lightblue",
  lwd=3,pts=FALSE,colors="darkorange",
  cex.points=1.4,delimit_map=TRUE)

OK. Let’s try a different example – this one is for Galapagos tortoises using data obtained from Poulakakis et al. (2020).

For this example, I’ll use a higher resolution map from the package mapdata.

library(mapdata)

data(tortoise.tree)
data(tortoise.geog)
tortoise.phymap<-phylo.to.map(tortoise.tree,
  tortoise.geog,plot=FALSE,direction="rightwards",
  database="worldHires")

## objective: 64

## objective: 52
## objective: 52
## objective: 52
## objective: 52
## objective: 52
## objective: 52
## objective: 52
## objective: 52
## objective: 46
## objective: 46
## objective: 46
## objective: 44
## objective: 44

tortoise.phymap

## Object of class "phylo.to.map" containing:
## 
## (1) A phylogenetic tree with 15 tips and 14 internal nodes.
## 
## (2) A geographic map with range:
##      -85.47N, 83.62N
##      -180W, 180W.
## 
## (3) A table containing 15 geographic coordinates (may include
##     more than one set per species).
## 
## If optimized, tree nodes have been rotated to maximize alignment
## with the map when the tree is plotted in a rightwards direction.

plot(tortoise.phymap,direction="rightwards",
  xlim=c(-92.5,-89),ylim=c(-2,1),pts=FALSE,
  ftype="i",fsize=0.8,map.bg="lightgreen",
  lty="solid",colors="darkblue",
  delimit_map=TRUE)

Lastly, we can still project the tree directly onto the map. Let’s try it here.

plot(tortoise.phymap,type="direct",
  xlim=c(-92.7,-88.5),ylim=c(-1.8,0.8),pts=FALSE,
  ftype="i",fsize=0.8,map.bg="lightgreen",
  map.fill="lightblue",lty="solid",
  colors="darkblue",delimit_map=TRUE)

I mean, this isn’t exactly cutting edge stuff – but I’m happy to be able to offer some small improvements for something that is so widely used.