In bcgov/hydrolook: Hydrometric Network Diagnostic

---

title: "r tidyhydat::hy_stations(station_number = params$stns)$STATION_NAME - r params$stns"

library(hydrolook)
library(tidyhydat)
library(tidyverse)
library(lubridate)
library(kableExtra)
library(knitr)
library(sf)
library(bcmaps)
library(cowplot)
library(zoo)
library(ggrepel)

# Copyright 2017 Province of British Columbia
# 
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
# 
# http://www.apache.org/licenses/LICENSE-2.0
# 
# Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and limitations under the License.


## Knitr options dependent on type of document
if(params$table_format == "latex"){
  opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, 
                      fig.width = 8, fig.height = 8, fig.path = file.path("report/station_report/"))
}

if(params$table_format == "html"){
  opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, 
                      fig.width = 8, fig.height = 8)
}

options(knitr.table.format = params$table_format)

## For testing
#params <- tibble(stns = "08MF005")
## A station not in hydat
#params <- tibble(stns = "07FC001")

## Create station header
stn_header <- allstations %>%
  filter(STATION_NUMBER == params$stns)

Information

This report, generated on

cat(format(Sys.Date(), '%d %B, %Y'))

provides basic information on a hydrometric station.

rl_data <- realtime_dd(station_number = params$stns)

rl_data_day <- rl_data %>%
  realtime_daily_mean() %>%
  filter(Parameter == "Flow")

HIST_FLOW <- hy_daily_flows(station_number = params$stns)

Daily Summary

range <- hy_stn_data_range(station_number = params$stns) %>%
  filter(DATA_TYPE == "Q") %>%
  select(YEAR_FROM, YEAR_TO, RECORD_LENGTH)

tbl_flow <- HIST_FLOW %>%
  mutate(Q7_day = rollapply(Value, 7, mean, na.rm = TRUE, partial = TRUE, fill = NA, align = "right")) %>% 
  filter(month(Date) %in% month(Sys.Date()-1)) %>%
  filter(day(Date) %in% day(Sys.Date()-1)) %>%
  bind_rows(rl_data_day %>%
              ungroup() %>%
              mutate(Q7_day = rollmean(Value, 7, fill = NA, align = "right")) %>%
              filter(Date %in% (Sys.Date()-1))
              ) %>% 
  filter(!is.na(Value)) %>%
  mutate(Rank = min_rank(Value)) %>%
  mutate(Rank7day = min_rank(Q7_day)) %>%
  filter(year(Date) == year(Sys.Date())) %>%
  select(Date, Value, Rank, Q7_day, Rank7day, Symbol)

col_vector <- c("Date","Daily Discharge","Daily Rank", "7 Day Discharge", "7 Day Rank","Symbol")

colnames(tbl_flow) <- col_vector


if(nrow(tbl_flow) == 0){
  tbl_flow <- as.data.frame(t(matrix(rep(NA, length(col_vector)), byrow = FALSE)))
  colnames(tbl_flow) <- col_vector
}

kable(cbind(tbl_flow,range), format = "latex", booktabs = T) %>%
  kable_styling(font_size = 6, latex_options = c("HOLD_position")) %>%
  add_footnote(c("All discharge presented in m^3/s","7 day rolling average includes previous 7 days"), notation = "alphabet")

Plot Summaries

stns_sf <- hy_stations(station_number = params$stns) %>%
  st_as_sf(., coords = c("LONGITUDE", "LATITUDE"), 
           crs = 4326, 
           agr = "constant") %>%
  transform_bc_albers() %>%
  mutate(
    CENTROID = map(geometry, st_centroid),
    COORDS = map(CENTROID, st_coordinates),
    COORDS_X = map_dbl(COORDS, 1),
    COORDS_Y = map_dbl(COORDS, 2)
  ) %>%
  as_tibble() %>%
  st_as_sf()




## Drainages for the stations
plt_wsc <- wsc_drainages() %>%
  st_join(stns_sf) %>%
  filter(!is.na(STATION_NUMBER))

##Rivers
plt_5m <- watercourses_5M() %>%
  st_join(plt_wsc) %>%
  filter(!is.na(STATION_NUMBER))

## Larger watershed
plt_nr <- nr_regions() %>%
  st_join(plt_wsc) %>%
  filter(!is.na(MAJOR_DRAINAGE_AREA_NAME)) %>%
  mutate(REG_NAME = gsub("Natural Resource Region","", REGION_NAME)) %>%
  st_intersection(bc_bound()) 

## Find the biggest city then bind it to the station sf object
biggest_city_stn <- bc_cities() %>%
  st_join(plt_nr %>%
            group_by(OCEAN_DRAINAGE_AREA_NAME) %>%
            summarise(), left = FALSE) %>%
  filter(POP_2000 == max(POP_2000)) %>%
  rename(geometry = SHAPE) %>%
  mutate(
    CENTROID = map(geometry, st_centroid),
    COORDS = map(CENTROID, st_coordinates),
    COORDS_X = map_dbl(COORDS, 1),
    COORDS_Y = map_dbl(COORDS, 2)
  ) %>%
  as_tibble() %>%
  st_sf() %>%
  select(NAME, COORDS_X, COORDS_Y, geometry) %>%
  rbind(stns_sf %>%
          select(NAME = STATION_NAME, COORDS_X, COORDS_Y, geometry)
        )



stn_map <- ggplot() +
  geom_sf(data = plt_wsc, aes(fill = MAJOR_DRAINAGE_AREA_NAME), alpha = 0.8) +
  geom_sf(data = plt_nr, aes(fill = REGION_NAME), alpha = 0.2) +
  geom_sf(data = plt_5m, aes(colour = name_en), size = 2) +
  geom_sf(data = biggest_city_stn, size = 2, aes(shape = NAME)) +
  geom_text_repel(biggest_city_stn, mapping = aes(COORDS_X, COORDS_Y, label = NAME), 
                           size = 3, min.segment.length = 0) +
  scale_colour_viridis_d(name = "River Feature") +
  guides(fill = FALSE, shape = FALSE) +
  labs(x = "", y = "") +
  theme_minimal() +
  theme(legend.position = "bottom")

rl_plt <- rl_data_day %>%
  ggplot(aes(x = Date, y = Value)) +
  geom_path() +
  geom_point() +
  theme_minimal() +
  labs(title = "Last Thirty Days", y = "Discharge (m^3/s)") 

## A dataframe 
pct_flow <- HIST_FLOW %>%
  mutate(year_day = yday(Date)) %>%
  filter(year_day %in% yday(seq.Date(from = (Sys.Date()-30), 
                                     to = Sys.Date(), by = "day"))) %>%
  group_by(year_day, STATION_NUMBER) %>%
  mutate(prctile = ecdf(Value)(Value)) %>%
  mutate(Date_no_year = dmy(paste0(day(Date),"-",month(Date),"-",year(Sys.Date())))) %>%
  ungroup()


hist_plt <- ggplot(pct_flow, aes(x = Date_no_year, y = Value)) +
  geom_point(aes(colour = prctile)) +
  geom_line(data = rl_data_day, aes(x = Date), colour = "red") +
  geom_point(data = rl_data_day, aes(x = Date, y = Value, 
                                     shape = factor(year(Date))), colour = "red") +
  scale_colour_viridis_c(name = "Percentile \n", option = "viridis") +
  scale_shape_discrete(name = "Year") +
  theme_minimal() +
  labs(title = "Historical flow relative to current year",
       subtitle = "Current year flows are displayed in red",
       x = "Date", y = "Discharge (m^3/s)")

#plot_grid(stn_map, rl_plt, hist_plt, hist_dist_plt)

toprow <- plot_grid(stn_map, rl_plt)

plot_grid(toprow, hist_plt, ncol = 1)

bcgov/hydrolook documentation built on May 2, 2021, 11:23 p.m.