---
title: "Using hydrotoolbox with Canadian data - Part 1"
author: "Kevin Shook and Ezequiel Toum"
date: "`r Sys.Date()`"
output: 
 rmarkdown::html_vignette:
   toc: true
vignette: >
  %\VignetteIndexEntry{Using hydrotoolbox with Canadian data - Part 1}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r libraries, eval = FALSE}
library(tidyhydat)
library(dplyr)
library(hydrotoolbox)
```

The package [**tidyhydat**](https://CRAN.R-project.org/package=tidyhydat)
provides a very easy way to access Canadian hydrometric data. As both **tidyhydat** and `hydrotoolbox` methods work with
the **tidyverse**, they work very well with each other.

## Installing tidyhydat and HYDAT

Once you have installed **tidyhydat**, you will need to
download the HYDAT sqlite database of Canadian hydrometric
data. The database is downloaded as soon as you first attempt
to query it. It can also be downloaded manually:

```{r download, eval = FALSE}
download_hydat()
```

Note that the database is updated from time to time, and you
will be prompted to download the most recent version.

## Combining tidyhydat with hydrotoolbox

Once the database is downloaded, you can query it for data
by station ID number. In this example
we add the daily flows from **tidyhydat** to the `qd` slot.

```{r station, eval = FALSE}
# Smith Creek near Marchwell, in Saskatchewan
station_number <- "05ME007"    

daily_vals <- 
  hy_daily_flows(station_number)
  

station_data <- 
  hm_create() %>%
  hm_set(qd = daily_vals[ , c("Date", "Value")])
```

You can also add the station metadata to the station object. In Canada,
many basin have variable contributing fractions, caused by the
existence of depressional storage. In these basins, the effective
area is that portion of the basin which contributes flow at least
one year in two. Some Canadian basins are designated as being
part of the Reference Hydrometric Basin Network, which consists
of stations having long records without much human interference. 

```{r metadata, eval = FALSE}
meta_data <- hy_stations(station_number)

station_data <- hm_set(station_data, 
                       id = meta_data$STATION_NUMBER[1],
                       station = meta_data$STATION_NAME[1], 
                       country = "Canada",
                       lat = meta_data$LATITUDE[1],
                       long = meta_data$LONGITUDE[1],
                       province = meta_data$PROV_TERR_STATE_LOC[1],
                       active = (meta_data$HYD_STATUS[1] == "ACTIVE"),
                       basin_area = meta_data$DRAINAGE_AREA_GROSS[1],
                       basin_eff = meta_data$DRAINAGE_AREA_EFFECT[1],
                       other_1 = paste("RHBN:", meta_data$RHBN[1]))
```

Having read in the data, we can now plot the flows

```{r daily_plot, fig.width = 6, fig.height = 4, eval = FALSE}
hm_plot(obj = station_data,
        slot_name = "qd",
        col_name = list("Value"))
```

The daily flows can be aggregated to monthly or yearly values

```{r monthly_plot, fig.width = 6, fig.height = 4, eval = FALSE}
monthly <- hm_agg(station_data, 
                  slot_name = "qd", 
                  fun = "mean", 
                  period = "monthly",
                  col_name = "Value")

monthly %>%
  hm_plot(slot_name = "qd",
          col_name = list("Value_mean")
          )
```

```{r annual_plot, fig.width = 6, fig.height = 4, warning = FALSE, eval = FALSE}

yearly <- hm_agg(station_data, 
                 slot_name = "qd", 
                 fun = "mean", 
                 period = "annually", 
                 col_name = "Value")

yearly %>% 
  hm_plot(slot_name = "qd", 
          col_name = list("Value_mean"))
```

## Using an own made function to build the object

Although this version of the `hydrotoolbox` package does
not offer the ability to build
automatically a hydro-meteorological station from HYDAT database, you can save a lot of 
time by recycling the following function:

```{r my_fun, eval = FALSE}
# before running this function, the packages
# tidyhydat and hydrotoolbox should be attached

# station_number: character with station ID
build_hydat <- function(station_number){
  # get daily flow series
  q <-
    hy_daily_flows(station_number) %>%
    as.data.frame()
  
  # get station metadata
  meta_table <- 
    hy_stations(station_number) %>%
    as.data.frame()
  
  # set values
  out_hm <- 
    hm_create() %>%
    hm_set(qd = q[ , c("Date", "Value")], 
           id = meta_table$STATION_NUMBER[1],
           station = meta_table$STATION_NAME[1], 
           country = "Canada",
           lat = meta_table$LATITUDE[1],
           long = meta_table$LONGITUDE[1],
           province = meta_table$PROV_TERR_STATE_LOC[1],
           active = (meta_table$HYD_STATUS[1] == "ACTIVE"),
           basin_area = meta_table$DRAINAGE_AREA_GROSS[1],
           basin_eff = meta_table$DRAINAGE_AREA_EFFECT[1],
           other_1 = paste("RHBN:", meta_table$RHBN[1]))
}
```

once the function is loaded in the *Global Environment*,
we set up the station

```{r fun_applied, eval = FALSE}
# we construct the Smith Creek near Marchwell 
# but in a single code line

smith_creek <- build_hydat(station_number = "05ME007")

smith_creek %>% 
  hm_show(slot_name = c("id", "station", "country", "qd"))
```

Since the builder function is the only one that differs
from what was developed for SNIH data, we recommend
(re)visiting this vignette (`vignette ('snih_arg', package = 'hydrotoolbox')`) to explore some of the available methods.