In USGS-R/HASP: Hydrologic AnalySis Package

library(knitr)
opts_chunk$set(echo = FALSE,
               warning = FALSE,
               message = FALSE,
               fig.width = 7,
               fig.height = 6)
library(dplyr)
library(tidyr)
library(ggplot2)
library(HASP)

options(knitr.kable.NA = '')

Site Information

library(HASP)

siteID <- "{{{ site }}}"
type <- "{{{ output_type }}}"

site_metadata <- site_summary(siteID, markdown = type == "html")

pcodes <- dataRetrieval::whatNWISdata(siteNumber = siteID,
                                      service = c("gw", "dv"))

data_info <- data_available(siteID)
kable(data_info)

library(dataRetrieval)

# Figure out what data to retrieve
available_pcodes <- unique(pcodes$parm_cd[pcodes$data_type_cd == "dv"])
available_statCds <- unique(pcodes$stat_cd[pcodes$data_type_cd == "dv"])

# Use the parameter with the highest count
# This can be changed and then r
parameterCd <- unique(pcodes$parm_cd[which(pcodes$count_nu ==
                                      max(pcodes$count_nu,na.rm = TRUE) & 
                                      pcodes$data_type_cd == "dv")])

statCd <- pcodes$stat_cd[which(pcodes$count_nu ==
                                      max(pcodes$count_nu,na.rm = TRUE) & 
                                      pcodes$data_type_cd == "dv")][1]

if(length(parameterCd) == 0){
  parameterCd <- unique(pcodes$parm_cd[which(pcodes$count_nu ==
                                      max(pcodes$count_nu,na.rm = TRUE) & 
                                      pcodes$data_type_cd == "gw")])[1]
}

pCode_info <- dataRetrieval::readNWISpCode(parameterCd)
yaxis <- pCode_info$parameter_nm

flip <- parameterCd == "72019"

# Daily data:
if(any(pcodes$data_type_cd == "dv")){
  gw_level_dv <- dataRetrieval::readNWISdv(siteID,
                            available_pcodes,
                            statCd = available_statCds)
} else {
  gw_level_dv <-  L2701_example_data$Discrete
  gw_level_dv <- gw_level_dv[-1:-nrow(gw_level_dv),]
}
# Field GWL measured:
if(any(pcodes$data_type_cd == "gw")){
  gwl_data <- dataRetrieval::readNWISgwl(siteID)
} else {
  gwl_data <-  L2701_example_data$Discrete
  gwl_data <- gwl_data[-1:-nrow(gwl_data),]
}

Water Level Data and Analysis

Monthly frequency

monthly_frequency_plot(gw_level_dv,
                       gwl_data,
                       parameter_cd = parameterCd,
                       stat_cd = statCd, 
                       plot_title = siteID, 
                       flip = flip,
                       y_axis_label = yaxis)

month_table <- monthly_frequency_table(gw_level_dv,
                                       gwl_data,
                                       parameter_cd = parameterCd,
                                       stat_cd = statCd,
                                       flip = flip)

if(type == "html"){
  knitr::kable(month_table, digits = 1) %>%
      kableExtra::kable_styling()%>%
      kableExtra::row_spec(as.integer(strftime(Sys.Date(), format = "%m")),
                           bold = TRUE,
                           hline_after = TRUE)
} else {
  knitr::kable(month_table, digits = 1)
}

Weekly frequency

weekly_frequency_plot(gw_level_dv, 
                      gwl_data,
                      parameter_cd = parameterCd, 
                      stat_cd = statCd,
                      plot_title = siteID,
                      flip = flip,
                      y_axis_label = yaxis)

Weekly frequency analysis of daily maximum water level record. Only showing the first 10 rows for this example:

weekly_table <- weekly_frequency_table(gw_level_dv, 
                                       gwl_data, 
                                       parameter_cd = parameterCd, 
                                       stat_cd = statCd,
                                       flip = flip) 

todays_week <- as.integer(strftime(Sys.Date(), format = "%V"))
centered_index <- (todays_week-4):(todays_week+5)
if(any(centered_index <= 0)){
  centered_index[centered_index <= 0] <- centered_index[centered_index <= 0] + nrow(weekly_table)
}
if(type == "html"){
  knitr::kable(weekly_table[centered_index,], digits = 1) %>%
      kableExtra::kable_styling()%>%
      kableExtra::row_spec(5,
                           bold = TRUE,
                           hline_after = TRUE)
} else {
  knitr::kable(weekly_table[centered_index,], digits = 1)
}

Daily 2-year

if(nrow(gw_level_dv) > 0){
  daily_gwl_plot(gw_level_dv, 
                 gwl_data,
                 parameter_cd = parameterCd,
                 plot_title = siteID,
                 stat_cd = statCd,
                 flip = flip,
                 historical_stat = "mean",
                 month_breaks = TRUE,
                 y_axis_label = yaxis)
} else {
  cat("Daily plot only useful for sites with daily data.")
}

Statistics of maximum daily water level record (DOY = day of year). Only showing the first 10 rows for this example:

if(nrow(gw_level_dv) > 0){
  daily_table <- daily_frequency_table(gw_level_dv, 
                                       NULL,
                                       parameter_cd = parameterCd,
                                       stat_cd = statCd) 

  todays_day <- as.POSIXlt(Sys.Date())$yday
  centered_day <- (todays_day-4):(todays_day + 5)

  if(any(centered_day <= 0)){
    centered_day[centered_day <= 0] <- centered_day[centered_day <= 0] + nrow(daily_table)
  }

  if(type == "html"){
    knitr::kable(daily_table[centered_day,]) %>%
        kableExtra::kable_styling()%>%
        kableExtra::row_spec(5,
                             bold = TRUE,
                             hline_after = TRUE)
  } else {
    knitr::kable(daily_table[centered_day,])
  }
}

Daily value trends

The following plot is based on daily data:

gwl_plot_all(gw_level_dv, 
             NULL, 
             y_label = yaxis,
             parameter_cd = parameterCd, 
             stat_cd = statCd,
             flip = flip,
             plot_title = siteID,
             add_trend = TRUE)

Summary statistics for maximum daily water level measurements:

if(nrow(gw_level_dv) > 0){

  siteDV <- daily_gwl_summary(gw_level_dv,
                              NULL,
                              parameter_cd = parameterCd,
                              stat_cd = statCd) 

  site_out <- knitr::kable(t(siteDV), digits = 1)

  if(type == "html"){
    site_out <- site_out %>%
    kableExtra::kable_styling(full_width = FALSE)
  }

  site_out
} else {
  cat("No daily data to aggregate.")
}

Results of trend analysis on maximum daily water levels:

if(nrow(gw_level_dv) > 0){
    data_list <- HASP:::set_up_data(gw_level_dv = gw_level_dv, 
                           gwl_data = NULL, 
                           parameter_cd = parameterCd,
                           date_col = NA, 
                           value_col = NA,
                           approved_col = NA,
                           stat_cd = statCd)

  gw_level_dv1 <- data_list$gw_level_dv
  gwl_data1 <- data_list$gwl_data

  all_data <- gw_level_dv1 %>%
    dplyr::bind_rows(gwl_data1) %>%
    dplyr::filter(grepl("A", Approve))

  trend_results <- trend_test(all_data,
                                gwl_data = NULL,
                                date_col = c("Date"),
                                value_col = c("Value"), 
                                approved_col = c("Approve"))

  kable(trend_results, digits = 1)
} else {
  cat("No daily data to aggregate.")
}

Field GWL values

gwl_plot_all(NULL, gwl_data,
             y_label = yaxis,
             flip = flip,
             parameter_cd = parameterCd,
             plot_title = siteID)

Summary statistics for manual water level measurements

quantiles <- gwl_data %>% 
  rename(value = sl_lev_va) %>% 
  site_data_summary() 

q_out <- knitr::kable(t(quantiles), digits = 1)

if(type == "html"){
  q_out <- q_out %>%
  kableExtra::kable_styling(full_width = FALSE)
}

q_out