AshplotR

1.0 Introduction

This script can be used to create bi-plots and Harker plots of all tephra major elemental data. It includes: tips for flitering data, interactive exploration of the data, more complex plots, using kernel densities to delimit chemical envelopes, and some discussion on identifying useful elements for discrimination via multivariate analysis. Guidance for saving plots as .svg or .pdf files is also provided. All the code used to produce this document is downloadable as an RMarkdown file to the right of this page and should appear for each specific plot generated. To run this you will require r and RStudio. We are happy for you to use, share and improve our code, but would ask that you cite this as a point of origin for the code.

Load packages

We first load in packages required to generate the plots and manipulate the data. You may need to install these packages in r or update them to work correctly.

library(ggplot2)
library(factoextra)
library(readr)
library(psych)
library(compositions)
library(gganimate)
library(cowplot)
library(tidyverse)
library(gt)
library(plotly)
library(svglite)
library(png)
library(here)
library(ggdensity)
library(ggrepel)
library(viridis)

2.0 Loading your data

Data should be placed within a ‘Data’ folder created in your project directory. As this is a Markdown script the ‘here’ package is used to ensure that project links do not break when moving scripts around, however you may not need to do this if you have a preferred data storage approach. As long as you start a new R project and drop this script into it then set up ‘Data’ and ‘Script’ folders then it will work. Data files should be in .csv format and contain the group of the tephra as column 1 under the title “id” (see console output below for format). Loaded here are some data for core ODP-980 in the North Atlantic and a terrestrial site from the UK called MT. These data were published by Candy et al. (2021) and the paper can be found here.

chem <- here::here("Data","ODP_norm.csv") %>% read.csv()

2.1 Check your data has loaded correctly

We can first just inspect the data in the console using the ‘head’ command, and then summarise it using the describe function of the psych package. These commands produce a couple of basic tables which should allow you to check whether you have correctly loaded the data you wish to analyse.

head(chem)

##            id     SiO2       TiO2    Al2O3       FeO        MnO        MgO
## 1 50.65-50.69 52.79290 3.62976406 13.17806 13.984674 0.26214963 3.98265779
## 2 50.65-50.69 71.85240 0.30401510 12.67428  4.665059 0.16773247 0.06289968
## 3 50.65-50.69 76.69421 0.08264463 11.67355  1.508264 0.05165289 0.02066116
## 4 50.65-50.69 76.76455 0.15455951 11.89078  2.071097 0.04121587 0.05151983
## 5 50.65-50.69 49.75570 3.87825733 12.85627 15.054967 0.24429967 4.80456026
## 6 50.65-50.69 71.67439 0.31942298 12.84905  4.420402 0.16486347 0.06182380
##         CaO     Na2O       K2O       P2O5 Total IRD
## 1 8.2173825 2.520669 0.7864489 0.65537407   100 YES
## 2 1.1741273 5.493238 3.5852815 0.02096656   100 YES
## 3 0.3099174 3.584711 6.0743802 0.00000000   100 YES
## 4 0.5461103 4.399794 4.0803709 0.01030397   100 YES
## 5 9.1307003 3.145358 0.6616450 0.46824104   100 YES
## 6 1.2570840 5.461103 3.7712519 0.02060793   100 YES

describe(chem)

##       vars   n   mean   sd median trimmed  mad    min    max range  skew
## id*      1 240   8.44 4.44   9.00    8.58 4.45   1.00  15.00 14.00 -0.25
## SiO2     2 240  66.45 9.91  72.00   67.31 2.84  48.10  78.07 29.98 -0.83
## TiO2     3 240   1.10 1.19   0.45    0.90 0.27   0.08   3.88  3.80  1.15
## Al2O3    4 240  11.86 2.08  12.70   12.03 1.11   7.07  15.62  8.55 -0.80
## FeO      5 240   7.69 4.70   7.22    7.37 5.46   1.42  17.13 15.71  0.56
## MnO      6 240   0.21 0.08   0.23    0.21 0.07   0.03   0.60  0.57 -0.06
## MgO      7 240   1.43 2.14   0.12    1.09 0.19  -0.03   7.54  7.58  1.08
## CaO      8 240   3.27 3.79   1.00    2.78 0.93   0.25  12.12 11.87  0.93
## Na2O     9 240   4.80 1.28   5.18    4.86 1.32   2.16   6.94  4.78 -0.39
## K2O     10 240   3.03 1.63   3.75    3.13 0.97   0.20   6.83  6.62 -0.57
## P2O5    11 240   0.16 0.25   0.03    0.11 0.03  -0.01   1.37  1.38  2.11
## Total   12 240 100.00 0.00 100.00  100.00 0.00 100.00 100.00  0.00   NaN
## IRD*    13 240   1.15 0.36   1.00    1.06 0.00   1.00   2.00  1.00  1.95
##       kurtosis   se
## id*      -1.16 0.29
## SiO2     -1.10 0.64
## TiO2     -0.26 0.08
## Al2O3    -0.81 0.13
## FeO      -1.11 0.30
## MnO       1.43 0.01
## MgO      -0.54 0.14
## CaO      -0.95 0.24
## Na2O     -1.24 0.08
## K2O      -1.14 0.11
## P2O5      5.06 0.02
## Total      NaN 0.00
## IRD*      1.80 0.02

2.2 Better tables in r.

The main problem with the console output from r is that produces tables that are not pretty to view. This can largely be addressed by the ‘gt’ package which produces far nicer output in an html format. The following script produces a summary table of mean values and number of analyses in your data.

layer	n	SiO2.	sd1	TiO2.	sd2	Al2O3.	sd3	FeOt	sd4	MnO.	sd5	MgO.	sd6	CaO.	sd7	Na2O.	sd8	K2O.	sd9	P2O5.	sd10	Total.	sd11
50.65-50.69	22	65.98	11.52	1.29	1.49	12.24	1.51	7.40	5.36	0.17	0.09	1.77	2.46	3.82	4.12	4.20	1.20	2.94	1.99	0.19	0.24	100	0
50.93-50.97	14	72.64	2.53	0.35	0.15	11.41	2.58	4.83	2.32	0.19	0.09	0.09	0.13	0.70	0.35	5.59	0.68	4.17	0.60	0.03	0.03	100	0
51.09-51.13	15	58.55	11.65	2.35	1.69	12.22	1.45	11.81	6.20	0.23	0.06	2.85	2.32	5.88	4.33	3.86	1.39	1.97	1.88	0.27	0.21	100	0
51.25-51.29	7	73.67	1.78	0.29	0.11	12.83	0.53	3.05	0.80	0.12	0.05	0.15	0.11	1.18	0.47	4.77	0.73	3.91	0.78	0.03	0.02	100	0
52.91-52.95	7	57.96	9.78	2.31	1.32	12.62	1.74	11.53	4.08	0.25	0.04	2.75	1.89	5.96	3.59	4.18	1.16	1.87	1.56	0.56	0.57	100	0
53.13-53.17	14	71.61	3.99	0.50	0.37	11.88	2.20	5.25	2.54	0.19	0.06	0.30	0.53	1.38	1.53	5.19	0.67	3.61	0.88	0.10	0.22	100	0
54.01-54.05	16	66.29	8.41	1.04	0.97	11.83	2.89	7.88	3.42	0.23	0.06	1.26	1.93	2.94	3.31	5.18	1.19	3.18	1.46	0.16	0.18	100	0
54.13-54.21	18	62.43	10.16	1.38	1.01	12.42	1.97	9.21	4.20	0.22	0.06	2.43	2.49	5.09	4.35	4.38	1.39	2.25	1.66	0.19	0.18	100	0
54.37-54.41	15	61.27	10.76	1.35	0.99	12.42	1.44	10.05	5.19	0.22	0.05	2.79	2.71	5.53	4.57	4.11	1.51	2.12	1.97	0.13	0.10	100	0
54.53-54.57	19	66.01	10.21	1.10	1.12	11.35	2.64	7.90	3.78	0.23	0.06	1.63	2.42	3.36	4.18	5.15	1.38	3.12	1.59	0.16	0.21	100	0
54.57-54.61	20	70.90	5.51	0.50	0.56	11.48	2.30	5.77	3.13	0.22	0.11	0.33	1.18	1.43	2.18	5.55	0.90	3.78	1.04	0.04	0.08	100	0
54.79-54.83	21	68.39	10.41	0.95	1.21	11.79	1.73	6.56	4.87	0.18	0.09	1.25	2.22	2.80	3.86	4.80	1.19	3.19	1.58	0.10	0.14	100	0
55.11-55.15	14	67.78	9.43	1.02	1.33	11.86	1.83	7.05	4.81	0.21	0.07	1.02	1.89	2.57	3.43	5.08	1.19	3.32	1.53	0.11	0.18	100	0
55.31-55.35	21	62.97	10.71	1.42	1.18	11.57	2.32	10.05	4.77	0.25	0.07	1.97	2.13	4.44	4.05	4.51	1.39	2.45	1.76	0.38	0.41	100	0
55.79-55.87	17	70.24	7.86	0.75	1.16	11.16	2.24	6.28	4.05	0.20	0.08	0.63	1.59	1.87	2.86	5.20	1.14	3.59	1.29	0.09	0.17	100	0

2.3 Make TAS calculations

To calculate for the TAS plot the data require normalisation of the data. The data used here has been normalised already so this step is largely redundant, but this is not the case for all tephra data and you might prefer to utilise unormalised data for your analyses, especially if you are interested in the effects of analytical totals on your data. This code does not modify your orignal data so can be safely run for both normalised and unormalised data. It is required for the TAS plot so should be run regardless of the state of your data.

TAS_x <- (100/ chem [["Total"]]) * chem [["SiO2"]]
TAS_y <- (100/ chem [["Total"]]) * (chem [["K2O"]] + chem [["Na2O"]])

chem <- cbind(chem, TAS_x, TAS_y)

3.0 Filtering data

It is possible your datasheet contains multiple types of tephra chemistry which would benefit from filtering before plotting. I suggest this is carried out initially by using silica values as a divisor. Rhyolites >68% silica would be a case for this. However, it might be worth plotting everything first before then selecting a filter to ensure vital data is not being excluded. Here we continue with all data being plotted before filtering later in the script when making comparisons.

#chem <- filter(chem, SiO2 >68) ## code is commented order to retain all data at this stage.

4.0 Plotting your data

First we generate a series of plots from all possible pairs of elements.

p1 <- ggplot(chem, aes(x = TAS_x, y = TAS_y, fill = id, shape = id)) # this uses the TAS calculations made earlier in step 2.3.
p2 <- ggplot(chem, aes(x = SiO2, y = TiO2, fill = id, shape = id))
p3 <- ggplot(chem, aes(x = SiO2, y = Al2O3, fill = id, shape = id))
p4 <- ggplot(chem, aes(x = SiO2, y = FeO,  fill = id, shape = id))
p5 <- ggplot(chem, aes(x = SiO2, y = MnO,  fill = id, shape = id))
p6 <- ggplot(chem, aes(x = SiO2, y = MgO,  fill = id, shape = id))
p7 <- ggplot(chem, aes(x = SiO2, y = CaO,  fill = id, shape = id))
p8 <- ggplot(chem, aes(x = SiO2, y = Na2O,  fill = id, shape = id))
p9 <- ggplot(chem, aes(x = SiO2, y = K2O,  fill = id, shape = id))
p47 <- ggplot(chem, aes(x = SiO2, y = K2O, fill = id, shape = id))
p10 <- ggplot(chem, aes(x = SiO2, y = P2O5,  fill = id, shape = id))
p11 <- ggplot(chem, aes(x = TiO2, y = Al2O3,  fill = id, shape = id))
p12 <- ggplot(chem, aes(x = TiO2, y = FeO,  fill = id, shape = id))
p13 <- ggplot(chem, aes(x = TiO2, y = MnO,  fill = id, shape = id))
p14 <- ggplot(chem, aes(x = TiO2, y = MgO,  fill = id, shape = id))
p15 <- ggplot(chem, aes(x = TiO2, y = CaO,  fill = id, shape = id))
p16 <- ggplot(chem, aes(x = TiO2, y = Na2O,  fill = id, shape = id))
p17 <- ggplot(chem, aes(x = TiO2, y = K2O,  fill = id, shape = id))
p18 <- ggplot(chem, aes(x = TiO2, y = P2O5,  fill = id, shape = id))
p19 <- ggplot(chem, aes(x = Al2O3, y = FeO,  fill = id, shape = id))
p20 <- ggplot(chem, aes(x = Al2O3, y = MnO,  fill = id, shape = id))
p21 <- ggplot(chem, aes(x = Al2O3, y = MgO,  fill = id, shape = id))
p22 <- ggplot(chem, aes(x = Al2O3, y = CaO,  fill = id, shape = id))
p23 <- ggplot(chem, aes(x = Al2O3, y = Na2O,  fill = id, shape = id))
p24 <- ggplot(chem, aes(x = Al2O3, y = K2O,  fill = id, shape = id))
p25 <- ggplot(chem, aes(x = Al2O3, y = P2O5,  fill = id, shape = id))
p26 <- ggplot(chem, aes(x = FeO, y = MnO,  fill = id, shape = id))
p27 <- ggplot(chem, aes(x = FeO, y = MgO,  fill = id, shape = id))
p28 <- ggplot(chem, aes(x = FeO, y = CaO,  fill = id, shape = id))
p29 <- ggplot(chem, aes(x = FeO, y = Na2O,  fill = id, shape = id))
p30 <- ggplot(chem, aes(x = FeO, y = K2O,  fill = id, shape = id))
p31 <- ggplot(chem, aes(x = FeO, y = P2O5,  fill = id, shape = id))
p32 <- ggplot(chem, aes(x = MnO, y = MgO,  fill = id, shape = id))
p33 <- ggplot(chem, aes(x = MnO, y = CaO,  fill = id, shape = id))
p34 <- ggplot(chem, aes(x = MnO, y = Na2O,  fill = id, shape = id))
p35 <- ggplot(chem, aes(x = MnO, y = K2O,  fill = id, shape = id))
p36 <- ggplot(chem, aes(x = MnO, y = P2O5,  fill = id, shape = id))
p37 <- ggplot(chem, aes(x = MgO, y = CaO,  fill = id, shape = id))
p38 <- ggplot(chem, aes(x = MgO, y = Na2O,  fill = id, shape = id))
p39 <- ggplot(chem, aes(x = MgO, y = K2O,  fill = id, shape = id))
p40 <- ggplot(chem, aes(x = MgO, y = P2O5,  fill = id, shape = id))
p41 <- ggplot(chem, aes(x = CaO, y = Na2O,  fill = id, shape = id))
p42 <- ggplot(chem, aes(x = CaO, y = K2O,  fill = id, shape = id))
p43 <- ggplot(chem, aes(x = CaO, y = P2O5,  fill = id, shape = id))
p44 <- ggplot(chem, aes(x = Na2O, y = K2O,  fill = id, shape = id))
p45 <- ggplot(chem, aes(x = Na2O, y = P2O5,  fill = id, shape = id))
p46 <- ggplot(chem, aes(x = K2O, y = P2O5,  fill = id, shape = id))

4.1 Custom tephra theme

To ensure consistency and to provide a shortcut of coding it everytime we set up a custom tephra plotting theme. This means all of our plots will have the same sized fonts etc…

Tephra.theme <-theme_bw()+theme(axis.title.x=element_text(size=16),
                                axis.title.y=element_text(size=16), 
                                axis.text.x=element_text(size=16),
                                axis.text.y=element_text(size=16))+
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

4.2 Custom colour palette

The hex numbers below specify a custom colour palette for plotting. These can be altered but were slected in order to help those with colour vision deficiencies, and this should be considered when making changes. At present 35 layers can be plotted using this pallete. If you have more than 35 groups of data you will either have to specify additional colours or consider splitting your dataset. The latter might be more appropriate for improving the final clarity in the plots.

chempalette <- c("#209A24","#A769EE","#B4C428","#e30303","#000000","#EF139D","#3A700F","#5288FF","#E64701","#7F9AF5","#F36D20","#0B3075","#C18522","#9634A6","#1D7D46","#80ddb8","#7CBC92","#B10E89","#384B27","#F896FD","#C6A766","#20628E","#9C0323","#C3AFE0","#572D24","#FF79A2","#3E6D6F","#E18C63","#6F2B40","#DF9FB4", "#291700", "#7e005b", "#00597a" , "#7b003a",  "#470600")

The palette is known as “chempalette”.

5.0 Biplots

We first generate the simple bi-plots of each element. This initial explorative approach facilitates data checking for unusual analyses and permits the identification of elements with significant variability.

TAS plot Le Maitre et al. (2002)

SiO2 VS TiO2

SiO2 VS Al2O3

SiO2 VS FeO

SiO2 VS MnO

SiO2 VS MgO

SiO2 VS CaO

SiO2 VS Na2O

SiO2 VS K2O

SiO2 VS P2O5

TiO2 VS Al2O3

TiO2 vs FeO

TiO2 vs MnO

TiO2 VS MgO

TiO2 VS CaO

TiO2 VS Na2O

TiO2 VS K2O

TiO2 VS P2O5

Al2O3 VS FeO

Al2O3 VS MnO

Al2O3 VS MgO

p21 <- p21 + geom_point(aes(fill = id), size=4, alpha=0.8)+xlab("Al2O3 (wt%)")+ylab("MgO (wt%)")+scale_shape_manual(values=c(21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25))+
  scale_fill_manual(values = chempalette) +
  Tephra.theme
p21

Al2O3 VS CaO

Al2O3 VS Na2O

Al2O3 VS K2O

Al2O3 VS P2O5

FeO VS MnO

FeO vs MgO

FeO vs CaO

FeO VS Na2O

FeO VS K2O

FeO VS P2O5

MnO VS MgO

MnO VS CaO

MnO VS Na2O

MnO VS K2O

MnO VS P2O5

MgO VS CaO

MgO VS Na2O

MgO VS K2O

MgO VS P2O5

CaO VS Na2O

CaO VS K2O

CaO VS P2O5

Na2O VS K2O

Na2O VS P2O5

K2O VS P2O5

5.1 Saving plots

We propose exporting your individual plots in .svg format. This is a vector file format that preserves your data without size and scaling issues. It is also readable and editable by software like Adobe Illustrator, Affinity Designer, or free alternatives like Inkscape. Using the file structure presented here the all plots can be named appropriately and will be saved in the Output -> Plots folder. The width and height and dpi outputs can be altered but override the preview functions in Rstudio and we recommend that you view the svg file directly in you preferred image editor to understand how it might look.

ggsave( filename = file.path("Script", "Output","Plots","fig1.svg"), plot = p1, width = 12, height = 6, dpi = 300)

6.0 Multivariate analysis

Multivariate analysis via ordination is often a helpful way of exploring your data. It permits the identification of useful elements for discrimination and facilitates dimension reduction. Unfortunately, percentage data can have problematic covariation and should be ‘freed’ from their compositional sum where possible. This can be achieved via undertaking a centred log-ratio transform. This is performed using the code given below

chem_multi <- chem[ -c(1,6,11,12,13,14,15)] ## check columns to drop. Currently dropping id, MnO, P2O5, Total, and TAS calculations to perform the clr transform.
chem_multi[is.na(chem_multi)] <- 0.005 # replace NAs with half the detection limit
chem_multi_clr <- clr(chem_multi)
cmclr <- as.data.frame(chem_multi_clr)
head(cmclr) # This permits us the opportunity to view the centred log-ratio dataset.

##       SiO2       TiO2     Al2O3        FeO        MgO        CaO       Na2O
## 1 2.135807 -0.5414021 0.7479836 0.80739231 -0.4486203  0.2756820 -0.9060452
## 2 3.332329 -2.1329626 1.5972904 0.59781585 -3.7084990 -0.7817596  0.7612333
## 3 3.996737 -2.8362948 2.1142366 0.06787029 -4.2225892 -1.5145389  0.9335885
## 4 3.716439 -2.4914800 1.8514592 0.10377467 -3.5900923 -1.2292383  0.8572538
## 5 2.026855 -0.5248841 0.6735616 0.83143798 -0.3107045  0.3313724 -0.7343422
## 6 3.317036 -2.0963366 1.5981721 0.53113307 -3.7385644 -0.7263028  0.7425532
##          K2O
## 1 -2.0707973
## 2  0.3345523
## 3  1.4609906
## 4  0.7818840
## 5 -2.2932962
## 6  0.3723095

We can summarise the importance of the resulting principal components by summarizing how much of the total variance is captured by each principal component.

6.1 PCA summary

Information	PC1	PC2	PC3	PC4	PC5	PC6	PC7	PC8
Standard deviation	2.452866	1.286943	0.5993116	0.1933593	0.1054135	0.07061895	0.06210361	4.673909e-16
Proportion of Variance	0.743770	0.204740	0.0444000	0.0046200	0.0013700	0.00062000	0.00048000	0.000000e+00
Cumulative Proportion	0.743770	0.948510	0.9929100	0.9975300	0.9989100	0.99952000	1.00000000	1.000000e+00

If this has captured a significant amount of variation it would be worth plotting the output in order to identify the important variables in determining the changes observed.

6.3 PCA scores scree plot

This is a visual representation of the summary table showing the percentage of the total variance explained by each dimension.

6.4 Variables PCA plot

6.5 Loadings dimension 1

6.6 Loadings dimension 2

6.7 PCA - Biplot

6.8 PCA - grouped

## Scale for 'shape' is already present. Adding another scale for 'shape', which
## will replace the existing scale.

7.0 Interactive figures

The function ggplotly within the plotly package is extremely useful for interactively examining data. ggplotly plots allow you to zoom into areas of interest and turn layers on and off as well as saving low resolution .png files of zoomed sections of your plots. We refer to these as iPlots (interactive plots) in order to distinguish them from static plots presented above. These are excellent for generating quick plots to support your inferences about your data. The plots can be exported as .png files.

iTAS plot

iFeO-CaO

iK-lines

iTiO2 vs Al2O3

iK2O vs CaO

8.0 Saving all figures to a .pdf

l <- list(p1, p2, p3, p4, p5, p6, p7, p8, p9, p47, p10, p11, p12, p13, p14, p15, p16, p17, p18, p19, p20, p21, p22, p23, p24, p25, p26, p27, p28, p29, p30, p31, p32, p33, p34, p35, p36, p37, p38, p39, p40, p41, p42, p43, p44, p45, p46)
pdf("allplotsout.pdf", width = 12, height = 6)
invisible(lapply(l, print))
dev.off()

## png 
##   2

9.0 Harker plot

It is frequently useful to plot all of the collected elemental data against SiO[2]. This permits an assessment of groupings within the data and which elements might be useful discriminators. Additionally, wt% K2O vs SiO2 can further classify the tephra to low, mid and high Potassium values. The dividing lines included here are from xxxxxxxx.

g1 <- plot_grid(p1 + theme(legend.position = "none"), 
                p2 + theme(legend.position = "none"),
                p3 + theme(legend.position="none"),
                p4 + theme(legend.position="none"),
                p5 + theme(legend.position="none"),
                p6 + theme(legend.position="none"),
                p7 + theme(legend.position="none"),
                p8 + theme(legend.position="none"), 
                p47 + theme(legend.position="none"), 
                p10 + theme(legend.position="none"),
                labels = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J"), align = "b")
legend <- get_legend(p1)
legplot <- plot_grid(legend)
px <- plot_grid( g1, rel_widths = c(4, 1))
l2 <- list(px,legplot)
pdf("Script/Output/Plots/Harkerplot.pdf", width = 30, height = 20)
invisible(lapply(l2, print))
dev.off()

## png 
##   2

ggsave( filename = file.path("Script","Output","Plots","Harker.svg"), plot = px, width = 30, height = 20, dpi = 300)

print(g1)

10.0 Other useful plots

10.1 Creating an inset within a figure

p100 <- ggplot(chem, aes(x = TAS_x, y = TAS_y, fill = id, shape = id))
p100 <- p100 + geom_point(aes(fill = id), size=4, alpha=0.8)+xlab("SiO2 (wt%)")+ylab("Na2O+K2O (wt%)")+
  scale_shape_manual(values=c(21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25))+
  scale_fill_manual(values = chempalette) +geom_segment(x=77,y=0,xend=69,yend=8)+ 
  geom_segment(x=69, y=8,xend=69,yend=13)+geom_segment(x=41, y=7,xend=52.5,yend=14)+
  geom_segment(x=52.5,y=14,xend=57.6,yend=11.7)+geom_segment(x=45,y=0,xend=45,yend=5)+
  geom_segment(x=45,y=5,xend=63,yend=14.56)+geom_segment(x=52,y=0,xend=52,yend=5)+
  geom_segment(x=52,y=5,xend=63,yend=7)+geom_segment(x=57,y=0,xend=57,yend=5.9)+
  geom_segment(x=63,y=0,xend=63,yend=7)+geom_segment(x=45,y=5,xend=52,yend=5)+
  geom_segment(x=52,y=5,xend=49.4,yend=7.3)+geom_segment(x=57,y=5.9,xend=53,yend=9.3)+
  geom_segment(x=63,y=7,xend=57.6,yend=11.7)+geom_segment(x=49.4,y=7.3,xend=45,yend=9.4)+
  geom_segment(x=53,y=9.3,xend=48.4,yend=11.5)+geom_segment(x=41,y=0,xend=41,yend=7)+
  geom_segment(x=52.5,y=14,xend=49,yend=15.5)+geom_segment(x=41,y=3,xend=45,yend=3)+
  geom_segment(x=63,y=7,xend=69,yend=8)+
  scale_x_continuous(limits=c(74, 79), expand = c(0, 0)) + ## set x axis of insert
  scale_y_continuous(limits=c(6, 10.5), expand = c(0, 0)) + ## set y axis of insert
  theme_bw()+theme(axis.title.x=element_text(size=9), axis.title.y=element_text(size=9), 
                   axis.text.x=element_text(size=9), axis.text.y=element_text(size=9))+ ## change the text size as neccessary
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none")
fig_in <- p1 + annotation_custom(ggplotGrob(p100), xmin = 40, xmax = 56, ymin = 7, ymax = 16) + ## place the insert
   annotate("rect", xmin = 74, xmax = 79, ymin = 6, ymax = 10.5,
  alpha = .2) ## use this to highlight the area of the insert plot

ggsave( filename = file.path("Script", "Output","Plots","insert_plot.svg"), plot = fig_in, width = 10, height = 6, dpi = 300)

print(fig_in)

10.2 Creating a HDR plot to simplify your data.

Simplifying the data we create to make more visually compelling figures can be achieved using highest density regions. The following code permits this as an addition to the previously created plots and as new plots. It can be added easily to any of the previously created plots by changing the p code. Visibility of the density regions can be applied to just reference data or entire datasets. Opacity is varied using the alpha command and the probability range can be altered in the ‘probs’ argument. You may need to vary the x and y axes in order to fit the full range in to the plot. These types of plots are best used for giving better eruptive chemical envelopes (fields) than hand drawing around existing data.

p28den <- p28 +
geom_hdr(stat= "hdr", method = "kde", na.rm = TRUE, inherit.aes = TRUE, alpha  = 0.3, xlim = c(0,100), ylim = c(0,30), smooth = TRUE, n = 1000, probs = c(0.75))
print(p28den)

10.3 Faceted plots to deal with lots of data

The kde plots are interesting but messy and still have lots of overplotting. This is not ideal for data visualisation. While it may be helpful to vary the colours or opacity, there is only so much that this can acheive. If you have particularly messy data like the ODP-980 dataset I’d suggest creating plot ‘facets’ using the following code.

p28den_facet <- p28 +
geom_hdr(method = "kde", alpha  = 0.5, xlim = c(0,100), ylim = c(0,30), show.legend = FALSE, smooth = TRUE, n = 1000, probs = c(0.95)) +
  facet_wrap(~ id)


print(p28den_facet)

The facet removes some of the complexity and permits the kde to be visualised more clearly, but arguably you lose some of the key information about total ranges of chemical data. In the second faceted plot (below) we remove the kernel density information and use the total dataset as a background to the facets in order to understand the total range and our own individual samples. This requires a little bit of data manipulation in order to create two datasets, one for the total data spread and another for the facets.

chem2 <- dplyr::select(chem, -Total) # creates the second dataset fpr plotting
chem2 <- rename(chem2, id2=id) # removes confusion over column names
f19 <- ggplot(chem, aes(x = FeO, y = CaO,  fill = id, shape = id)) # creates a new plot of FeO and CaO.

f19 <- f19 + 
    geom_point(data = chem2, inherit.aes = FALSE, aes(x=FeO, y = CaO, shape = id2), colour = "grey70") + ## adds all the data to the plot using a grey colour scheme
  geom_point(aes(fill = id), size=4, alpha=0.8)+xlab("FeO (wt%)")+ylab("CaO (wt%)")+scale_shape_manual(values=c(21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25))+ # adds all the data again but retain the groups
  scale_fill_manual(values = chempalette) + ## adds standard colours
  Tephra.theme + # adds standard theme
  coord_cartesian(xlim = c(0,20), ylim = c(0,15)) + # scale the plots to display the data effectively
  facet_wrap(~ id) # facets the diagram using id in the main dataset.

print(f19)

ggsave( filename = file.path("Script", "Output","Plots","faceted plot.pdf"), plot = f19, width = 30, height = 20, units = "cm", dpi = 300)

11.0 Comparing unknown samples

All of the above code can be used on unknown samples and for your own data exploration. It can also be used to explore a reference dataset in order to understand what chemical plots can be used to best differentiate tephra from different volcanoes and even different eruptions. However, some times it is better to keep the exploration and analysis phases separate. This section looks at importing a specific unknown tephra and making direct comparisons to a reference dataset. It uses the ODP-980 data presented above as a record of volcanism in the North Atlantic during MIS-11c and then compares this with our unknown tephra from a site in England called Marks Tey. First we load and summarise the data as before but then we look at filtering and presenting the data more clearly.

This code loads our unknown data from a normalised dataset called MT_norm.csv. We then present all of the data in a table.

unknowns <- here::here("Data","MT_norm.csv") %>% read.csv()
gt(unknowns)

id2	SiO2	TiO2	Al2O3	FeO	MnO	MgO	CaO	Na2O	K2O	P2O5	Total
MT-15.47	77.46191	0.12271193	11.66786	2.147459	0.07158196	0.00000000	0.7874016	4.591472	3.149606	0.00000000	100
MT-15.47	77.52786	0.12158055	11.24620	2.249240	0.05065856	0.02026343	0.7294833	4.528875	3.515704	0.01013171	100
MT-15.47	77.59528	0.13955985	11.70156	2.190016	0.07514761	0.03220612	0.7407407	4.551798	2.973698	0.01073537	100
MT-15.47	77.23816	0.11152793	11.84224	2.220420	0.06083342	0.00000000	0.6894454	4.430701	3.406671	0.01013890	100
MT-15.47	77.40131	0.12433945	11.94695	2.124132	0.06216972	0.07253134	0.6320589	4.403689	3.212102	0.01036162	100
MT-15.47	77.05858	0.12507817	11.79904	1.928289	0.09380863	0.01042318	0.7504690	4.732124	3.481343	0.01042318	100
MT-15.47	76.96308	0.12480499	11.93968	2.132085	0.08320333	0.00000000	0.7696308	4.492980	3.473739	0.02080083	100
MT-15.47	77.23074	0.12597103	11.11694	2.466933	0.05248793	0.05248793	0.8188117	4.797397	3.327735	0.00000000	100
MT-15.47	77.54255	0.14893617	11.34043	2.297872	0.08510638	0.05319149	0.6702128	4.808511	3.063830	0.00000000	100
MT-15.47	75.02112	0.08449514	13.45585	1.373046	0.06337136	-0.02112378	0.7498944	4.002957	5.259823	0.01056189	100
MT-15.47	77.76506	0.12485694	11.50765	1.966497	0.07283321	0.04161898	0.6867131	4.557278	3.267090	0.01040474	100
MT-15.47	77.22352	0.12736149	11.83401	2.133305	0.06368075	0.02122692	0.7004882	4.733602	3.173424	0.00000000	100
MT-15.47	78.07481	0.12433945	11.36670	2.041239	0.08289296	0.00000000	0.6527821	4.403689	3.253549	0.01036162	100
MT-15.47	77.86376	0.13236941	11.55687	2.026270	0.05091131	0.01018226	0.6720293	4.470013	3.197230	0.01018226	100
MT-15.47	77.59115	0.13314216	11.68578	1.915199	0.06145022	0.00000000	0.6964359	4.475625	3.430971	0.01024170	100
MT-15.47	77.21947	0.13723213	11.89697	2.216827	0.07389423	-0.02111264	0.8022802	4.496991	3.156339	0.01055632	100
MT-15.47	77.39612	0.12528712	11.67258	2.213406	0.09396534	0.01044059	0.6264356	4.740029	3.111297	0.01044059	100
MT-15.47	76.98995	0.12568077	11.86636	2.199413	0.08378718	0.01047340	0.6912442	4.639715	3.382907	0.02094680	100
MT-15.47	77.37722	0.13584117	11.58830	2.361546	0.06269592	0.03134796	0.8045977	4.420063	3.218391	0.01044932	100

Like before this is useful to use to understand that our data has loaded correctly, but a summary table might be more useful. This is performed using the code from section 2.2.

unknown_sum <- unknowns %>% group_by(id2) %>% dplyr::summarise(                                            
                                                       count=n(),
                                                       SiO2.= mean(SiO2),
                                                       sd1 = sd(SiO2),
                                                       TiO2. = mean(TiO2),
                                                       sd2 = sd(TiO2),
                                                       Al2O3. = mean(Al2O3),
                                                       sd3 = sd(Al2O3),
                                                       FeOt = mean(FeO),
                                                       sd4 = sd(FeO),
                                                       MnO.= mean(MnO),
                                                       sd5 = sd(MnO),
                                                       MgO.= mean(MgO),
                                                       sd6 = sd(MgO),
                                                       CaO.= mean(CaO),
                                                       sd7 = sd(CaO),
                                                       Na2O. = mean(Na2O),
                                                       sd8 = sd(Na2O),
                                                       K2O. = mean(K2O),
                                                       sd9 = sd(K2O),
                                                       P2O5. = mean(P2O5),
                                                       sd10 = sd(P2O5),
                                                       Total.=mean(Total),
                                                       sd11 = sd(Total))
unknown_sum <- unknown_sum %>% mutate(across(where(is.numeric), ~ round(., digits = 2)))
unknown_sum %>% rename(layer = id2, n = count , SiO2 = SiO2., TiO2 = TiO2.)

## # A tibble: 1 x 24
##   layer     n  SiO2   sd1  TiO2   sd2 Al2O3.   sd3  FeOt   sd4  MnO.   sd5  MgO.
##   <chr> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 MT-1~    19  77.3  0.62  0.13  0.01   11.7  0.48  2.12  0.23  0.07  0.01  0.02
## # ... with 11 more variables: sd6 <dbl>, CaO. <dbl>, sd7 <dbl>, Na2O. <dbl>,
## #   sd8 <dbl>, K2O. <dbl>, sd9 <dbl>, P2O5. <dbl>, sd10 <dbl>, Total. <dbl>,
## #   sd11 <dbl>

gt(unknown_sum)

id2	count	SiO2.	sd1	TiO2.	sd2	Al2O3.	sd3	FeOt	sd4	MnO.	sd5	MgO.	sd6	CaO.	sd7	Na2O.	sd8	K2O.	sd9	P2O5.	sd10	Total.	sd11
MT-15.47	19	77.29	0.62	0.13	0.01	11.74	0.48	2.12	0.23	0.07	0.01	0.02	0.03	0.72	0.06	4.54	0.19	3.37	0.48	0.01	0.01	100	0

We evidently have a relatively homogenous highly silicic tephra with average silica values of 77.29 %wt and FeOt values below 3%wt. We could undertake further plots to check this out, but it might also be an examination of the tabular data is sufficient for this purpose. This helps us to filter our possible correlatives. We can also restrict our possible correlatives using other stratigraphic and chronological data. This might be plausible time windows that the tephra could occur or other factors. The ODP-980 data have some additional information that might aid some discrimination. The publication lists sections of Ice rafted debris (IRD) which reflect mixtures of ash layers and therefore not primary tephra deposits. We will exclude these here in order to only deal with potential primary ash layers.

chem_red <- filter(chem, IRD == "NO") #We might want to filter our data by age or other stratigraphic information. here we are excluding layers identified as ice rafted debris. This can be achieved using additional columns in the input data sheet. All of the plotting code will ignore these unless told otherwise.

chem_red <- filter(chem_red, SiO2 >68 & FeO <4) # Our unknown sample is a rhyolite so we can now employ our filtering by value, we can also exclude some of the high iron examples in the ODP-980 dataset.

c28 <- ggplot(chem_red, aes(x = SiO2, y = CaO, fill = id, shape = id)) # create a new plot using the reduced chemical comparison data
c28 <- c28 + 
  geom_hdr(method = "kde", alpha  = 0.4, xlim = c(0,100), ylim = c(0,30), show.legend = TRUE, smooth = TRUE, n = 500, probs = c(0.95), color = "black") + # Add the ODP-980 data as a kde field for comparison this uses 95% probability ranges for comparison.
   scale_fill_viridis(discrete = TRUE) + # We will use the Viridis scale for now as this gives relatively good differentiation.
    geom_point(data = unknowns, inherit.aes = FALSE, aes(x=SiO2, y = CaO, shape = id2, fill = id2), size=4, alpha=0.7, show.legend = FALSE) + # we take our unknown data and plot it as points
  xlab("SiO2 (wt%)") + # Label our axes
  ylab("CaO (wt%)") + 
  scale_shape_manual(values=c(1:5,7,11, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25)) +
  Tephra.theme # add the tephra theme.
print(c28)

We can see that only samples 54.79-54.83 and 51.09-51.13 have any overlap with our unknown at 95% probability. It is worth noting here that this is in part down to our filtering of the data in the previous step. We need to be sure that we are not excluding important data when undertaking this step otehrwise our resulting kdes can be strongly altered. Therefore it is also sensible to plot the raw point data for comparison also. The interactive plot below uses the same elements but permits exploration of the data. From this we can see that 54.79-54.83 is the best match for our sample on these elements.

u28 <- ggplot(chem_red, aes(x = SiO2, y = CaO,  fill = id, shape = id))
u28 <- u28 + geom_point(aes(fill = id), size=4, alpha=0.8)+xlab("SiO2 (wt%)")+ylab("CaO (wt%)")+scale_shape_manual(values=c(21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25, 21:25))+
      geom_point(data = unknowns, inherit.aes = FALSE, aes(x=SiO2, y = CaO, shape = id2, fill = id2), size=4, alpha=0.8) +
  scale_fill_manual(values = chempalette) +
  Tephra.theme

ggplotly(u28, width = 1000, height = 800)

References

R. W. Le Maitre (editor), A. Streckeisen, B. Zanettin, M. J. Le Bas, B. Bonin, P. Bateman, G. Bellieni, A. Dudek, S. Efremova, J. Keller, J. Lamere, P. A. Sabine, R. Schmid, H. Sorensen, and A. R. Woolley, Igneous Rocks: A Classification and Glossary of Terms, Recommendations of the International Union of Geological Sciences, Subcommission of the Systematics of Igneous Rocks. Cambridge University Press, 2002.

Papers using versions of this code.

Walsh, A.A., Blockley, S.P., Milner, A.M. and Martin-Puertas, C., 2023. Updated age constraints on key tephra markers for NW Europe based on a high-precision varve lake chronology. Quaternary Science Reviews, 300, p.107897. read it here

Carter‐Champion, A., Abrook, A.M., Pike, J.H., Matthews, I.P., Palmer, A.P. and Lowe, J.J., 2022. Ice‐sheet deglaciation and Loch Lomond Readvance in the eastern Cairngorms: implications of a Lateglacial sediment record from Glen Builg. Journal of Quaternary Science, 37(8), pp.1332-1347. read it here

Walsh, A.A., Blockley, S.P., Milner, A.M., Matthews, I.P. and Martin-Puertas, C., 2021. Complexities in European Holocene cryptotephra dispersal revealed in the annually laminated lake record of Diss Mere, East Anglia. Quaternary Geochronology, 66, p.101213. read it here

Palmer, A.P., Matthews, I.P., Lowe, J.J., MacLeod, A. and Grant, R., 2020. A revised chronology for the growth and demise of Loch Lomond Readvance (‘Younger Dryas’) ice lobes in the Lochaber area, Scotland. Quaternary Science Reviews, 248, p.106548.read it here

Lowe, J., Matthews, I., Mayfield, R., Lincoln, P., Palmer, A., Staff, R. and Timms, R., 2019. On the timing of retreat of the Loch Lomond (‘Younger Dryas’) Readvance icefield in the SW Scottish Highlands and its wider significance. Quaternary Science Reviews, 219, pp.171-186.read it here