Monday, July 13, 2015

Adding gridlines to a phylo.to.map plot

A phytools blog reader asks:

“Is it possible to add latitude and longitude tick marks with this function phylo.to.map?”

The answer is - yes. Well, at least we can adjust the margins to include x & y axes (which are on the scale of lat. & long.), and (if we want to) we can even overlay lines of latitude & longitude on the plot. Here's how:

## first simulate a tree & some data:
library(phytools)
tree<-pbtree(n=26,scale=100,tip.label=LETTERS)
lat<-fastBM(tree,sig2=10,bounds=c(-36,-16),a=-20)
long<-fastBM(tree,sig2=80,bounds=c(115,150),a=130)
loc<-cbind(lat,long)
loc
##         lat     long
## A -30.00298 136.8887
## B -29.52615 139.2913
## C -21.94540 117.6450
## D -21.44360 137.3063
## E -30.69070 129.6878
## F -27.91373 138.7563
## G -21.24415 119.6728
## H -19.85676 138.0375
## I -29.70913 143.8215
## J -32.41438 130.5098
## K -34.28836 127.8772
## L -32.97290 135.9813
## M -26.05441 145.2799
## N -23.89187 124.5674
## O -17.84134 140.5776
## P -28.67361 138.5719
## Q -16.77305 128.0718
## R -26.44700 140.7354
## S -24.02262 139.6478
## T -35.53304 136.3906
## U -30.89431 139.9094
## V -19.50736 149.8240
## W -21.92831 149.1358
## X -18.30595 138.6660
## Y -19.17043 140.1686
## Z -22.59701 132.3479

Now, with these (made-up) data in hand, let's create our plot & customize it:

## now plot projection
phylo.to.map(tree,loc,ylim=c(-40,-10), xlim=c(110,155),
    mar=c(5.1,4.1,2.1,2.1))
## objective: 180
## objective: 172
## objective: 172
## objective: 172
## objective: 168
## objective: 168
## objective: 168
## objective: 154
## objective: 140
## objective: 140
## objective: 140
## objective: 140
## objective: 140
## objective: 138
## objective: 138
## objective: 134
## objective: 134
## objective: 134
## objective: 134
## objective: 134
## objective: 134
## objective: 126
## objective: 126
## objective: 126
## objective: 126
## create the axes
long.lines<-seq(110,160,by=10)
lat.lines<-seq(-40,-10,by=10)
axis(1,at=long.lines)
axis(2,at=seq(-40,-10,by=5))
## add gridlines
for(i in 1:length(lat.lines))
    lines(c(par()$usr[1],160),rep(lat.lines[i],2),lty="dotted",
        col="#00000040")
for(i in 1:length(long.lines))
    lines(rep(long.lines[i],2),c(par()$usr[3],-10),lty="dotted",
        col="#00000040")
lines(c(par()$usr[1],160),rep(-23,43333,2),lwd=2,col="#00000080")
text(105,-23.5,"Tropic of Capricorn",cex=0.7,pos=4)
title(xlab="longitude")
title(ylab="latitude                         ")

plot of chunk unnamed-chunk-2

Well, you get the idea.

1 comment:

  1. I'm trying to plot the graph through data without simulations. At the end of the plot generation process an error occurs. (Error in xy.coords (x, y) 'x' and 'y' lengths differ).
    How can I do this plot?
    Follow the step by step I'm running:

    > Require (phytools)
    > Library (rworldmap)
    > Tree <- read.tree ("/ home / cleysinho / begomodw / view / file / FastTreeMPTueJul2116: 34: 222015.tree")
    > Source ("/ home / cleysinho / begomodw / view / file / phylo.to.map.R")
    > Lat_log = read.table ("FASTATueJul2117_41_462015.txt_lat_log.txt", sep = "\ t")
    > Data <- (lat_log [2: 3])
    > Rownames (date) <- lat_log [, 1]
    > COLNAMES (date) <- c ("lat", "long")
    > Date
    lat long
    FJ538207 -38.41610 -63.61667
    KJ742419 -38.41610 -63.61667
    JX863082 -23.89924 -51.22822
    JX863081 -23.74301 -51.31422
    KC683374 -23.44250 -58.44383
    AY090558 -20.75640 -42.87458
    NC_004639 -14.23500 -51.92528
    JX871385 -14.23500 -51.92528
    JX871384 -14.23500 -51.92528
    JX871388 -14.23500 -51.92528
    EU710750 -20.85773 -42.80155
    KC706538 -20.81462 -42.87776
    EU710749 -22.43139 -43.42895
    JN848772 6.42375 -66.58973
    KJ174330 18.73569 -70.16265
    KC004097 -14.23500 -51.92528
    JF694482 -12.57974 -41.70073

    > Phylo.to.map (tree, date)
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    objective: 0
    Error in xy.coords (x, y) 'x' and 'y' lengths differ

    Thanks, Cleydson

    ReplyDelete