Piper Diagram

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)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0     ✔ purrr   0.3.2
## ✔ tibble  2.1.3     ✔ dplyr   0.8.3
## ✔ tidyr   0.8.3     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
library(readxl)
library(broom)
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.4-4, (SVN revision 833)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/gdal
##  GDAL binary built with GEOS: FALSE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.6/Resources/library/rgdal/proj
##  Linking to sp version: 1.3-1
library(sp)
library(cowplot)
## 
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
##   default ggplot2 theme anymore. To recover the previous
##   behavior, execute:
##   theme_set(theme_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.6

Converting 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             142

Plotting

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 Polygons
well_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