ggplot2 pollen diagram (and dendrogram)

In R, excellent packages provide ways to plot pollen or more generally stratigraphic diagrams. The two main tools come from the rioja package with “strat.plot” and the analogue package with “Stratiplot”. Althought those two functions are very comprehensive (you can include a dendrogram, pollen zones, etc.), easy to use, and highly customizable; I was still wondering if there is a way in R to plot a simple pollen diagram using only general plot syntax an preferably ggplot2. I came up with this simple solution that involve only ggplot2 syntax. Probably there is a possibility to add a dendrogram with grid.arrange and to overlay additionnal information but right now this approach seems sufficient to cover my very basic plotting needs:

library(ggplot2)
library(grid)
library(neotoma)
library(analogue)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.4-1
## analogue version 0.17-0

First extract data from Neotoma using neotoma package, we will have a look at “Le Grand Etang de Suze-La-Rousse” in Provençal Drome by: “Argant, J. 1990. Climat et environnement au Quaternaire dans le bassin du Rhône d’après les données palynologiques. Documents du Laboratoire de Géologie de Lyon 3:1-199.”

suze <- get_site(sitename = 'Le Grand Etang%')
suze_pollen=get_dataset(suze)
suze_data=get_download(suze_pollen)

# head(suze_data[[1]]$counts)

core.pct <- data.frame(tran(suze_data[[1]]$counts, method = 'percent'))

age <- suze_data[[1]]$sample.meta$age
core.pct <- chooseTaxa(core.pct, max.abun = 10) 

Stratiplot(age ~ ., core.pct, sort = 'wa', type = 'poly',
 ylab ='Years Before Present')
           

Now we create a simple data.frame with information for plotting a pollen diagram in ggplot2:

df=data.frame(yr=rep(age,ncol(core.pct)),
 per=as.vector(as.matrix(core.pct)),
 taxa=as.factor(rep(colnames(core.pct),each=nrow(core.pct))))

We define a theme without gridlines, etc. and elements that will decrease readability of the diagram

theme_new <- theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), # remove grids
 panel.background = element_blank(), axis.line = element_line(colour = "black"),
 strip.text.x = element_text(size=10, angle=90, vjust=0), # Taxa names 
 strip.background = element_blank(),
 strip.text.y = element_text(angle = 0),
 legend.position="none",panel.border = element_blank(),
 axis.text.x=element_text(angle=45,hjust=1)) # Axis tick label angle

And finally the pollen diagram:

ggplot(df)+
 geom_line(aes(yr,per))+
 geom_area(aes(yr,per))+
 scale_x_reverse(breaks =seq(0,100000,1000))+
 scale_y_continuous(breaks =seq(0,100,10))+
 xlab("Age (cal. BP)")+ylab("%")+
 coord_flip()+
 theme_new+
 facet_grid(~df$taxa,scales = "free", space = "free")

Add an exageration line (the quick and dirty way):

df$exag=df$per*10
id=unique(df$taxa)
for(i in 1:length(id)){
 uplim=which(df$exag[df$taxa==id[i]]>max(df$per[df$taxa==id[i]]))
 df$exag[df$taxa==id[i]][uplim]<-max(df$per[df$taxa==id[i]])
} 

diag=ggplot(df)+
 geom_line(aes(yr,per))+
 geom_area(aes(yr,exag,fill=taxa))+
 geom_area(aes(yr,per))+
 scale_x_reverse(breaks =seq(0,1e8,1000))+
 scale_y_continuous(breaks =seq(0,100,10))+
 xlab("Age (cal. BP)")+ylab("%")+
 coord_flip(xlim=c(3000,14000))+
 theme_new+
 facet_grid(~df$taxa,scales = "free", space = "free")

diag

 

Because only ggplot2 syntax is used you can play around and add geom_point, geom_bar layouts etc. and customize your diagram the way you want with the theme() function.

UPDATE: Adding a dendrogram is not straightforward but definitively doable. We will use a combinason of rioja and ggdendro package such as:

diss <- dist(sqrt(core.pct/100)^2)
clust <- chclust(diss, method="coniss")
# bstick(clust) # optionnally look at number of optimal zones

Now we will use ggdendro and modify its output in order to inject ages that will replace sample numbers:

dendro <- as.dendrogram(clust)
ddata <- dendro_data(dendro, type="rectangle")

## MOD the segment data frame for age instead of value order -----------------------
ddata$segments->yo
yo$xx=NA
yo$xxend=NA

for(i in 1:length(age)){
 yo$xx[which(yo$x == i)]=age[i]
 yo$xxend[which(yo$xend == i)]=age[i]
 da=age[i+1]-age[i]
 yo$xx[yo$x>i & yo$x<i+1]<-age[i]+(yo$x[yo$x>i & yo$x<i+1]-i)*da
 yo$xxend[yo$xend>i & yo$xend<i+1]<-age[i]+(yo$xend[yo$xend>i & yo$xend<i+1]-i)*da
} # Please dont ask....
ddata$segments$x<-yo$xx
ddata$segments$xend<-yo$xxend

# -----------------------------------------------------------------------------------

We finally can plot the dendrogram along with the diagram with gridExtra

dendro=ggplot(segment(ddata)) +
 geom_segment(aes(x=x, y=y, xend=xend, yend=yend)) +
 coord_flip(xlim=c(3000,14000)) +
 ylab("Distance")+
 scale_x_reverse(breaks =seq(0,1e8,1000))+
 theme_new+theme(axis.title.y=element_blank(),
 axis.text.y=element_blank(),
 axis.ticks.y=element_blank(),
 axis.line.y = element_blank())


g1 = ggplotGrob(diag)
g2 = ggplotGrob(dendro)
g2$heights[5:7]<-g1$heights[6:8] # Important to align plots!
grid.arrange(g1, g2, widths =c(1,0.3),ncol=2,newpage = TRUE)

You can follow the overall logic to add any type of additional plot as long as its a ggplot2 object!