Here I show my efforts to generate an interactive Piper digram. It was inspired by code from http://srmulcahy.github.io/2012/12/04/ternary-plots-r.html https://github.com/markolipka/ggplot_Piper and functions in the smwrGraphs package.
The function toPercent takes as input a data frame with the cations and anion concentrations (miliequivalent) and converts them to percentages.
toPercent <- function (d) {
totalCations <- d$Ca + d$Mg + d$Na + d$K
d$Ca <- 100 * (d$Ca/totalCations)
d$Mg <- 100 * (d$Mg/totalCations)
d$Na <- 100 * (d$Na/totalCations)
d$K <- 100 * (d$K/totalCations)
totalAnions <- d$Cl + d$SO4 + d$CO3 + d$HCO3
d$Cl <- 100 * (d$Cl/totalAnions)
d$SO4 <- 100 * (d$SO4/totalAnions)
d$CO3 <- 100 * (d$CO3/totalAnions)
d$HCO3 <- 100 * (d$HCO3/totalAnions)
return(d)
}Next the function transform_piper_data takes percentage data and converts them to x, y coordinates for plotting.
transform_piper_data <- function(Mg, Ca, Cl,SO4, name=NULL){
if(is.null(name)){
name = rep(1:length(Mg),3)
} else {
name = rep(name,3)
}
y1 <- Mg * 0.86603
x1 <- 100*(1-(Ca/100) - (Mg/200))
y2 <- SO4 * 0.86603
x2 <-120+(100*Cl/100 + 0.5 * 100*SO4/100)
new_point <- function(x1, x2, y1, y2, grad=1.73206){
b1 <- y1-(grad*x1)
b2 <- y2-(-grad*x2)
M <- matrix(c(grad, -grad, -1,-1), ncol=2)
intercepts <- as.matrix(c(b1,b2))
t_mat <- -solve(M) %*% intercepts
data.frame(x=t_mat[1,1], y=t_mat[2,1])
}
np_list <- lapply(1:length(x1), function(i) new_point(x1[i], x2[i], y1[i], y2[i]))
npoints <- do.call("rbind",np_list)
data.frame(observation=name,x=c(x1, x2, npoints$x), y=c(y=y1, y2, npoints$y))
}Now we are ready to create a ggplot or optionally an interactive plotly chart via the follwoing function. The required data frame piper.data has columns labeled x, y (the transformed coordinates for the plot) and for my own purposes columns called date, watertype (for grouping plotted points), and observation (a sample id).
library(dplyr)
library(smwrGraphs)
library(ggplot2)
library(plotly)
ggplot_piper <- function(piper.data,output = c("ggplot","plotly")) {
grid1p1 <<- data.frame(x1 = c(20,40,60,80), x2= c(10,20,30,40),y1 = c(0,0,0,0), y2 = c(17.3206,34.6412,51.9618, 69.2824))
grid1p2 <<- data.frame(x1 = c(20,40,60,80), x2= c(60,70,80,90),y1 = c(0,0,0,0), y2 = c(69.2824, 51.9618,34.6412,17.3206))
grid1p3 <<- data.frame(x1 = c(10,20,30,40), x2= c(90,80,70,60),y1 = c(17.3206,34.6412,51.9618, 69.2824), y2 = c(17.3206,34.6412,51.9618, 69.2824))
grid2p1 <<- grid1p1
grid2p1$x1 <- grid2p1$x1+120
grid2p1$x2 <- grid2p1$x2+120
grid2p2 <<- grid1p2
grid2p2$x1 <- grid2p2$x1+120
grid2p2$x2 <- grid2p2$x2+120
grid2p3 <<- grid1p3
grid2p3$x1 <- grid2p3$x1+120
grid2p3$x2 <- grid2p3$x2+120
grid3p1 <<- data.frame(x1=c(100,90, 80, 70),y1=c(34.6412, 51.9618, 69.2824, 86.603), x2=c(150, 140, 130, 120), y2=c(121.2442,138.5648,155.8854,173.2060))
grid3p2 <<- data.frame(x1=c(70, 80, 90, 100),y1=c(121.2442,138.5648,155.8854,173.2060), x2=c(120, 130, 140, 150), y2=c(34.6412, 51.9618, 69.2824, 86.603))
p <- ggplot2::ggplot() +
## left hand ternary plot
ggplot2::geom_segment(ggplot2::aes(x=0,y=0, xend=100, yend=0)) +
ggplot2::geom_segment(ggplot2::aes(x=0,y=0, xend=50, yend=86.603)) +
ggplot2::geom_segment(ggplot2::aes(x=50,y=86.603, xend=100, yend=0)) +
## right hand ternary plot
ggplot2::geom_segment(ggplot2::aes(x=120,y=0, xend=220, yend=0)) +
ggplot2::geom_segment(ggplot2::aes(x=120,y=0, xend=170, yend=86.603)) +
ggplot2::geom_segment(ggplot2::aes(x=170,y=86.603, xend=220, yend=0)) +
## Upper diamond
ggplot2::geom_segment(ggplot2::aes(x=110,y=190.5266, xend=60, yend=103.9236)) +
ggplot2::geom_segment(ggplot2::aes(x=110,y=190.5266, xend=160, yend=103.9236)) +
ggplot2::geom_segment(ggplot2::aes(x=110,y=17.3206, xend=160, yend=103.9236)) +
ggplot2::geom_segment(ggplot2::aes(x=110,y=17.3206, xend=60, yend=103.9236)) +
## Add grid lines to the plots
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p1, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p2, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p3, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p1, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p2, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p3, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid3p1, linetype = "dashed", size = 0.25, colour = "grey50") +
ggplot2::geom_segment(ggplot2::aes(x=x1, y=y1, yend=y2, xend=x2), data=grid3p2, linetype = "dashed", size = 0.25, colour = "grey50") +
### Labels and grid values
ggplot2::geom_text(ggplot2::aes(c(20,40,60,80),c(-5,-5,-5,-5), label=c(80, 60, 40, 20)), size=3) +
ggplot2::geom_text(ggplot2::aes(c(35,25,15,5),grid1p2$y2, label=c(80, 60, 40, 20)), size=3) +
ggplot2::coord_equal(ratio=1) +
ggplot2::geom_text(ggplot2::aes(c(215,205,195,185),grid2p3$y2, label=c(20, 40, 60, 80)), size=3) +
ggplot2::geom_text(ggplot2::aes(c(140,160,180,200),c(-5,-5,-5,-5), label=c(20, 40, 60, 80)), size=3) +
ggplot2::geom_text(ggplot2::aes(grid3p1$x1-5,grid3p1$y1, label=c(80, 60, 40, 20)), size=3) +
ggplot2::geom_text(ggplot2::aes(grid3p1$x2+5,grid3p1$y2, label=c(20, 40, 60, 80)), size=3) +
ggplot2::geom_text(ggplot2::aes(grid3p2$x1-5,grid3p2$y1, label=c(20, 40, 60, 80)), size=3) +
ggplot2::geom_text(ggplot2::aes(grid3p2$x2+5,grid3p2$y2, label=c(80, 60, 40, 20)), size=3) +
ggplot2::theme_bw() +
ggplot2::theme(panel.grid.major = ggplot2::element_blank(), panel.grid.minor = ggplot2::element_blank(),
panel.border = ggplot2::element_blank(), axis.ticks = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank(), axis.text.y = ggplot2::element_blank(),
axis.title.x = ggplot2::element_blank(), axis.title.y = ggplot2::element_blank(),
legend.title = element_blank()) +
ggplot2::geom_point(ggplot2::aes(x,y, colour=factor(watertype),
shape = watertype,
text = paste(observation,
'</br></br>Date: ',
date)), data=piper.data)
if (output == "ggplot"){
p <- p +
ggplot2::geom_text(ggplot2::aes(17,50, label="Mg^'2+'"), angle=60, size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(77.5,50, label="Na^'+'~+~K^'+'"), angle=-60, size=4,parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(50,-10, label="Ca^'2+'"), size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(170,-10, label="Cl^'-'"), size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(205,50, label="SO[4]^'-'"), angle=-60, size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(142,50, label="Alkalinity~as~{HCO[3]^'-'"), angle=60, size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(72.5,150, label="SO[4]^'-'~+~Cl^'-'"), angle=60, size=4, parse=TRUE) +
ggplot2::geom_text(ggplot2::aes(147.5,150, label="Ca^'2+'~+~Mg^'2+'"), angle=-60, size=4, parse=TRUE)
}
if (output == "plotly"){
#this fixes an issue that plotly can't render geom_text() with the angle option set properly
p <- plotly::ggplotly(p,
tooltip = c("text")
)
p <- p %>% plotly::layout(
annotations=list(text=c("Mg<sup>2+</sup>",
"Na<sup>+</sup> + K<sup>+</sup>",
"Ca<sup>2+</sup>",
"Cl<sup>-</sup>",
"SO<sub>4</sub><sup>-</sup>",
"Alkalinity as HCO<sub>3</sub><sup>-</sup>",
"SO<sub>4</sub><sup>-2</sup> + Cl<sup>-</sup>",
"Ca<sup>2+</sup> + Mg<sup>2+</sup>"),
x = c(17,77.5,50,170,205,142.5,72.5,147.5),
y = c(50,50,-10,-10,50,50,150,150),
textangle = c(-60,60,0,0,60,-60,-60,60),
"showarrow"=F, font=list(size = 12, color = "black")
))
}
return(p)
}