Stream Discharge Visualization Utilizing Tidyverse

library(dataRetrieval)
library(ggplot2)
library(scales)
library(lubridate)
library(tidyverse)

# Log scale tick marks
logTicks <- function(n = 5, base = 10){
  divisors <- which((base / seq_len(base)) %% 1 == 0)
  mkTcks <- function(min, max, NumBase, divisor){
    f <- seq(divisor, base, by = divisor)
    return(unique(c(base^min, as.vector(outer(f, base^(min:max), `*`)))))
  }
  function(x) {
    rng <- range(x, na.rm = TRUE)
    lrng <- log(rng, base = base)
    min <- floor(lrng[1])
    max <- ceiling(lrng[2])
    tck <- function(divisor){
      t <- mkTcks(min, max, base, divisor)
      t[t >= rng[1] & t <= rng[2]]
    }
    tcks <- lapply(divisors, function(d) tck(d))
    l <- vapply(tcks, length, numeric(1))
    i <- which.min(abs(n - l))
    if(l[i] < 2){
      ticks <- pretty(x, n = n, min.n = 2)
    } else {
      ticks <- tcks[[i]]
    }
    return(ticks)
  }
}

# Custom ggplot theme
themebox <- function(base_family = "sans", ...){
  theme_bw(base_family = base_family, ...) +
    theme(
      plot.title = element_text(size = 18, face = "bold", hjust = 0.5),
      panel.background = element_rect(fill = "gray95", colour = "black"),
      panel.grid.major = element_line(size = 0.95, linetype = 'solid', colour = "white"),
      panel.grid.minor = element_line(size = 0.6, linetype = 'solid', colour = "white"),
      axis.ticks.length = unit(.15, "cm"),
      axis.ticks = element_line(colour = 'black', size = .60),
      axis.title = element_text(face = "bold", size = 18),
      axis.text.x = element_text(margin = unit(c(0.3,0.3,0.3,0.3), "cm"), colour = "black", size = 12),
      axis.text.y = element_text(margin = unit(c(0.3,0.3,0.3,0.3), "cm"), colour = "black", size = 18),
      axis.line = element_line(size = 0.5, colour = "black"),
      legend.position = "bottom",
      legend.box = "horizontal",
      legend.box.margin = margin(-15, -10, -1, -10),
      legend.title = element_text(size = 9),
      plot.margin = margin(0.25, 1, 0.01, 0.1, "cm")
    )
}

# Setup
Today <- as.Date(Sys.Date())
Year <- format(Today, "%Y")

Dashboard <- data.frame(
  site_no = c('01463500', '01438500', '01474500', '01481000'),
  site_name = c('TRENTON, NJ', 'MONTAGUE, NJ', 'PHILADELPHIA, PA', 'CHADDS FORD, PA'),
  river_name = c('DELAWARE RIVER', 'DELAWARE RIVER', 'SCHUYLKILL RIVER', 'BRANDYWINE CREEK'),
  file_name = c('Trenton', 'Montague', 'Schuylkill', 'ChaddsFord')
)

# Loop through each site
for (i in 1:nrow(Dashboard)) {
  
  png(file = paste0(Dashboard$file_name[i], ".png"), width = 4000, height = 2000, res = 300)
  
  Stats <- tryCatch({
    readNWISstat(
      siteNumbers = Dashboard$site_no[i],
      parameterCd = "00060",
      statReportType = "daily",
      statType = c("P10", "P25", "P50", "P75")
    )
  }, error = function(e) {
    message(paste("Error retrieving stats for site", Dashboard$site_no[i]))
    dev.off()
    next
  })
  
  if (nrow(Stats) == 0) {
    message(paste("No stats data for site", Dashboard$site_no[i]))
    dev.off()
    next
  }
  
  Stats$MonthDay <- with(Stats, paste(month_nu, day_nu, sep = "-"))
  
  DateFill <- seq(Today - 21, Today, by = "day")
  NewTimeSeries <- data.frame(DateFill) %>%
    mutate(
      Year = year(DateFill),
      Month = month(DateFill),
      Day = mday(DateFill),
      MonthDay = paste(Month, Day, sep = "-")
    )
  
  Update <- right_join(Stats, NewTimeSeries, by = "MonthDay")
  Update$Ymd <- as.Date(with(Update, paste(Year, Month, Day, sep = "-")), "%Y-%m-%d")
  Update$dateTime <- as.POSIXct(Update$Ymd, tz = "UTC")
  
  NOW <- tryCatch({
    readNWISuv(
      siteNumbers = Dashboard$site_no[i],
      parameterCd = "00060",
      startDate = Today - 21,
      endDate = Today
    )
  }, error = function(e) {
    message(paste("Error retrieving current data for site", Dashboard$site_no[i]))
    dev.off()
    next
  })
  
  OUTPUT <- right_join(Update, NOW, by = "dateTime") %>%
    select(21, 23, 24, 10:14)
  
  colnames(OUTPUT) <- c("dateTime", "site_no", "X_00060_00000", "count_nu", "p10_va", "p25_va", "p50_va", "p75_va")
  
  data_long <- gather(OUTPUT, PercentileName, PercentileValue, p10_va:p75_va, factor_key = TRUE)
  
  All <- ggplot(data_long, aes(x = dateTime, y = X_00060_00000)) +
    geom_point(aes(y = PercentileValue, color = factor(PercentileName)), size = 5, shape = 17) +
    stat_smooth(aes(y = X_00060_00000, fill = 'Discharge'), method = "lm",
                formula = y ~ poly(x, 21), se = FALSE, size = 2, color = "black") +
    scale_y_log10(
      name = 'Discharge (cubic feet per second)',
      breaks = logTicks(n = 4),
      minor_breaks = logTicks(n = 10)
    ) +
    scale_x_datetime(
      breaks = date_breaks("2 day"),
      labels = date_format("%b\n%d\n%Y"),
      name = NULL,
      expand = c(0, 0)
    ) +
    ggtitle(paste0("USGS ", Dashboard$site_no[i], " ", Dashboard$river_name[i], " AT ", Dashboard$site_name[i])) +
    themebox() +
    labs(color = paste0("(", mean(data_long$count_nu, na.rm = TRUE), " years)")) +
    scale_fill_manual("", values = c("#171A21")) +
    scale_color_manual(
      labels = c("P10", "P25", "P50", "P75"),
      values = c("#FA2121", "#FCD600", "#01A558", "#00AFEA")
    ) +
    guides(fill = guide_legend(order = 1), colour = guide_legend(title.position = "right"))
  
  print(All)
  dev.off()
}

graphics.off()

tren.mont.sch.chadds.dashboard.zip
File Size:	1606 kb
File Type:	zip

Download File

---