# Model Archive
# Method finalized March 28, 2017
# Minor modifications and additional comments made 
# August 10, 2017, August 26, 2017, and September 7, 2017
# The U.S. Geological Survey Peak-Flow File Data Verification Project, 2008-16
#
# Must be connected to the Internet to download packages and retrieve 
# latitude and longitude data
#
# Necessary R packages
# If already installed on system, install.packages() commands
# can be skipped. Listed separatly for those that may need to
# install only some of them.

install.packages("dplyr")
install.packages("plyr")
install.packages("dataRetrieval")
install.packages("maps")
install.packages("mapproj")

# attach the libraries
library(dplyr)
library(dataRetrieval)
library(maps)
library(mapproj)

# Data retrieval
# U.S. Geological Survey data release supporting this report is available at  
# http://dx.doi.org/10.5066/F7NG4NQ4. The data release contains three files and 
# their associated metadata files that provide a public summary of the changes 
# made to the U.S. Geological Survey peak flow file since formal Nationwide 
# checking began in 2008. The following code assumes one has downloaded the 
# November 19, 2008, and November 14, 2016, data files and placed them in the 
# same directory in which this script runs.  

pff2008 <- read.csv("Peak Flow File November-19-2008.csv",
                   colClasses = c("character", "numeric", 
                   rep("character", 2), "numeric", "character", 
                   "numeric"), stringsAsFactors = FALSE)

pff2016 <- read.csv("Peak Flow File November-14-2016.csv",
                   colClasses = c("character", "numeric", 
                   rep("character", 2), "numeric", "character", 
                   "numeric"), stringsAsFactors = FALSE)
pff2016 <- subset(pff2016, !is.na(peak_va))

allpkdat <- merge(pff2008[, c(1:3, 5:8)], pff2016[, c(1:3, 5:8)],
                 by = c("site_no", "WaterYr"))


# The following code produces counts presented in the report.

# peaks in 2008
prettyNum(dim(pff2008)[1], big.mark = ",")

# peaks in 2016
prettyNum(dim(pff2016)[1], big.mark = ",")

# More details for 2016 peak counts
sites <- unique(pff2016$site_no)

mycounts <- pff2016 %>% group_by(site_no) %>% tally()

paste("There are ", prettyNum(dim(pff2016)[1], big.mark = ","), 
      "entries of peak streamflow for", 
      prettyNum(length(sites), big.mark = ","), 
      "sites with periods of record from 1 to", 
      max(mycounts$n), "years")

PercentOf2008Peaks <- function(dataset) {
  percent <- formatC(dim(dataset)[1] / dim(pff2008)[1] * 100, 
                     digits = 2, drop0trailing = FALSE, format = "f")
  
  message(paste("There were", prettyNum(dim(dataset)[1], big.mark = ","), 
                "changes representing", percent, "percent of the 2008 peaks."))
}

# percent of 2008 peaks with changed peak value
cpv <- subset(allpkdat, peak_va.x != peak_va.y)
PercentOf2008Peaks(cpv)

# percent of 2008 peaks where 2008 is larger than 2016
cpv1 <- subset(cpv, peak_va.x > peak_va.y)
PercentOf2008Peaks(cpv1)

# percent of 2008 peaks where 2008 is larger than 2016
cpv2 <- subset(cpv, peak_va.x < peak_va.y)
PercentOf2008Peaks(cpv2)

#  changed peak date and percent calculation
cpd <- subset(allpkdat, peak_dt.x != peak_dt.y)
PercentOf2008Peaks(cpd)

# changed peak qualification code 
cpqc <- subset(allpkdat, peak_cd.x != peak_cd.y)
PercentOf2008Peaks(cpqc)

# changed peak gage height
cgh <- subset(allpkdat, (is.na(gage_ht.x) & !is.na(gage_ht.y)) | 
                (!is.na(gage_ht.x) & is.na(gage_ht.y)) | 
                (gage_ht.x != gage_ht.y))
PercentOf2008Peaks(cgh)

# changed peak gage-height qualification code 
cghqc <- subset(allpkdat, gage_ht_cd.x != gage_ht_cd.y)
PercentOf2008Peaks(cghqc)

# Obtain latitude, longitude, and station name

sites <- unique(allpkdat$site_no)
siteInfo <- data.frame(site_no = sites, name = "Site name",
                      lat = -9999, lng = -9999, 
                      stringsAsFactors = FALSE)

# update to readNWISsite function in dataRetrieval
# which was failing because it did not use https
readNWISsite <- function (siteNumbers) 
{
  siteNumber <- paste(siteNumbers, collapse = ",")
  urlSitefile <- paste0("https://waterservices.usgs.gov/nwis/site/?format=rdb&siteOutput=Expanded&sites=", 
                        siteNumber)
  data <- importRDB1(urlSitefile, asDateTime = FALSE)
  data[, grep("_va", names(data))] <- sapply(data[, grep("_va", 
                                                         names(data))], as.numeric)
  return(data)
}

# look up site information in groups of 100
# much faster than 1 at a time,
# but process can fail if all are attempted at once
groups <- seq(from = 1, to = length(sites), by = 100)
groups <- c(groups, length(sites))
i <- 1
staInfo <- readNWISsite(sites[i])
siteInfo$lat[i] <- staInfo$dec_lat_va
siteInfo$lng[i] <- staInfo$dec_long_va
siteInfo$name[i] <- staInfo$station_nm
   
for (i in 2:length(groups)) {
   staInfo <- readNWISsite(sites[(groups[(i - 1)] + 1):groups[i]])
   siteInfo$lat[(groups[(i - 1)] + 1):groups[i]] <- staInfo$dec_lat_va
   siteInfo$lng[(groups[(i - 1)] + 1):groups[i]] <- staInfo$dec_long_va
   siteInfo$name[(groups[(i - 1)] + 1):groups[i]] <- staInfo$station_nm
}


allpkdat <- left_join(allpkdat, siteInfo, by = "site_no")

# incorporate percentage difference map changes

allpkdat$Percent.Diff <- -999999
pck <- allpkdat$peak_va.x == allpkdat$peak_va.y
allpkdat$Percent.Diff[pck] <- 0
pck <- allpkdat$peak_va.x == 0 & allpkdat$peak_va.y != allpkdat$peak_va.x
allpkdat$Percent.Diff[pck] <- NA
pck <- allpkdat$Percent.Diff < 0 & !is.na(allpkdat$Percent.Diff)
allpkdat$Percent.Diff[pck] <- (allpkdat$peak_va.y[pck] - 
                                 allpkdat$peak_va.x[pck]) / allpkdat$peak_va.x[pck] * 100

NationalPeakChangeMap <- function(dataset, mytitle, ptcol = "blue", ptpch = 16, 
                                  ptbg = "grey50") {
  require(maps)
  require(mapproj)
  proj.type <- "lambert"
  proj.orient <- c(90, 0, -100)
  par(mar = c(0.5, 0.5, 0, 0))
  map("world", regions = c("usa", "canada", "mexico", "puerto rico", "cuba", "haiti",
                           "dominican republic"), col = 1, project = proj.type, 
      param = c(20, 50), orient = proj.orient,
      xlim = c(-160, -55), ylim = c(17, 75))
  map("state", lty = 1, col = "grey60", project = proj.type, param = c(20, 50),
      orient = proj.orient, add = TRUE)
  map("usa", interior = FALSE, col = 1, xlim = c(-160, -55), ylim = c(17.5, 75),
      add = TRUE, project = "")
  points(mapproject(dataset$lng, dataset$lat), pch = ptpch, col = ptcol, 
         cex = 0.9, bg = ptbg)
  numpks <- prettyNum(dim(dataset)[[1]], big.mark = ",")
  numsites <- prettyNum(length(unique(dataset$site_no)), big.mark = ",")
  
  myleg <- paste(numpks, " peaks at ", numsites, " sites\n", mytitle, sep = "")
  
  legend(mapproject(-150, 42), myleg, pch = c(ptpch), col = c(ptcol), 
         pt.bg = ptbg, bty = "n",
         yjust = 1, xjust = 0.5, cex = 0.9)
  mtext(side = 1, 
        "Base from maps package for R (Becker and others, 2016)\nLambert projection", 
        adj = 0, cex = 0.5)
}

NationalPeakChangeTable <- function(dataset, mytitle) {
  numpks <- prettyNum(dim(dataset)[[1]], big.mark = ",")
  numsites <- prettyNum(length(unique(dataset$site_no)), big.mark = ",")
  myleg <- paste(numpks, " peaks at ", numsites, " sites\n", mytitle, sep = "")
  
  if (dim(dataset)[[1]] <= 30 & dim(dataset)[[1]] > 0) {
    res.df <- dataset[, c(1, 13, 2:12)]
    colnames(res.df) <- c("Site", "Site name", "Water year", "Date in 2008", 
                          "Peak in 2008", "Code in 2008", "Gage height in 2008", 
                          "Gage height code in 2008", "Date in 2016", "Peak in 2016", 
                          "Code in 2016", "Gage height in 2016", 
                          "Gage height code in 2016")
    knitr::kable(res.df, caption = myleg, row.names = FALSE, 
                 format.args = list(big.mark = ","))
  }
}

CountChange <- function(data, percent1, percent2, lt = FALSE, mytitle = "") {
  if (lt == FALSE ) {
    dat <- subset(allpkdat, Percent.Diff > percent1 & Percent.Diff <= percent2)
    NationalPeakChangeMap(dat, mytitle, ptcol = "firebrick1", ptpch = 17)
    NationalPeakChangeTable(dat, mytitle)
  } else {
    dat <- subset(allpkdat, Percent.Diff < percent1 & Percent.Diff >= percent2)
    NationalPeakChangeMap(dat, mytitle, ptcol = "firebrick1", ptpch = 25, 
                          ptbg = "firebrick1")
    NationalPeakChangeTable(dat, mytitle)
  }
}

PlotTimeLine <- function(x, y, filename, mycol, ys, ylabs) {
  fn <- paste(filename, ".pdf", sep = "")
  pdf(fn, width = 8, height = 3.5, pointsize = 9)
  par(las = 1, tck = 0.02, cex = 0.8, mar = c(4, 4, 1, 1) + 0.1, mfrow = c(1, 1))
  plot(x, y, type = "n", axes = FALSE, 
       xlab = "Water Year", ylab = "Count", ylim = ys, 
       xlim = c(1780, 2020), xaxs = "i", yaxs = "i")
  abline(h = ylabs[-1], col = "grey80", lty = 3)
  yrs <- seq(from = 1780, to = 2020, by = 20)
  points(x, y, type = "h", col = "grey60")
  points(x, y, pch = 18, col = mycol, 
         cex = 0.9)
  axis(1, at = yrs)
  axis(3, at = yrs, labels = FALSE)
  axis(2, at = ylabs, labels = prettyNum(ylabs, big.mark = ","))
  axis(4, at = ylabs, labels = FALSE)
  dev.off()
}

# Figures 4-11

pdf("Fig04SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 0, 5, 
            mytitle = "Percent change in \n peak value > 0 and <= 5")
CountChange(allpkdat, 0, -5, lt = TRUE, 
            mytitle = "Percent change in \n peak value < 0 and >= -5")
dev.off()

pdf("Fig05SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 5, 10, 
            mytitle = "Percent change in \n peak value > 5 and <= 10")
CountChange(allpkdat, -5, -10, lt = TRUE, 
            mytitle = "Percent change in \n peak value < -5 and >= -10")
dev.off()

pdf("Fig06SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 10, 25, 
            mytitle = "Percent change in \n peak value > 10 and <= 25")
CountChange(allpkdat, -10, -25, lt = TRUE, 
            mytitle = "Percent change in \n peak value < -10 and >= -25")
dev.off()

pdf("Fig07SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 25, 50, 
            mytitle = "Percent change in \n peak value > 25 and <= 50")
CountChange(allpkdat, -25, -50, lt = TRUE, 
            mytitle = "Percent change in \n peak value < -25 and >= -50")
dev.off()

pdf("Fig08SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 50, 100, 
            mytitle = "Percent change in \n peak value > 50 and <= 100")
CountChange(allpkdat, -50, -100, lt = TRUE, 
            mytitle = "Percent change in \n peak value < -50 and >= -100")
dev.off()

pdf("Fig09SitesWithChangeToPeakValue.pdf", onefile = FALSE, width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 100, 500, 
            mytitle = "Percent change in \n peak value > 100 and <= 500")
CountChange(allpkdat, 500, 1000, 
            mytitle = "Percent change in \n peak value > 500 and <= 1,000")
dev.off()

pdf("Fig10SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 1000, 5000, 
            mytitle = "Percent change in \n peak value > 1,000 \n and <= 5,000")
CountChange(allpkdat, 5000, 10000, 
            mytitle = "Percent change in \n peak value > 5,000 \n and <= 10,000")
dev.off()

pdf("Fig11SitesWithChangeToPeakValue.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountChange(allpkdat, 10000, 50000, 
            mytitle = "Percent change in \n peak value > 10,000 \n and <= 50,000")
CountChange(allpkdat, 50000, Inf, 
            mytitle = "Percent change in \n peak value > 50,000")
dev.off()

PkChng <- subset(allpkdat, peak_va.x != peak_va.y)
PkChngCount <- PkChng %>% group_by(WaterYr) %>% tally()

# Figure 12
PlotTimeLine(PkChngCount$WaterYr, PkChngCount$n, "Fig12TimeLinePkChng",
             mycol = "firebrick1", ys = c(-1, 200), 
             ylabs = c(0, 50, 100, 150, 200))

NationalCodeChangeTable <- function(dataset, mytitle) {
  numpks <- prettyNum(dim(dataset)[[1]], big.mark = ",")
  numsites <- prettyNum(length(unique(dataset$site_no)), big.mark = ",")
  myleg <- paste(numpks, " peaks at ", numsites, " sites \n", mytitle, sep = "")
  if (dim(dataset)[[1]] <= 30 & dim(dataset)[[1]] > 0) {
    res.df <- dataset[, c(1, 13, 2:12)]
    colnames(res.df) <- c("Site", "Site name", "Water Year", "Date in 2008", 
                          "Peak in 2008", "Code in 2008", "Gage height in 2008", 
                          "Gage height code in 2008", "Date in 2016", 
                          "Peak in 2016", "Code in 2016", "Gage height in 2016", 
                          "Gage height code in 2016")
    knitr::kable(res.df, caption = myleg, row.names = FALSE)
  }
}

CountCodeChange1 <- function(data, code) {
  dat <- data[grepl(code, data$peak_cd.x) & !grepl(code, data$peak_cd.y), ]
  mytitle <- paste(" with code", code, "removed")
  NationalPeakChangeMap(dat, mytitle, ptcol = "purple", ptpch = 15)
  NationalCodeChangeTable(dat, mytitle)
}

CountCodeChange2 <- function(data, code) {
  dat <- data[!grepl(code, data$peak_cd.x) & grepl(code, data$peak_cd.y), ]
  mytitle <- paste(" with code", code, "added")
  NationalPeakChangeMap(dat, mytitle, ptcol = "darkgreen", ptpch = 18)
  NationalCodeChangeTable(dat, mytitle)
}

# Figures 13-26

pdf("Fig13SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "1")
CountCodeChange2(allpkdat, "1")
dev.off()

pdf("Fig14SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "2")
CountCodeChange2(allpkdat, "2")
dev.off()

pdf("Fig15SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "3")
CountCodeChange2(allpkdat, "3")
dev.off()

pdf("Fig16SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "4")
CountCodeChange2(allpkdat, "4")
dev.off()

pdf("Fig17SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "5")
CountCodeChange2(allpkdat, "5")
dev.off()

pdf("Fig18SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "6")
CountCodeChange2(allpkdat, "6")
dev.off()

pdf("Fig19SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "7")
CountCodeChange2(allpkdat, "7")
dev.off()

pdf("Fig20SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "8")
CountCodeChange2(allpkdat, "8")
dev.off()

pdf("Fig21SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "9")
CountCodeChange2(allpkdat, "9")
dev.off()

pdf("Fig22SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "A")
CountCodeChange2(allpkdat, "A")
dev.off()

pdf("Fig23SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "B")
CountCodeChange2(allpkdat, "B")
dev.off()

pdf("Fig24SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "C")
CountCodeChange2(allpkdat, "C")
dev.off()

pdf("Fig25SitesWithChangeToPeakCode.pdf",  width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "D")
CountCodeChange2(allpkdat, "D")
dev.off()

pdf("Fig26SitesWithChangeToPeakCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountCodeChange1(allpkdat, "E")
CountCodeChange2(allpkdat, "E")
dev.off()

PkCdChng <- subset(allpkdat, peak_cd.x != peak_cd.y)
PkCdChngCount <- PkCdChng %>% group_by(WaterYr) %>% tally()

# Figure 27

PlotTimeLine(PkCdChngCount$WaterYr, PkCdChngCount$n, "Fig27TimeLinePkCdChng",
             mycol = "darkgreen", ys = c(-2, 800), 
             ylabs = c(0, 200, 400, 600, 800))

CountPkDt <- function(data) {
  dat <- subset(data, peak_dt.x != peak_dt.y)
  mytitle <- " with changes to peak date"
  NationalPeakChangeMap(dat, mytitle)
  
  dat <- subset(data, peak_dt.x != peak_dt.y & peak_cd.x != peak_cd.y)
  mytitle <- " with changes to peak date \n and peak qualification code"
  NationalPeakChangeMap(dat, mytitle)
}

# Figure 28

pdf("Fig28SitesWithChanges2PkDt.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountPkDt(allpkdat)
dev.off()

PkDtChng <- subset(allpkdat, peak_dt.x != peak_dt.y)
PkDtChngCount <- PkDtChng %>% group_by(WaterYr) %>% tally()

# Figure 29

PlotTimeLine(PkDtChngCount$WaterYr, PkDtChngCount$n, "Fig29TimeLinePkDtChng",
             mycol = "blue", ys = c(0, 50), ylabs = c(0, 10, 20, 30, 40, 50))

CountPkValPkDtPkCode <- function(data) {
  dat <- subset(data, peak_va.x != peak_va.y & peak_dt.x != peak_dt.y & 
                  peak_cd.x != peak_cd.y)
  mytitle <- " with changes to peak value, \n peak date, and peak qualification code"
  NationalPeakChangeMap(dat, mytitle)
}

# Figure 30 

pdf("Fig30SitesWithChanges2PkValPkDtPkCd.pdf", width = 4.5, height = 4.25, pointsize = 7)
CountPkValPkDtPkCode(allpkdat)
dev.off()

CountGageHeightAddedRemoved <- function(data) {
  dat <- subset(data, is.na(gage_ht.x) & !is.na(gage_ht.y))
  mytitle <- " with gage height added"
  NationalPeakChangeMap(dat, mytitle, ptcol = "darkgreen", ptpch = 23, 
                        ptbg = "darkgreen")
  
  dat <- subset(data, !is.na(gage_ht.x) & is.na(gage_ht.y))
  mytitle <- " with gage height removed"
  NationalPeakChangeMap(dat, mytitle, ptcol = "purple", ptpch = 15)
}

# Figure 31

pdf("Fig31SitesWithChangeToGageHeight.pdf", width = 4.5, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountGageHeightAddedRemoved(allpkdat)
dev.off()

CountGHValGHCode <- function(data) {
  dat <- subset(data, !is.na(gage_ht.x) & !is.na(gage_ht.y) & 
                  gage_ht.x != gage_ht.y)
  mytitle <- " for which gage height existed \n in 2008 and in 2016, but differed"
  NationalPeakChangeMap(dat, mytitle)
  
  dat <- subset(data, gage_ht.x != gage_ht.y & gage_ht_cd.x != gage_ht_cd.y)
  mytitle <- " with changes to gage height \n and gage-height qualification code"
  NationalPeakChangeMap(dat, mytitle)
}

# Figure 32

pdf("Fig32SitesWithChanges2GHValGHCd.pdf", width = 4.5, height = 4.25, pointsize = 7)
par(mfrow = c(1, 2))
CountGHValGHCode(allpkdat)
dev.off()

GgHtChng <- subset(allpkdat, (is.na(gage_ht.x) & !is.na(gage_ht.y)) | 
                     (is.na(gage_ht.y) & !is.na(gage_ht.x))  | gage_ht.x != gage_ht.y)
GgHtChngCount <- GgHtChng %>% group_by(WaterYr) %>% tally()

# Figure 33

PlotTimeLine(GgHtChngCount$WaterYr, GgHtChngCount$n, "Fig33TimeLineGgHtChng",
             mycol = "purple", ys = c(-1, 400), ylabs = c(0, 100, 200, 300, 400))

CountGageHeightCodeChange <- function(data, code) {
  dat <- data[grepl(code, data$gage_ht_cd.x) & !grepl(code, data$gage_ht_cd.y), ]
  mytitle <- paste(" with gage-height qualification \n code", code, "removed")
  NationalPeakChangeMap(dat, mytitle, ptcol = "purple", ptpch = 15)
  
  dat <- data[!grepl(code, data$gage_ht_cd.x) & grepl(code, data$gage_ht_cd.y), ]
  mytitle <- paste(" with gage-height qualification \n code", code, "added")
  NationalPeakChangeMap(dat, mytitle, ptcol = "darkgreen", ptpch = 23, 
                        ptbg = "darkgreen")
}

# Figures 34-39

pdf("Fig34SitesWithChangeToGageHeightCode.pdf", width = 9, 
    height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "1")
dev.off()

pdf("Fig35SitesWithChangeToGageHeightCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "2")
dev.off()

pdf("Fig36SitesWithChangeToGageHeightCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "3")
dev.off()

pdf("Fig37SitesWithChangeToGageHeightCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "4")
dev.off()

pdf("Fig38SitesWithChangeToGageHeightCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "5")
dev.off()

pdf("Fig39SitesWithChangeToGageHeightCode.pdf", width = 9, height = 4.25, pointsize = 7)
par(mfrow = c(1,2))
CountGageHeightCodeChange(allpkdat, "6")
dev.off()


GgHtCdChng <- subset(allpkdat, gage_ht_cd.x != gage_ht_cd.y)
GgHtCdChngCount <- GgHtCdChng %>% group_by(WaterYr) %>% tally()

# Figure 40

PlotTimeLine(GgHtCdChngCount$WaterYr, GgHtCdChngCount$n, "Fig40TimeLineGgHtCdChng",
             mycol = "darkgreen", ys = c(-1, 400), ylabs = c(0, 100, 200, 300, 400))

AllPkChng <- subset(allpkdat, (peak_va.x != peak_va.y) | (is.na(gage_ht.x) & !is.na(gage_ht.y)) | 
                      (!is.na(gage_ht.x) & is.na(gage_ht.y)) | (gage_ht.x != gage_ht.y) | 
                      (peak_dt.x != peak_dt.y) | (peak_cd.x != peak_cd.y) | 
                      (gage_ht_cd.x != gage_ht_cd.y))
AllPkChngCount <- AllPkChng %>% group_by(WaterYr) %>% tally()

# Figure 41

PlotTimeLine(AllPkChngCount$WaterYr, AllPkChngCount$n, "Fig41TimeLineAllPkChng", 
             mycol = "lightsalmon2", ys = c(-2, 1200), 
             ylabs = c(0, 200, 400, 600, 800, 1000, 1200))


cat("\n", format(Sys.time(), "%a %d %b %Y %X"))
print(sessionInfo(), locale = TRUE)