# transform_piper_data for data transformation:
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))
}
# ggplot_piper for piper diagram background:
ggplot_piper <- function() {
library(ggplot2)
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 <- ggplot() +
# Left hand ternary plot:
geom_polygon(data=polygons2, aes(x=x,y=y,fill=Major.Water.Type,
group=Major.Water.Type)) +
geom_segment(aes(x=0,y=0, xend=100, yend=0)) +
geom_segment(aes(x=0,y=0, xend=50, yend=86.603)) +
geom_segment(aes(x=50,y=86.603, xend=100, yend=0)) +
# Right hand ternary plot:
geom_polygon(data=polygons3, aes(x=x,y=y,fill=Major.Water.Type,
group=Major.Water.Type)) +
geom_segment(aes(x=120,y=0, xend=220, yend=0)) +
geom_segment(aes(x=120,y=0, xend=170, yend=86.603)) +
geom_segment(aes(x=170,y=86.603, xend=220, yend=0)) +
# Upper diamond:
geom_polygon(data=polygons, aes(x=x,y=y,fill=Major.Water.Type,
group=Major.Water.Type)) +
geom_segment(aes(x=110,y=190.5266, xend=60, yend=103.9236)) +
geom_segment(aes(x=110,y=190.5266, xend=160, yend=103.9236)) +
geom_segment(aes(x=110,y=17.3206, xend=160, yend=103.9236)) +
geom_segment(aes(x=110,y=17.3206, xend=60, yend=103.9236)) +
# Add grid lines to the plots:
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p1,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p2,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid1p3,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p1,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p2,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid2p3,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid3p1,
linetype = "dashed", size = 0.5, colour = "grey50") +
geom_segment(aes(x=x1, y=y1, yend=y2, xend=x2), data=grid3p2,
linetype = "dashed", size = 0.5, colour = "grey50") +
# Labels and grid values
geom_text(aes(c(20,40,60,80),c(-5,-5,-5,-5), label=c(80, 60, 40, 20)), size=2) +
geom_text(aes(c(35,25,15,5),grid1p2$y2, label=c(80, 60, 40, 20)), size=2) +
geom_text(aes(c(95,85,75,65),grid1p3$y2, label=c(80, 60, 40, 20)), size=2) +
coord_equal(ratio=1)+
geom_text(aes(17,50, label="Mg^2"), angle=60, size=4, parse=TRUE) +
geom_text(aes(82.5,50, label="Na + K"), angle=-60, size=4) +
geom_text(aes(50,-10, label="Ca^2"), size=4, parse=TRUE) +
geom_text(aes(170,-10, label="Cl^-phantom()"), size=4, parse=TRUE) +
geom_text(aes(205,50, label="SO^4"), angle=-60, size=4, parse=TRUE) +
geom_text(aes(137.5,50, label="Alkalinity~as~HCO^3"), angle=60, size=4, parse=TRUE) +
geom_text(aes(72.5,150, label="SO^4~+~Cl^-phantom()"), angle=60, size=4, parse=TRUE) +
geom_text(aes(147.5,150, label="Ca^2~+~Mg^2"), angle=-60, size=4, parse=TRUE) +
geom_text(aes(c(155,145,135,125),grid2p2$y2, label=c(20, 40, 60, 80)), size=2) +
geom_text(aes(c(215,205,195,185),grid2p3$y2, label=c(20, 40, 60, 80)), size=2) +
geom_text(aes(c(140,160,180,200),c(-5,-5,-5,-5), label=c(20, 40, 60, 80)), size=2) +
geom_text(aes(grid3p1$x1-5,grid3p1$y1, label=c(80, 60, 40, 20)), size=2) +
geom_text(aes(grid3p1$x2+5,grid3p1$y2, label=c(20, 40, 60, 80)), size=2) +
geom_text(aes(grid3p2$x1-5,grid3p2$y1, label=c(20, 40, 60, 80)), size=2) +
geom_text(aes(grid3p2$x2+5,grid3p2$y2, label=c(80, 60, 40, 20)), size=2) +
theme_bw() +
theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.border = element_blank(), axis.ticks = element_blank(),
axis.text.x = element_blank(), axis.text.y = element_blank(),
axis.title.x = element_blank(), axis.title.y = element_blank())
return(p)
}
# Transform the data into piper based coordinate:
piper_data <- transform_piper_data(Ca = ions$Ca, Mg = ions$Mg, Cl = ions$Cl,
SO4= ions$SO4, name = ions$Sample.ID)
# Upper Diamond
ids <- factor(c("Sodium Bicarbonate", "Sodium Chloride",
"Calcium Bicarbonate", "Calcium Sulfate"))
values <- data.frame(Major.Water.Type = ids, value = c(1,2,3,4))
positions <- data.frame(Major.Water.Type=rep(ids, each = 4),
x=c(110,85,110,135,
135,110,135,160,
85,60,85,110,
110,85,110,135),
y=c(17.3206,60.6221, 103.9236,60.6221,
60.6221, 103.9236, 147.2251, 103.9236,
60.6221,103.9236,147.2251,103.9236,
103.9236,147.2251,190.5266,147.2251))
polygons <- merge(values, positions)
# Left Ternary Plot
ids2 <- factor(c("5", "6", "7", "8"))
values2 <- data.frame(Major.Water.Type = ids, value = c(5,6,7,8))
positions2 <- data.frame(Major.Water.Type=rep(ids2, each = 3),
x=c(50,0,25,
50,25,75,
100,50,75,
75,25,50),
y=c(0,0,43.3015,
0,43.3015,43.3015,
0,0,43.3015,
43.3015,43.3015,86.603))
polygons2 <- merge(values2, positions2)
# Right Ternary Plot
ids3 <- factor(c("9", "10", "11", "12"))
values3 <- data.frame(Major.Water.Type = ids, value = c(9,10,11,12))
positions3 <- data.frame(Major.Water.Type=rep(ids2, each = 3),
x=c(170,120,145,
170,145,195,
220,170,195,
195,145,170),
y=c(0,0,43.3015,
0,43.3015,43.3015,
0,0,43.3015,
43.3015,43.3015,86.603))
polygons3 <- merge(values3, positions3)
# Then, add a line for each 'plot' (upper, left, right) within the ggplot_piper
#function above
