Using hydrotoolbox with Canadian data - Part 1
In hydrotoolbox: Hydrological Tools for Handling Hydro-Meteorological Data Records

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(tidyhydat)
library(dplyr)
library(hydrotoolbox)

The 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:

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.

# 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.

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

hm_plot(obj = station_data,
        slot_name = "qd",
        col_name = list("Value"))

The daily flows can be aggregated to monthly or yearly values

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")
          )

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:

# 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

# 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.