library(iL04)
library(ape)

data(PA)
data(PA.coo)
data(PA.map)

attr(PA, "Labels") <- sprintf("%s %s", PA.coo[labels(PA),3], labels(PA))
row.names(PA.coo) <- sprintf("%s %s", PA.coo[,3], row.names(PA.coo))

# one of:  nj  bionj  fastme.bal  fastme.ols
t <- nj(PA) 

mds <- cmdscale(cophenetic(t), 3)

asp <- mean.asp(PA.map[, 2])
mdsmap(PA.coo[, 1:2], mds, labels = PA.coo$idx, main = "Pennsylvania, USA", asp = asp, map = PA.map)

dev.set(1)

for (i in 1:3) mds[,i] <- (mds[,i] - min(mds[,i])) / (max(mds[,i]) - min(mds[,i]))
cols <- rgb(mds[,1], mds[,2], mds[,3])

# for help, see: ?plot.phylo
plot(t, type="unrooted", no.margin=TRUE, lab4ut="azial", tip.color=cols, edge.color="#A0A0ff")