BBN_data.Rmd

Bayesian Belief Network: Hydrological Datasets

Flow and water level data as well as IHA variables and ice cycle variables. 

Run "Eflows_FUNCTIONS.R" file before running this script

```{r}

library(IHA)
library(tidyverse)
library(tidyhydat)
library(zoo)
library(lubridate)
library(ggplot2)
library(dataRetrieval)
library(knitr)

```

##FLOW DATA

#1 extract flow data and add column for water year and day of the year

```{r}

#get stns to analyze
stn_all <- tidyhydat::hy_daily_flows(station_number = c('07HA001')) 


#write.csv(stn_all, "PeaceRiver_PR_07HA001_RawFlowData_1915-2020.csv")

stn_all <- stn_all %>%
group_by(STATION_NUMBER) %>%
#make a complete set of days for each year 
complete(Date = seq.Date(as.Date("1989/10/1"), as.Date("2021/09/30"), by="day"))

#add water year
stn_all <- addWaterYear(stn_all) %>%
mutate(waterYear = as.character(waterYear))

#add day of the year from "Eflows_FUNCTIONS.R" function 
stn_all_split <- split(stn_all, stn_all$STATION_NUMBER )
stn_all_split_doy <- lapply(stn_all_split, calc_day_of_wyear)
stn_all <- bind_rows(stn_all_split_doy, .id = "STATION_NUMBER")

rm(stn_all_split, stn_all_split_doy)

#other station options:
#hy_stations(station_number = '07DA001')
#peace point: 07KC001
#slave river at fitzgerald: 07NB001
#birch river below alice creek: 07KE001
#athabasca at ft mcmurry: 07DA001

```

#Export raw flow data to csvs
Note: these csvs are untidied flow data. No estimates, no deleting years with >14 days of missing data, etc. See below for tidied data

```{r}
#before do this make sure to get rid of years we don't want, this wasn't done above
#df1 <- stn_all %>% filter(STATION_NUMBER == "07KC001")
#df2 <- stn_all %>% filter(STATION_NUMBER == "07NB001 ")
#df3 <- stn_all %>% filter(STATION_NUMBER == "07KE001 ")
#df4 <- stn_all %>% filter(STATION_NUMBER == "07DA001 ")

#write.csv(df1, "PeaceRiver_07KC001_Flow.csv")
#write.csv(df2, "SlaveRiver_07NB001_Flow.csv")
#write.csv(df3, "BirchRiver_07KE001_Flow.csv")
#write.csv(stn_all, "AthabascaRiver_07DA001_RawFlowData.csv")

```

#2 Make three separate datasets in prep for IHA calculations: 
-full year
-open water data (April-October) 
-winter/closed water (November-March)

```{r}
stn_all <- stn_all # full dataset

#open water data (change numbers in filter argument to select diff months)
stn_openwater <- stn_all %>% 
filter(month(Date) %in% c(4,5,6,7,8,9,10))  

#closed water data (change numbers in filter argument to select diff months)
stn_closed <- stn_all %>%
  filter(month(Date) %in% c(11,12,1,2,3))

```

#3 Prep for IHA calculations: Get rid of years that have more than 14 consecutive days of missing data, and fill in missing values
-full dataset
```{r}

#data prep for full dataset
stn_cln <- stn_all %>%
  mutate(Date = as.Date(Date)) %>%
  filter(waterYear > 1989) %>%
  filter(waterYear != 2021 & waterYear != 2022 & waterYear != 2023)

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_full <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_dmy <- stns_ready_full %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

unique(stns_dmy$waterYear)

```

-open water data

```{r}

#data prep for open water data
stn_cln <- stn_openwater %>%
  mutate(Date = as.Date(Date)) %>%
  filter(waterYear > 1989) %>%
  filter(waterYear != 2021 & waterYear != 2022 & waterYear != 2023)

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_openwater <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_dmy_openwater <- stns_ready_openwater %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

unique(stns_dmy_openwater$waterYear)

```

-Winter/closed water data

```{r}

#data prep for winter/closed water
stn_cln <- stn_closed %>%
  mutate(Date = as.Date(Date)) %>%
  filter(waterYear > 1989) %>%
  filter(waterYear != 2021 & waterYear != 2022 & waterYear != 2023)

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_closed <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_dmy_closedwater <- stns_ready_closed %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

unique(stns_dmy_closedwater$waterYear)

```

#4 IHA calculations
-full year
-open water
-closed water

```{r}
#IHA calculations: -full year
#make a list of dfs- split by station number
lst_stns <- split(stns_dmy, stns_dmy$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
#make a "watershed" column for graphing later
IHA_fullyr <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")

###########################

#IHA calculations: -open water
#make a list of dfs- split by station number
lst_stns <- split(stns_dmy_openwater, stns_dmy_openwater$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
#make a "watershed" column for graphing later
IHA_openwater <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")

###########################

#IHA calculations: -closed water
#make a list of dfs- split by station number
lst_stns <- split(stns_dmy_closedwater, stns_dmy_closedwater$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
#make a "watershed" column for graphing later
IHA_closedwater <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")
```


```{r}


write.csv(IHA_fullyr, "PeaceRiver_PR_07HA001_IHAfullyear_1990-2020.csv")
write.csv(IHA_openwater, "PeaceRiver_PR_07HA001_IHAopenwater_1990-2020.csv")
write.csv(IHA_closedwater, "PeaceRiver_PR_07HA001_IHAclosedwater_1990-2020.csv")




```

#5 

#6 Calculate Quantiles  
```{r}

#make a list of dfs- split by station number
lst_stns <- split(stns_dmy, stns_dmy$STATION_NUMBER )
#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)
#Combine all IHA outputs into one df (unnlist the dfs)
IHA_fullyr <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")


IHA_fullyr_test <- IHA_fullyr %>%
  select('STATION_NUMBER', '1 Day Max', '7 Day Min') %>%
  filter(STATION_NUMBER == '07DA001')

#df_spl <- split(IHA_fullyr_1daymax_test, IHA_fullyr_1daymax_test$STATION_NUMBER)
View(IHA_fullyr_1daymax_test)


#####################
stn_nu <- list()
for (i in unique(IHA_fullyr_1daymax_test$STATION_NUMBER)) {
  df_subset <- IHA_fullyr_1daymax_test[IHA_fullyr_1daymax_test$STATION_NUMBER == i,]
  #st <- df_subset$STATION_NUMBER
  #stn_nu[[i]] <- i
  onedaymax <- quantile(IHA_fullyr_1daymax_test$`1 Day Max`, probs= quants)
  one_day_max <- as.data.frame(onedaymax)
  #cbind(one_day_max, STATION_NUMBER)
  #STATION_NUMBER <- unlist(stn_nu)
  df_quant <- tibble::rownames_to_column(one_day_max, "Quantiles")
  #output <- cbind(df_quant, STATION_NUMBER)
  print(df_quant)
}

##################


#one station, two values
  onedaymax <- quantile(IHA_fullyr_test$`1 Day Max`, probs= quants)
  one_day_max <- as.data.frame(onedaymax)
  sevendaymin <- quantile(IHA_fullyr_test$`7 Day Min`, probs= quants)
  sevendaymin <- as.data.frame(sevendaymin)
  sevendaymin <- tibble::rownames_to_column(sevendaymin, "Quantiles")
  STATION_NUMBER <- IHA_fullyr_test$STATION_NUMBER[0:4]
  df<- cbind(STATION_NUMBER, sevendaymin, one_day_max )
  rownames(df) <- NULL
  print(df)

  
  
calc_quantiles <- function(data){
  quants = c(0,.25,.75,1)
  onedaymax <- quantile({{data}}$`1 Day Max`, probs= quants)
  one_day_max <- as.data.frame(onedaymax)
  sevendaymin <- quantile({{data}}$`7 Day Min`, probs= quants)
  seven_day_min <- as.data.frame(sevendaymin)
  sevendaym <- tibble::rownames_to_column(seven_day_min, "Quantiles")
  STATION_NUMBER <- {{data}}$STATION_NUMBER
  df <- cbind.data.frame(STATION_NUMBER, sevendaym, one_day_max )
  rownames(df) <- NULL
  return(df)
}


IHA_fullyr_test <- IHA_fullyr %>%
  select(`7 Day Min`, `1 Day Max`,'STATION_NUMBER' )


stnsplit <- split(IHA_fullyr_test, IHA_fullyr_test$STATION_NUMBER)
lst_stns_test <- lapply(stnsplit, calc_quantiles)
bnd <- bind_rows(lst_stns_test, .id = 'STATION_NUMBER' )
```

##CALCULATE QUANTILES AND EXPORT AS CSVS
```{r}
 
calc_quantiles <- function(data){
  quants = c(0,.25,.75,1)
  onedaymax <- quantile({{data}}$`1 Day Max`, probs= quants)
  one_day_max <- as.data.frame(onedaymax)
  sevendaymin <- quantile({{data}}$`7 Day Min`, probs= quants)
  seven_day_min <- as.data.frame(sevendaymin)
  sevendaym <- tibble::rownames_to_column(seven_day_min, "Quantiles")
  STATION_NUMBER <- {{data}}$STATION_NUMBER
  df <- cbind.data.frame(STATION_NUMBER, sevendaym, one_day_max )
  rownames(df) <- NULL
  return(df)
}


#full year
#calc_quantiles(IHA_fullyr)

IHA_fullyr_71  <- IHA_fullyr %>%
  select(`7 Day Min`, `1 Day Max`,'STATION_NUMBER' )

stnsplit <- split(IHA_fullyr_71, IHA_fullyr_71$STATION_NUMBER)
lst_stns_test <- lapply(stnsplit, calc_quantiles)
df <- bind_rows(lst_stns_test, .id = 'STATION_NUMBER' )
rm(stnsplit,lst_stns_test )
#write.csv(output1, "Node_States_FullYear_Flow.csv")

quantile(IHA_fullyr$`1 Day Max`, probs = quants)

#open water
calc_quantiles(IHA_openwater)
IHA_openwater_71 <- IHA_openwater %>%
  select(`7 Day Min`, `1 Day Max`,'STATION_NUMBER' )

stnsplit <- split(IHA_openwater_71, IHA_openwater_71$STATION_NUMBER)
lst_stns_test <- lapply(stnsplit, calc_quantiles)
output2 <- bind_rows(lst_stns_test, .id = 'STATION_NUMBER')

rm(stnsplit,lst_stns_test )

quantile(IHA_openwater$`7 Day Min`, probs = quants)
#write.csv(output2, "Node_States_OpenWater_Flow.csv")


#closed water
calc_quantiles(IHA_closedwater)
IHA_closed_71 <- IHA_closedwater %>%
    select(`7 Day Min`, `1 Day Max`,'STATION_NUMBER' )

stnsplit <- split(IHA_closed_71, IHA_closed_71$STATION_NUMBER)
lst_stns_test <- lapply(stnsplit, calc_quantiles)
df <- bind_rows(lst_stns_test, .id = 'STATION_NUMBER' )

quantile(IHA_closedwater$`1 Day Max`, probs = quants)

#write.csv(output1, "Node_States_ClosedWater_Flow.csv")

```

#7 Export IHA variables to csv with a separate csv by station number

```{r}
#select variables for each dataset

IHA_fullyr_s <- IHA_fullyr %>% 
  select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")

IHA_openwater_s <- IHA_openwater  %>% 
  select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")

IHA_closedwater_s <- IHA_closedwater  %>%
select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")


#write a csv for each station number- full year, IHA
for (i in unique(IHA_fullyr_s$STATION_NUMBER)) {
  dfsubset <- IHA_fullyr_s[IHA_fullyr_s$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_fullyr_flow_IHA.csv"))
}

#write a csv for each station number- open water, IHA
for (i in unique(IHA_openwater_s$STATION_NUMBER)) {
  dfsubset <- IHA_openwater_s[IHA_openwater_s$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_openwater_flow_IHA.csv"))
}

#write a csv for each station number- closed water, IHA
for (i in unique(IHA_closedwater_s$STATION_NUMBER)) {
  dfsubset <- IHA_closedwater_s[IHA_closedwater_s$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_closed_flow_IHA.csv"))
}





```


#8 Calculate ice variables
-full dataset 
```{r}

#data prep 
#change the date format, and fill any values that are missing with preceding values 
stns_yearmonthday <- stns_ready_full %>%
  mutate(Date = as.Date(Date)) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup')

# make a list of dfs- split by station number
lst_stns <- split(stns_yearmonthday, stns_yearmonthday$STATION_NUMBER )

#apply ice variable function to all dfs in the list, by station 
stn_g1 <- lapply(lst_stns, Group_1_ice_cover)
stn_g2 <- lapply(lst_stns, Group_2_freeze_thaw)
#stn_g3 <- lapply(lst_stns, Group_3_freshet)

#Unnlist the outputs
output1 <- bind_rows(stn_g1, .id = "STATION_NUMBER") 
output2 <- bind_rows(stn_g2, .id = "STATION_NUMBER")
#output3 <- bind_rows(stn_g3, .id = "STATION_NUMBER")

#Deal with each output df having different lengths and not being able to combine them
#find max and min years for each result df
min_yro1 <- min(output1$waterYear)
max_yro1 <- max(output1$waterYear)

min_yro2 <- min(output2$waterYear)
max_yro2 <- max(output2$waterYear)

#min_yro3 <- min(output3$waterYear)
#max_yro3 <- max(output3$waterYear)

#expand data set to include all years for all output dfs
output1_ex <- output1 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro1:max_yro1) #fills NA for values for missing years 

output2_ex <- output2 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro2:max_yro2) #fills NA for values for missing years 

#output3_ex <- output3 %>%
#mutate(waterYear = as.integer(waterYear)) %>%
#complete(STATION_NUMBER, waterYear = min_yro3:max_yro3) #fills NA for values for missing years 

df_ice_full <- cbind(output1_ex, output2_ex) #output3_ex

df_ice_full_final <- df_ice_full[!duplicated(as.list(df_ice_full))] #removed duplicated columns (delete this if columns other than station number are the same)


```

-open water
```{r}

#data prep 
#change the date format, and fill any values that are missing with preceding values 
stns_yearmonthday <- stns_ready_openwater %>%
  mutate(Date = as.Date(Date)) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup')

# make a list of dfs- split by station number
lst_stns <- split(stns_yearmonthday, stns_yearmonthday$STATION_NUMBER )

#apply ice variable function to all dfs in the list, by station 
stn_g1 <- lapply(lst_stns, Group_1_ice_cover)
stn_g2 <- lapply(lst_stns, Group_2_freeze_thaw)
#stn_g3 <- lapply(lst_stns, Group_3_freshet)

#Unnlist the outputs
output1 <- bind_rows(stn_g1, .id = "STATION_NUMBER") 
output2 <- bind_rows(stn_g2, .id = "STATION_NUMBER")
#output3 <- bind_rows(stn_g3, .id = "STATION_NUMBER")

#Deal with each output df having different lengths and not being able to combine them
#find max and min years for each result df
min_yro1 <- min(output1$waterYear)
max_yro1 <- max(output1$waterYear)

min_yro2 <- min(output2$waterYear)
max_yro2 <- max(output2$waterYear)

#min_yro3 <- min(output3$waterYear)
#max_yro3 <- max(output3$waterYear)

#expand data set to include all years for all output dfs
output1_ex <- output1 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro1:max_yro1) #fills NA for values for missing years 

output2_ex <- output2 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro2:max_yro2) #fills NA for values for missing years 

#output3_ex <- output3 %>%
#mutate(waterYear = as.integer(waterYear)) %>%
#complete(STATION_NUMBER, waterYear = min_yro3:max_yro3) #fills NA for values for missing years 

df_ice_opwater <- cbind(output1_ex, output2_ex) #output3_ex

df_ice_opwater_final <- df_ice_opwater[!duplicated(as.list(df_ice_opwater))] #removed duplicated columns (delete this if columns other than station number are the same)


```

-closed water
```{r}

#data prep 
#change the date format, and fill any values that are missing with preceding values 
stns_yearmonthday <- stns_ready_closed %>%
  mutate(Date = as.Date(Date)) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup')

# make a list of dfs- split by station number
lst_stns <- split(stns_yearmonthday, stns_yearmonthday$STATION_NUMBER )

#apply ice variable function to all dfs in the list, by station 
stn_g1 <- lapply(lst_stns, Group_1_ice_cover)
stn_g2 <- lapply(lst_stns, Group_2_freeze_thaw)
#stn_g3 <- lapply(lst_stns, Group_3_freshet)

#Unnlist the outputs
output1 <- bind_rows(stn_g1, .id = "STATION_NUMBER") 
output2 <- bind_rows(stn_g2, .id = "STATION_NUMBER")
#output3 <- bind_rows(stn_g3, .id = "STATION_NUMBER")

#Deal with each output df having different lengths and not being able to combine them
#find max and min years for each result df
min_yro1 <- min(output1$waterYear)
max_yro1 <- max(output1$waterYear)

min_yro2 <- min(output2$waterYear)
max_yro2 <- max(output2$waterYear)

#min_yro3 <- min(output3$waterYear)
#max_yro3 <- max(output3$waterYear)

#expand data set to include all years for all output dfs
output1_ex <- output1 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro1:max_yro1) #fills NA for values for missing years 

output2_ex <- output2 %>%
mutate(waterYear = as.integer(waterYear)) %>%
complete(STATION_NUMBER, waterYear = min_yro2:max_yro2) #fills NA for values for missing years 

#output3_ex <- output3 %>%
#mutate(waterYear = as.integer(waterYear)) %>%
#complete(STATION_NUMBER, waterYear = min_yro3:max_yro3) #fills NA for values for missing years 

df_ice_clwater <- cbind(output1_ex, output2_ex) #output3_ex

df_ice_clwater_final <- df_ice_clwater[!duplicated(as.list(df_ice_clwater))] #removed duplicated columns (delete this if columns other than station number are the same)

```

#9 Export ice variables to csv

```{r}

#write a csv for each station number- full year, ice variables
for (i in unique(df_ice_full_final$STATION_NUMBER)) {
  dfsubset <- df_ice_full_final[df_ice_full_final$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_fullyr_icevar.csv"))
}

#write a csv for each station number- open water, ice variables
for (i in unique(df_ice_opwater_final$STATION_NUMBER)) {
  dfsubset <- df_ice_opwater_final[df_ice_opwater_final$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_openwater_icevar.csv"))
}

#write a csv for each station number- closed water, ice variables
for (i in unique(df_ice_clwater_final$STATION_NUMBER)) {
  dfsubset <- df_ice_clwater_final[df_ice_clwater_final$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_closedwater_icevar.csv"))
}


```

##WATER LEVEL DATA

#1 extract water level for stations and add column for water year and day of the year
```{r}

stns_wl <- tidyhydat::hy_daily_levels(station_number = '07DD011')

stn_all <- stns_wl %>%
group_by(STATION_NUMBER) %>%
#make a complete set of days for each year 
complete(Date = seq.Date(as.Date("1989/10/1"), as.Date("2021/09/30"), by="day"))

#add a regular year
stn_all <- stn_all %>%
    mutate(Date = ymd(Date),  Year = format_ISO8601(Date, precision = "y")) 

#add water year
stn_all <- addWaterYear(stn_all) %>%
mutate(waterYear = as.character(waterYear))

#add day of the year from "Eflows_FUNCTIONS.R" function 
stn_all_split <- split(stn_all, stn_all$STATION_NUMBER )
stn_all_split_doy <- lapply(stn_all_split, calc_day_of_wyear)
stn_all <- bind_rows(stn_all_split_doy, .id = "STATION_NUMBER")

rm(stn_all_split, stn_all_split_doy)

```

#Correct for elevation 

```{r}
#Note: WL needs datum elevation corrections but there are incongruities as to which have been done in tidyhydat already. 

#These rules were followed: 

#1) if both min and max wl values are around 200, assume corrections were done and do not correct

#2) if min < 10, and and max >200, assume no corrections were done and make all corrections

#3) if both min and max wl values are low (<200), assume the AEM 204 correction was done and add the GNSS correction to all wl values.

min(stn_all$Value, na.rm = TRUE) 
max(stn_all$Value, na.rm = TRUE) 


#WL conversion factor corrections:

#Case 1) - do nothing

#Case 2)
stn_all <- stn_all %>%
  mutate(Value = case_when(Year < 2014 ~  (Value - 204 + 204.256),
                          Year  >= 2014 ~ (Value + 204.256)))

#Case 3)  
stn_all <- stn_all %>%
  mutate(Value = Value + 207.246)
  
#outputs should be around 200
min(stn_all$Value, na.rm = TRUE)
max(stn_all$Value, na.rm = TRUE)


#LAKE CLAIRE NEAR OUTLET TO PRAIRIE RIVER 07KF002
#MAMAWI LAKE CHANNEL AT OLD DOG CAMP	07KF003
#LAKE ATHABASCA AT FORT CHIPEWYAN	 07MD001
#Slave River at Fort Fitzgerald  07NB001: (case #3)
#Peace River Below Chenal Des Quatres Fourches 07KC005 (case #2)
#Riviere Des Rochers Above Slave River 07NA001 (case #2)
#Athabasca River at Embarras Airport 07DD001 (case #1)
#Peace River at Peace Point 07KC001 (case #3)
#Lake Athabasca at Fort Chipewyan 07MD001 (case #2)
#Peace River at Peace Point 07KC001 (case #3)
#Peace River below Chenal des Quatre Fourches 07KC005 (case #2)
#ATHABASCA RIVER ABOVE JACKFISH CREEK	07DD007 (case #2)
#ATHABASCA RIVER NEAR OLD FORT	07DD011 (case #2) 


```

#Freeze up dates wL: merge ice variable df (output #2) with the corrected stn_all for WL df

```{r}

#rename so the columns have the same name
output2_rename <- output2 %>%
  rename("Date" = "Freeze_Date")

library(data.table)
setDT(output2_rename)
setDT(stn_all)

#merge and just keep the ice variable outputs with the freeze dates
merged <- merge(output2_rename, stn_all, by = "Date", all.x = TRUE)

#next calc quantiles for Values column (Wls)

quants = c(0,.25,.75,1)

quantile(merged$Value, probs = quants, na.rm =TRUE )

df <- merged %>%
na.omit("Value") 

df$waterYear.x

```

#2 export water level data to csvs
Note: these csvs are untidied flow data. No estimates, no deleting years with >14 days of missing data, etc. See below for tidied data

```{r}
#before do this make sure to get rid of years we don't want, this wasn't done above

#df1 <- stn_all %>% filter(STATION_NUMBER == "07KF002")
#df2 <- stn_all %>% filter(STATION_NUMBER == "07KF003")
#df3 <- stn_all %>% filter(STATION_NUMBER == "07MD001")

#write.csv(df1, "LakeClaire_07KF002_WL.csv")
#write.csv(df2, "MamawiLake_07KF003_WL.csv")
#write.csv(df3, "LakeAthabasca_07MD001_WL.csv")

```

#3 Make three separate datasets in prep for IHA calculations: 
-full year
-open water data (April-October) 
-winter/closed water (November-March)
```{r}

stn_all <- stn_all # full dataset

#open water data (change numbers in filter argument to select diff months)
stn_openwater <- stn_all %>% 
filter(month(Date) %in% c(4,5,6,7,8,9,10))  

#closed water data (change numbers in filter argument to select diff months)
stn_closed <- stn_all %>%
  filter(month(Date) %in% c(11,12,1,2,3))

```

#4 Prep for IHA calculations: Get rid of years that have more than 14 consecutive days of missing data, and fill in missing values
-full dataset

```{r}

#data prep for full dataset
stn_cln <- stn_all %>%
  mutate(Date = as.Date(Date)) %>%
  filter(waterYear > 1989) %>%
  filter(waterYear != 2023) 
  #filter(waterYear != 2021 & waterYear != 2022 & waterYear != 2023) 
  

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_full_wl <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_dmy_wl <- stns_ready_full_wl %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

 yrs_stn_full <- unique(stns_dmy_wl$waterYear)
 

```

-open water

```{r}

#data prep for open water dataset
stn_cln <- stn_openwater %>%
  mutate(Date = as.Date(Date)) %>%
 filter(waterYear > 1989) %>%
  filter(waterYear != 2023) 

 # filter(waterYear != 2021 & waterYear != 2022) 

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_opwater_wl <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_opwater_wl <- stns_ready_opwater_wl %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

 yrs_stn_openwater <- unique(stns_opwater_wl$waterYear)

```

-closed water

```{r}

#data prep for closed water dataset
stn_cln <- stn_closed %>%
  mutate(Date = as.Date(Date)) %>%
 filter(waterYear > 1989) %>%
    filter(waterYear != 2023) 

# filter(waterYear != 2021 & waterYear != 2022) 

#delete years that have >14 days missing data using calc_rle function from functions for this project
lst_stns_cln <- split(stn_cln, stn_cln$STATION_NUMBER )
lst_stns_cln_rle <- lapply(lst_stns_cln, calc_rle)

#unlist the station numbers
stns_ready_cl_wl <- bind_rows(lst_stns_cln_rle, .id = "STATION_NUMBER") 

#change date format and fill any values that are missing with preceding values 
stns_dmy_cl_wl <- stns_ready_cl_wl %>%
  mutate(Date = format(Date,"%d-%m-%Y")) %>%
  group_by(STATION_NUMBER) %>%
  fill(Value, .direction = 'downup') #fill missing values going down then up (ie at start of df)

yrs_stn_clwater <- unique(stns_dmy_cl_wl$waterYear)

```

#5 IHA calculations
-full year
-open water
-closed water

```{r}

#IHA calculations: -full year
#make a list of dfs- split by station number
lst_stns <- split(stns_dmy_wl, stns_dmy_wl$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
IHA_fullyr_wl <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")

###########################

#IHA calculations: -open water
#make a list of dfs- split by station number
lst_stns <- split(stns_opwater_wl, stns_opwater_wl$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
#only works for group 2 output
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
IHA_openwater_wl <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")

###########################

#IHA calculations: -closed water
#make a list of dfs- split by station number
lst_stns <- split(stns_dmy_cl_wl, stns_dmy_cl_wl$STATION_NUMBER )

#apply calc_IHA function to all dfs in the list
lst_stns_IHA <- lapply(lst_stns, calc_IHA)

#Combine all IHA outputs into one df (unnlist the dfs)
IHA_closedwater_wl <- bind_rows(lst_stns_IHA, .id = "STATION_NUMBER")

```

#Calculate Quantiles 
```{r}
#function

calc_quantiles <- function(data){
  quants = c(0,.25,.75,1)
  onedaymax <- quantile({{data}}$`1 Day Max`, probs= quants)
  one_day_max <- as.data.frame(onedaymax)
  #sevendaymin <- quantile({{data}}$`7 Day Min`, probs= quants)
  #seven_day_min <- as.data.frame(sevendaymin)
  #sevendaym <- tibble::rownames_to_column(seven_day_min, "Quantiles")
  one_day_max <- tibble::rownames_to_column(one_day_max, "Quantiles")
  STATION_NUMBER <- {{data}}$STATION_NUMBER
  return(one_day_max)
  #df <- cbind.data.frame(STATION_NUMBER, one_day_max )
  #rownames(df) <- NULL
  #return(df)
}


```

-full year
```{r}

#full year
df <- IHA_fullyr_wl %>%
    select(`1 Day Max`,'STATION_NUMBER' )

stnsplit <- split(df, df$STATION_NUMBER)
lst_stns_test <- lapply(stnsplit, calc_quantiles)
output4 <- bind_rows(lst_stns_test, .id = 'STATION_NUMBER' )

#write.csv(output4, "FullYear_WL07kc005.csv")
```

```{r}
#open water

df <- IHA_openwater_wl %>%
    select(`1 Day Max`,'STATION_NUMBER')

stnsplit <- split(df, df$STATION_NUMBER)
lst_stns <- lapply(stnsplit, calc_quantiles)
output5 <- bind_rows(lst_stns, .id = 'STATION_NUMBER' )


#write.csv(output5, "OpenWater_WL07kc005.csv")

```

```{r}
#closed water 
df <- IHA_closedwater_wl %>%
    select(`1 Day Max`,'STATION_NUMBER')

stnsplit <- split(df, df$STATION_NUMBER)
lst_stns <- lapply(stnsplit, calc_quantiles)
output6 <- bind_rows(lst_stns, .id = 'STATION_NUMBER' )


#write.csv(output6, "ClosedWater_WL07kc005.csv")


```

```{r}


write.csv(stn_cln, "Athabasca_River_OldFort_07DD011_OpenWater_WL_DatumCorrected_1990-2022.csv" )

write.csv(stn_all,  'PR_QuatresFourches_07KC005_Raw_WL_Data_DatumCorrected_1990-2020.csv')
write.csv(IHA_fullyr_wl,"PR_QuatresFourches_07KC005_fullyear_WL_DatumCorrected_1990-2020.csv")
write.csv(IHA_closedwater_wl,'PR_QuatresFourches_07KC005_closedwater_WL_DatumCorrected_1990-2020.csv')
write.csv(IHA_openwater_wl,'PR_QuatresFourches_07KC005_openwater_WL_DatumCorrected_1990-2020.csv')

```


#5 Export as csv for each station
```{r}

#select variables for each dataset

IHA_fullyrwl_S <- IHA_fullyr_wl %>% 
  select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")

IHA_openwater_wl_s <- IHA_openwater_wl  %>% 
  select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")

IHA_closedwater_wl_s <- IHA_closedwater_wl  %>%
select("STATION_NUMBER", "Year", "1 Day Max", "30 Day Min", "30 Day Max", "Low pulse length", "High pulse length", "Reversals")


#write a csv for each station number- full year, IHA
for (i in unique(IHA_fullyrwl_S$STATION_NUMBER)) {
  dfsubset <- IHA_fullyrwl_S[IHA_fullyrwl_S$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_fullyr_wl_IHA.csv"))
}

#write a csv for each station number- open water, IHA
for (i in unique(IHA_openwater_wl_s$STATION_NUMBER)) {
  dfsubset <- IHA_openwater_wl_s[IHA_openwater_wl_s$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_openwater_wl_IHA.csv"))
}

#write a csv for each station number- closed water, IHA
for (i in unique(IHA_closedwater_wl_s$STATION_NUMBER)) {
  dfsubset <- IHA_closedwater_wl_s[IHA_closedwater_wl_s$STATION_NUMBER == i,]
  write.csv(dfsubset, paste0(i, "_closed_wl_IHA.csv"))
}

```