Find all current hydrometric stations then subset by NRD

Station_selection.R

library(tidyverse)
library(sf)
library(tidyhydat)
library(bcmaps)

## Say get Thompson District
nr_district_to <- st_as_sf(nr_districts) %>%
  filter(REG_ORG_UN == "Thompson-Okanagan Natural Resource Region")

## Combine HYDAT stations and realtime stations
## Eliminate duplicated stations then convert to sf object
stns <- STATIONS(PROV_TERR_STATE_LOC = "BC") %>%
  filter(HYD_STATUS == "ACTIVE") %>%
  bind_rows(., realtime_network_meta(PROV_TERR_STATE_LOC = "BC")) %>%
  distinct(STATION_NUMBER, .keep_all = TRUE) %>%
  select(STATION_NUMBER, STATION_NAME, LATITUDE, LONGITUDE) %>%
  st_as_sf(coords = c("LONGITUDE","LATITUDE"), 
           crs = 4326) %>%
  transform_bc_albers()


### Find all stations in Thompson Okanagan Region
## Export it to .csv
stns %>%
  st_join(nr_district_to, left = FALSE) %>%
  select(-F_NAME, -SHAPE, -OBJECTID,-F_CODE) %>%
  st_set_geometry(NULL) %>%
  write_csv("thompson-okanagan-stations.csv")