###############################################################################################################################################
###	R script for plotting PBDB collections through time																						###
###	Written by V. Fischer in April–July 2020																								###
###	Associated paper: Fischer, Weis & Thuy. In Review. Refining the marine reptile turnover at the Early-Middle Jurassic transition. PeerJ.	###
###############################################################################################################################################

library(geoscale)

PBDB_raw <- read.csv("SUPP_pbdb_data.csv", skip=18, header=TRUE)
PBDB <- PBDB_raw[PBDB_raw$max_ma<202 & PBDB_raw$min_ma>144, ]
PBDB$time_bins <- as.factor(as.character(PBDB$time_bins))
count <- table(PBDB$time_bins)[-1]
time_raw <- read.csv("SUPP_jurassic_time.csv", header=TRUE)
mid_stages <- rowMeans(time_raw[, 2:3]); names(mid_stages) <- time_raw[, 1]
collections <- as.data.frame(cbind(mid_stages, count=count[names(mid_stages)]))

pdf(file="Collections.pdf", width=9, height=4)
geoscalePlot(collections$mid_stages, collections$count, type="l", pch=19, col='#01abe9', lwd=2.5, data.lim=c(0, max(collections$count)+10), age.lim=c(201, 145), erotate=0, arotate=0, cex.pt=1, cex.age=0.7, cex.ts=0.5, label="Number of PBDB collections")
rect(177, -10, 166.1, max(collections$count)+10, col='#D3D3D350', border=NA)
text(x=176.4, y=max(collections$count), labels="Turnover window \nuncertainty", adj=0)
dev.off()