File:Zagoskin lake pollen 1.svg

3. Pollen diagram, arranged by pollen %
5. "R" script, uses neotoma, rioja, ...
7. 25.10.2021 0000.0007
8. 1. Original code is
9. 1. "Making a pollen diagram from Neotoma", Richard Telford’s Blog
10. 1. https://quantpalaeo.wordpress.com/2019/10/10/making-a-pollen-diagram-from-neotoma/

1. load packages

library("tidyverse") library("neotoma") library("rioja")

1. devtools::install_github("richardjtelford/ggpalaeo") #NEW VERSION

library("ggpalaeo") #Need new version library("stringr")

1. 1. load by sitename. if this fails, try by dataset id!

by_sitename=1 by_dataset_id=0

1. 1. drawing style

draw_by_age=1 draw_by_depth=0 draw_by_depth_without_age=0

1. 1. minimum pollen count %

count_percent_limit1=3

top_margin=0.8

1. sitename1<-"Mammoth Dima site"
2. dataset_id=692

sitename1<-'Zagoskin Lake'

1. sitename1<-'Puyuk Lake'

1. sitename1<-"Plesheevo lake"
2. dataset_id=4396

1. 1. sitename1<-"Lake Nero"
   2. dataset_id=4303

1. sitename1<-"Bugutak"

1. sitename1<-"Bludlivaya River"

1. sitename1<-"Byllatskoye Exposure"

1. sitename1<-"Smorodinovoye Lake"

1. 1. sitename1<-"Lake Billyakh"

1. sitename1<-"Ledovyi Obryv Exposure, Northern Section"
2. dataset_id=1688

1. sitename1<-"Alut Lake"
2. sitename1<-"Mamontovy Khataya IC"

1. sitename1<-'Zagoskin Lake'
2. count_percent_limit1=5

1. sitename1<-'Elikchan 4 Lake'
2. count_percent_limit1=3

1. sitename1<-"Tanneurs"

1. sitename1<-"Marais de la Poupiniere" #NOK

1. sitename1<-'Bois des Gardes'

1. dataset_id=22759
2. sitename1<-'La Chaumette'

1. dataset_id=22759
2. sitename1<-'Lake Gosciaz'

1. dataset_id=42510
2. sitename1<-'Meerfelder Maar'

1. puyuk lake
2. dataset_id=1987 ## if you download by site name, this is overridden
   1. berelekh
   2. dataset_id=3063

1. sitename1<-'Berelyekh River, Indigirka Lowland'

1. sitename1<-'Meerfelder Maar'
2. count_percent_limit1=10

1. sitename1="BER2IND"

count_percent_limit1=5

sitename_x0=tolower(sitename1) sitename_x1<-str_replace(sitename_x0,' ', '_') outfilename1<-paste0("./",sitename_x1,"_pollen_n.svg")

if(by_sitename==1) {

print("By sitename ...")

site1 <- get_site(sitename = paste0(sitename1,'%')) site1_pollen=get_dataset(site1) site1_data0=get_download(site1_pollen)

#str(site1) #str(site1_pollen)

#str(site1_data0)

#str(site1_data01)

#str(site1_data2)

#dataset_indexi=1 #site1_data=site1_data01

site1_data<-site1_data0

#stop(-1)

#print (str(site1))

print (str(site1_pollen))

bambili_raw<-site1_data

#str(bambili_raw)

sq1<-site1_pollen2 sq2<-sq1[2] sq3<-sq2[1] sq4<-sq3[1] sq5<-sq4[1] sq6<-sq51 sq7<-as.data.frame(sq6['dataset.id'])

sxx1<-sq7['dataset.id'] sxx2<-as.character(sxx1) site_id1<-as.integer(sxx2)

site_id1

site_id_char1<-as.character(site_id1)

bambili_raw <- site1_data

}

if(by_dataset_id==1) { print("By dataset ID ...") ds1<-get_dataset(dataset_id)#check correct site print(ds1)

bambili_raw <- get_download(dataset_id)

}

eco_types <- get_table("EcolGroupTypes")

1. 1. by site number:
2. bambili_raw <- get_download(17391)

str(bambili_raw)

tabula0<-bambili_raw1

1. str(tabula0)

cnts<-counts(tabula0) #

meta = tabula0$sample.meta #sample depths/ages

taxa0 = tabula0$taxon.list

taxa = taxa0 %>%

mutate_all(as.character) 

    

bambili <- meta %>% select(age, depth) %>% bind_cols(cnts) %>% pivot_longer(cols = -c("age", "depth"), names_to = "species", values_to = "count") %>% left_join(taxa, by = c("species" = "taxon.name"))

bambili %>% count(variable.element)

eco_types %>% semi_join(bambili, by = c("EcolGroupID" = "ecological.group")) %>% select(EcolGroupID, EcolGroup)

str(eco_types)

bambili %>% filter(ecological.group == "SEED", count > 0) %>% select(species, count) %>% group_by(species) %>% summarise(n = n(), max = max(count))

1. 1. old: filter unwanted groups out

1. bambili = bambili %>%
2. filter(!ecological.group %in% c("AQVP", "UNID", "SEED"))
3. use `%in%` not `==`

1. check count sums

bambili %>% group_by(depth) %>% summarise(s = sum(count)) %>% arrange(s) %>% slice(1:5)

1. calculate percent

bambili = bambili %>%

 group_by(depth) %>% 
 mutate(percent = count/sum(count) * 100) 

1. remove rare taxa

bambili1 = bambili %>%

 group_by(species) %>% 
 filter(
   sum(percent > 0) >= count_percent_limit1, #must be in at least three samples
   max(percent) > count_percent_limit1) #must have a max percent > 3
   

bambili2 <- bambili1 %>% select(age, depth, species, percent) %>% pivot_wider(names_from = "species", values_from = "percent")

bambili_spp <- bambili2 %>% select(-age, -depth) %>% as.data.frame()#

bambili2 <- bambili1 %>% mutate(

1. make ecological.group a factor with TRSH first
2. ecological.group = factor(ecological.group, levels = c("TRSH", "UPHE", "VACR","AQVP", "AQBR","ALGA")),
3. ecological.group = factor(ecological.group, levels = c("TRSH","UPHE","AQVP","VACR","AQBR","FUNGI","ALGA","FUNG","AMOE","TARD")),

ecological.group = factor(ecological.group, levels = c("TRSH", "UPHE", "VACR", "AQBR", "AQVP", "TARD", "ALGA")),

mean_percent = mean(percent)) %>%

1. arrange by ecological.group and mean_percent (largest first)

arrange(ecological.group, desc(mean_percent)) %>% ungroup() %>%

1. make species into a factor so we can perserve the order

mutate(species = factor(species, levels = unique(species)))

1. reshape using tidyr::spread as pivot_wider (currently?) ignores factor order

bambili3 <- bambili2 %>% select(age, depth, species, percent) %>% spread(key = "species", value = "percent")

bambili_spp <- bambili3 %>% select(-age, -depth) %>% as.data.frame(bambili_spp)

1. set up for ecological group colours

ecological_groups <- bambili2 %>%

 distinct(species, ecological.group) %>% 
 pull(ecological.group)
 

1. ecological_colours <- c("red", "green", "orange")

ecological_colours <- c("forestgreen", "orange", "red", "violet", "blue", "darkblue","darkviolet")

1. 1. jn warning

bambili2_spp <- bambili2 %>% select(-age, -depth) %>% as.data.frame(bambili_spp)

1. stop(-1)

bambili_dist <- dist(sqrt(bambili_spp/100))#chord distance clust <- chclust(bambili_dist, method = "coniss")

1. bstick(clust)#five groups

if(draw_by_age==1) {

svg(filename=outfilename1, width=8, height=10, pointsize=16)

1. set up mgp (see ?par)

mgp <- c(2, 0.25, 0) par(tcl = -0.15, mgp = mgp)#shorter axis ticks - see ?par

pt <- strat.plot( title=sitename1, cex.title=1.5, d = bambili_spp,

1. yvar = bambili3$depth,

yvar = bambili3$age, y.rev = TRUE, #reverse direction of y-axis scale.percent = TRUE, #use constant scale for all taxa srt.xlabel = 90, #rotate x-label by 45 degrees cex.xlabel = 0.8, #smaller font mgp = mgp, xRight = 0.98, #right margin xLeft = 0.21, #left margin with space for 2nd axis

yTop = top_margin, #top margin yBottom = 0.1, #bottom margin plot.line=FALSE, plot.poly=TRUE, plot.bar=FALSE, col.line = ecological_colours[ecological_groups],#colours col.poly = ecological_colours[ecological_groups], #colours col.poly.line = ecological_colours[ecological_groups], #colours col.bar = ecological_colours[ecological_groups], #colours lwd.poly=1,

ylabel = "Date yr BP", clust = clust )

1. add zone boundaries

addClustZone(pt, clust = clust, nZone = 5)

1. add a secondary scale

secondary_scale(pt, yvar = bambili3$age,

               yvar2 = bambili3$depth, 
               ylabel2 = "Depth m",
               n = 10)

dev.off()

}

if(draw_by_depth==1) {

svg(filename=outfilename1, width=8, height=10, pointsize=16)

1. set up mgp (see ?par)

mgp <- c(2, 0.25, 0) par(tcl = -0.15, mgp = mgp)#shorter axis ticks - see ?par

pt <- strat.plot( title=sitename1, cex.title=1.5, d = bambili_spp, yvar = bambili3$depth,

1. yvar = bambili3$age,

y.rev = TRUE, #reverse direction of y-axis scale.percent = TRUE, #use constant scale for all taxa srt.xlabel = 90, #rotate x-label by 45 degrees cex.xlabel = 0.8, #smaller font mgp = mgp, xRight = 0.98, #right margin xLeft = 0.21, #left margin with space for 2nd axis

yTop = top_margin, #top margin yBottom = 0.1, #bottom margin plot.line=FALSE, plot.poly=TRUE, plot.bar=FALSE, col.line = ecological_colours[ecological_groups],#colours col.poly = ecological_colours[ecological_groups], #colours col.poly.line = ecological_colours[ecological_groups], #colours col.bar = ecological_colours[ecological_groups], #colours lwd.poly=1,

ylabel = "Depth m", clust = clust )

1. add zone boundaries

addClustZone(pt, clust = clust, nZone = 5)

1. add a secondary scale

secondary_scale(pt, yvar = bambili3$depth,

               yvar2 = bambili3$age, 
               ylabel2 = "Age BP",
               n = 10)

dev.off()

}

if(draw_by_depth_without_age==1) {

svg(filename=outfilename1, width=8, height=10, pointsize=16)

1. set up mgp (see ?par)

mgp <- c(2, 0.25, 0) par(tcl = -0.15, mgp = mgp)#shorter axis ticks - see ?par

pt <- strat.plot( title=sitename1, cex.title=1.5, d = bambili_spp, yvar = bambili3$depth,

1. yvar = bambili3$age,

y.rev = TRUE, #reverse direction of y-axis scale.percent = TRUE, #use constant scale for all taxa srt.xlabel = 90, #rotate x-label by 45 degrees cex.xlabel = 0.8, #smaller font mgp = mgp, xRight = 0.98, #right margin xLeft = 0.21, #left margin with space for 2nd axis

yTop = top_margin, #top margin yBottom = 0.1, #bottom margin plot.line=FALSE, plot.poly=TRUE, plot.bar=FALSE, col.line = ecological_colours[ecological_groups],#colours col.poly = ecological_colours[ecological_groups], #colours col.poly.line = ecological_colours[ecological_groups], #colours col.bar = ecological_colours[ecological_groups], #colours lwd.poly=1,

ylabel = "Depth m", clust = clust )

1. add zone boundaries

addClustZone(pt, clust = clust, nZone = 5)

1. add a secondary scale
2. secondary_scale(pt, yvar = bambili3$depth,
3. yvar2 = bambili3$age,
4. ylabel2 = "Age BP",
5. n = 10)

dev.off()

}

print(".")