knitr::opts_chunk$set(echo = FALSE)

# not all libraries are used in this script
library(tidyverse) #dplyr, ggplot2, tidyr, readr, purrr, tibble, stringr, forcats
library(reshape2)
library(lubridate)  #for specifying dates and times
library(pander)
library(xtable)
library(gridExtra)
library(ggthemes)
library(ggrepel)
library(gganimate)  #for data graphics animations
library(gifski)     #works with gganimate
library(png)
library(readxl)
library(rvest) # for scraping
library(httr)  # for web APIs
library(transformr)
library(data.table)
library(RColorBrewer)  #color palettes for graphs
library(dataRetrieval)  #usgs data retrieval tool

Historical data

Every month MyRWAS has sampled at 10 sites around the watershed for 20 years. Specific conductance data over time at each site is depicted in this figure. Trends at all sites seem to be positive, suggesting increasing chloride concentrations in freshwater lakes and streams over time.

# load exported water quality data

wq <- read.csv("/Users/Andy/Mystic River Dropbox/Projects/water quality/Chloride/2020/wq_export_20200115.csv")

# use lubridate package to rationalize Datetime field

wq$Datetime <- ymd_hms(wq$Datetime)

#filter for parameter specific conductance, only baseline data, and only freshwater sites, 
# and add a column for season

spcond_base <- wq %>%
  filter(CharacteristicID=="SPCOND") %>%
  filter(ProjectID=="BASE")%>%
  filter(WaterType=="Fresh")%>%
  mutate(season = ifelse(month(Datetime) %in% c(1,2,12), "winter", 
                                  ifelse(month(Datetime)%in% c(3,4,5), "spring",
                                         ifelse(month(Datetime)%in% c(7,8,9), "summer", "fall"))))
         

g <- ggplot(spcond_base, aes(Datetime, ResultValue, ymin=0, color=season), na.omit=TRUE)
g + geom_point() + geom_line() + facet_wrap(~LocationID, scales="free_y") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1)) +scale_color_brewer(palette="Spectral")

Mystic River conductivity Dec 2019-Jan 2020

In 2019 we deployed a conductivity sensor in Mystic River at 15 West Street, Medford. Flow and precipitation data drawn from USGS flow gage data from Aberjona River.

Note that conducitivity goes up in precipitation and snow-melt events, as salt applied to roads in winter is transported to the river.

#load conductivity sensor, create day and hour columns, join with data.table, and convert to long format

sensor <- read.csv("/Users/Andy/Mystic River Dropbox/Projects/water quality/Chloride/2020/COND1_20200107.csv")

sensor$Datetime <- mdy_hm(sensor$Datetime)

#load USGS data

library(dataRetrieval) # see documentation 
ABR <- "01102500" # site id
discharge <- "00060"  #parameter code
height <- "00065"
precip <- "00045"

# readNWISdv    siteNumber  NWIS daily data (parameterCd,startDate,endDate,)statCd  
# readNWISqw    siteNumber  NWIS water quality data
# readNWISuv    siteNumber  NWIS instantaneous value data
# readNWISrating    siteNumber  NWIS rating table for active streamgage
# type

#instantaneous value data for discharge 

instant_ABR <- readNWISuv(ABR,discharge,
                      "2019-12-13","2020-01-07")

#instant_ABR$dateTime <- as.Date(instant_ABR$dateTime)

g <- ggplot(instant_ABR, aes(dateTime, X_00060_00000))

g + geom_line()

# alternative way to load USGS data: upload from manually copied data request from website: 
# 
USGS <- read.csv("/Users/Andy/Mystic River Dropbox/Projects/water quality/Chloride/2020/ABR_flow.csv")

USGS$Datetime <- mdy_hm(USGS$Datetime)

#create data tables to use data.table functions
sensor_dt <- data.table(sensor)
USGS_dt <- data.table(USGS)

#join data tables on single nearest time match, to avoid multiple row duplication
# see: https://atrebas.github.io/post/2019-03-03-datatable-dplyr/
#using data.table package

joined_dt <- USGS_dt[sensor_dt, on = .(Datetime==Datetime), roll = "nearest"]

#change to long format

joined_dt_long <- joined_dt %>%
  gather(parameter, value, 5,7,9,11,12 )

#join in dplyr creates duplicate records
# joined <- left_join(spcond, USGS, by = c("date", "hour"))
# joined <- joined %>%
#   gather(parameter, value, 2,3,11,13)

#graph 

g <- ggplot(joined_dt_long, aes(Datetime, value, color=parameter)) + geom_line() + facet_wrap(~parameter, scales = "free_y") + scale_color_brewer(palette="Set2")
g
## Warning: Removed 2 rows containing missing values (geom_path).

##Data animation

## Warning: Removed 2 rows containing missing values (geom_path).