create_hydrometric_map.R

# Packages ----------------------------------------------------------------

library(sf)
library(dplyr)
library(tidyhydat)
library(readr)
library(ggplot2)
library(here)
library(raster)


# Read in data ------------------------------------------------------------

## Load spatial layers
data_path <- here::here("data")

canada_border <- raster::getData('GADM', country="CAN", path = data_path, level=1) 

proj_var <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs"


canada_border <- st_as_sf(canada_border) %>%
  st_transform(proj_var) 

## Load realtime stations and convert to spatial object
stns_sf <-  st_as_sf(
  suppressMessages(allstations),
  coords = c("LONGITUDE", "LATITUDE"),
  crs = 4326,
  agr = "constant"
) %>%
  filter(HYD_STATUS == "ACTIVE") %>%
  st_transform(proj_var)



# Manipulations -----------------------------------------------------------

dput(sort(unique(allstations$PROV_TERR_STATE_LOC)))

cw_table <- tibble(PROV_TERR_STATE_LOC = c("AB", "BC",  "MB", "NB", "NL", "NT",
                        "NS", "NU", "ON", "PE", "QC", "SK", "YT"),
NAME_1 = c("Alberta", "British Columbia", "Manitoba", "New Brunswick", 
  "Newfoundland and Labrador", "Northwest Territories", "Nova Scotia", 
  "Nunavut", "Ontario", "Prince Edward Island", "Québec", "Saskatchewan", 
  "Yukon"))


num_stations <- allstations %>%
  filter(HYD_STATUS == "ACTIVE") %>%
  group_by(PROV_TERR_STATE_LOC) %>%
  count() %>%
  right_join(cw_table) %>%
  ungroup() 

canada_border <- canada_border %>%
  left_join(num_stations) %>%
  mutate(area = st_area(.)/1000) %>%
  mutate(num_area = as.numeric(n/area))
  


# Plotting ----------------------------------------------------------------


ggplot() +
  geom_sf(data = canada_border, aes(fill = num_area), alpha = 0.5) +
  geom_sf(data = stns_sf, colour = "black", size = 0.8) +
  scale_fill_viridis_c(name = "Number of stations\nper size of province") +
  theme_void()