# Imports
library(dplyr)
library(tidyhydat)
library(purrr)

calc_thresh <- function(wsc_station_number){

  # Add tick to progress bar
  pb$tick()$print()

  # Grab data for station (all parameters)
  station_data <- hy_daily(
    station_number = wsc_station_number
  ) %>%
    # Group by parameter type for stat calculations
    group_by(Parameter) %>%
    # Get rid of missing values
    filter(!is.na(Value)) %>%
    summarise(
      "station_number" = wsc_station_number,
      # Number of values considered
      "n" = n(),
      # Percentiles
      "25%" = quantile(Value, 0.25),
      "50%" = quantile(Value, 0.5),
      "75%" = quantile(Value, 0.75),
      "90%" = quantile(Value, 0.9),
      "100%" = quantile(Value, 1)
    ) %>%
    mutate_if(is.double, round, 3) %>%
    ungroup()

  return(station_data)
}

# All active stations
active_ont_stns <- hy_stations(
  prov_terr_state_loc = "ON"
) %>%
  filter(
    HYD_STATUS == "ACTIVE"
  )

# Set progress bar
pb <- progress_estimated(length(active_ont_stns$STATION_NUMBER))

# Grab thresholds for all stations
all_stn_thresh <- purrr::map(
  active_ont_stns$STATION_NUMBER,
  calc_thresh
) %>%
  bind_rows() %>%
  # Only flow/level stats
  filter(
    Parameter %in% c("Flow", "Level")
  ) %>%
  # Rearrange columns
  select(
    station_number,
    everything()
  )