Making a Piper Plot
Piper plots are used for visualizing groundwater chemistry and can tell a useful story to help hydrogeologists better understand spatial relationships in water sources or water paths.
To run this code you will need several packages and a piece of code linked here.
In addition, you will need groundwater data. Please format your csv datafiles as shown below.
source("ggplot_Piper.R")Required Libraries
library(hydrogeo)
library(tidyverse)library(readxl)
library(broom)
library(rgdal)library(sp)
library(cowplot)Loading Data
data <- read.csv("water_chemistry_subset.csv", header = TRUE)
data %>%
filter(!is.na(Well.ID)) %>%
mutate(Well.ID = as.factor(Well.ID),
Formation = as.factor(Formation),
Well.Depth.Feet = as.factor(Well.Depth.Feet)) -> data
head(data)## Formation Well.ID Sample.Date Well.Depth.Feet Ca K Mg
## 1 HIGH CASCADE VOLCANICS 48N/1W-28Fl 7-18-78 632 2 2 1
## 2 HIGH CASCADE VOLCANICS 4BN/1W-28J2 7-18-78 1322 20 8 15
## 3 HIGH CASCADE VOLCANICS 48N/1W-34B 7-18-78 515 26 13 21
## 4 HIGH CASCADE VOLCANICS 4BN/1W-34Gl 7-18-78 383 14 20 7
## 5 HIGH CASCADE VOLCANICS 47N/1W-2J 7-20-78 1131 5 10 2
## 6 HIGH CASCADE VOLCANICS 47N/1W-2 3Hl 7-20-78 1065 8 10 4
## Na Cl SO4 HCO3 CO3
## 1 42 5 5 117 5.4
## 2 26 6 5 203 67.2
## 3 42 7 5 291 90.6
## 4 75 12 5 295 38.4
## 5 29 5 5 111 12.6
## 6 40 7 5 149 21.6Converting milliequivalents to percents
data <- toPercent(data)Check that sum adds to 100%
cation.sums <- apply(data[, c("Ca", "Mg", "Na", "K")], 1, FUN = sum)
anion.sums <- apply(data[, c("Cl", "SO4", "CO3", "HCO3")], 1, FUN = sum)Transforming Data
piper_data <- transform_piper_data(Ca = data$Ca,
Mg = data$Mg,
Cl = data$Cl,
SO4 = data$SO4,
name = data$Well.ID)Merging with previous data
piper_data <- merge(piper_data,
data[,c("Formation", "Well.ID", "Well.Depth.Feet")],
by.x = "observation",
by.y = "Well.ID")
head(piper_data)## observation x y Formation
## 1 46N/1E-80 89.77568 60.233569 HIGH CASCADE VOLCANICS
## 2 46N/1E-80 122.92588 2.815442 HIGH CASCADE VOLCANICS
## 3 46N/1E-80 65.00000 17.320600 HIGH CASCADE VOLCANICS
## 4 46N/1W-18Ql 48.90110 25.695396 LAKE DEPOSITS OR BUTTE VALLEY BASALT
## 5 46N/1W-18Ql 82.24130 83.442625 LAKE DEPOSITS OR BUTTE VALLEY BASALT
## 6 46N/1W-18Ql 125.68182 8.201042 LAKE DEPOSITS OR BUTTE VALLEY BASALT
## Well.Depth.Feet
## 1 230
## 2 230
## 3 230
## 4 142
## 5 142
## 6 142Plotting
fig.plot <- ggplot_piper() +
geom_point(data = piper_data, mapping = aes(x,y, color=Formation), size = 1, stroke=1) +
scale_shape_manual(values = c(21,23)) + labs(color = "Well") +
theme(legend.position = c(0.1, 0.9),
legend.title = element_text(face = "bold", size = 6),
legend.text = element_text(size = 6))
fig.plot
Combining Points With Data
well_locations <- readOGR(dsn = ".", layer = "water_chemistry_subset", verbose = FALSE)Loading Shapefiles
watershed <- readOGR(dsn = ".", layer = "border", verbose = FALSE)Transform Watershed and points
watershed2 <- fortify(watershed) %>% mutate(id = as.numeric(id))## Regions defined for each Polygonswell_locations2 <- as_data_frame(well_locations)## Warning: `as_data_frame()` is deprecated, use `as_tibble()` (but mind the new semantics).
## This warning is displayed once per session.Plotting Map
fig.map <- ggplot() +
geom_polygon(data = watershed2, mapping = aes(long,lat, group = group), fill = "grey50") +
geom_point(data = well_locations2, mapping = aes(coords.x1, coords.x2, color=Type)) +
scale_shape_manual(values = c(21,23)) + coord_equal() +
theme_void() + theme(legend.position = "none")
plot(fig.map)
Plot it
## Plotting
fig.plot <- plot_grid(fig.map,fig.plot, labels = c("(A)", "(B)"), rel_widths = c(1,2))
fig.plot