Rainfall_runoff_data_prep.R

## This script prepares inputs for rainfall-runoff modeling

# ===========================  Setup Script ====================================
# if not already run...
source('setup.R')

#=========================  Import necessary datasets ==========================

park <- "BRCA"

park_boundary <- getParkBoundary(park)
  
# Buffer area around the park ~ 100 km away:
aoi <- sf::st_buffer(park_boundary, dist = 0.3)

# Get NWIS sites in aoi
nwis <- listNWIS(aoi = aoi, dist = 0)
  
# Identify gages representative of natural conditions:
ref_gages <- nwis %>% 
  dplyr::filter(code == "00060") %>%
  dplyr::pull(site_no) %>%
    get_gagesII(id =.) %>%
    dplyr::filter(class == "Ref")

# Now join to get df of select gages
nwis_stream <- nwis %>%
    dplyr::filter(data_type_cd %in% c("dv","cv"),
                  code == "00060",
                  year(end_date) >= 1980, # 1980 onwards
                  site_no %in% ref_gages$staid) %>%
    dplyr::left_join(st_drop_geometry(ref_gages), 
                     by = c("site_no" ="staid"))
  
#Download data from reference stream sites
nwis_stream_daily <- 
    dataRetrieval::readNWISdv(siteNumbers = nwis_stream$site_no,
                              parameterCd = c('00060','00065')) %>%
    dplyr::mutate(flow_cfs = as.numeric(X_00060_00003),
                  Date = as.Date(Date) ) %>%
    dplyr::select(c(Date, site_no, flow_cfs)) %>% 
    dplyr::filter(year(Date) >= 1980) %>%
    pivot_wider(names_from = site_no, values_from = flow_cfs)
    
# Convert to monthly also
#nwis_stream_monthly <- nwis_stream_daily %>%
#    tidyr::pivot_longer(cols = -Date, 
#                        names_to = "site_no", 
#                        values_to = "flow_cfs") %>%
#    dplyr::mutate(y = year(Date), 
#                  m = month(Date)) %>%
#    dplyr::group_by(y, m, site_no) %>%
#    dplyr::summarize(mean_discharge = mean(flow_cfs, na.rm. = TRUE),
#                     .groups = "keep") %>%
#    dplyr::ungroup() %>%
#    dplyr::mutate(ym = lubridate::make_date(year = y, month = m, day = 1)) %>% 
#    dplyr::select(ym, site_no, mean_discharge) %>%
#    tidyr::pivot_wider(names_from = site_no, values_from = mean_discharge)
 
# Get watersheds associated with streamflow locations
nldi_watershed <- nwis_stream$site_no %>%
    map_dfr(~ nldi_finder(site_no = .) ) %>%
    dplyr::mutate(data = map(site_no, ~nldi_meta(site_no = .) ) ) %>%
    unnest(cols = data) %>%
    left_join(st_drop_geometry(nwis_stream), by = "site_no")

# Get flowlines associated with streamflow locations
nldi_flowlines <- mapNHDPlusHR(aoi = dplyr::summarize(nldi_watershed)) %>% 
    dplyr::summarize()
  
#Get Gw sites!
nwis_groundwater <- nwis %>%
    # locations with more than one days' worth of data:
    dplyr::filter(begin_date != end_date,
                  year(end_date) >= 2000,
                  n_obs > 10,
                  # groundwater sites only:
                  site_type_cd == "GW",
                  data_type_cd == "gw")
  
  
  # pull those sites' level data
nwis_groundwater_levels <- 
    dataRetrieval::readNWISgwl(nwis_groundwater$site_no) %>%
    dplyr::filter(parameter_cd == 72019,
                  year(lev_dt) >= 1980) %>%
    dplyr::mutate(ym = lubridate::ym(substr(lev_dt, 1, 7))) %>% 
    dplyr::group_by(ym, site_no) %>%
    dplyr::summarize(mean_lev_va = mean(lev_va, na.rm. = TRUE),
                     .groups = "keep") %>%
    dplyr::select(ym, site_no, mean_lev_va) %>%
    tidyr::pivot_wider(names_from = site_no, 
                      values_from = mean_lev_va) 

# import well data csv  
well_data <- read_csv('data/park/BRCA/manual/BRCA_Well_Data.csv', na = c("NaN", "NA", "")) %>%
  janitor::clean_names() %>%
  dplyr::mutate(static_in = as.numeric(static_in),
                meter_gpm = as.numeric(meter_gpm),
                total = as.numeric(total),
                date = mdy(date))

well_monthly <- well_data %>%
  dplyr::mutate(ym = ym(substr(date, 1, 7))) %>%
  dplyr::select(well,date,ym,static_in, meter_gpm, total) %>%
  group_by(well, ym) %>% 
  mutate(day_count = ifelse(is.na(static_in), 0, 1)) %>%
  dplyr::summarize(level_observations = sum(day_count, na.rm = TRUE),
                   level_mean = mean(static_in, na.rm = TRUE),
                   level_median = median(static_in, na.rm = TRUE),
                   level_sd = sd(static_in, na.rm = TRUE),
                   level_spread = max(static_in, na.rm = TRUE) - 
                     min(static_in, na.rm = TRUE),
                   meter_total = max(meter_gpm),
                   pump_total = max(total),
                   .groups = "keep")

# The state of Utah also tracks monthly water use of that system:
water_supply_id <- getWaterSuppliersUtah(aoi = park_boundary) %>%
  filter(grepl("National Park", WRNAME, ignore.case=TRUE)) %>%
  .$WRID

# Well 1
water_use_1 <- getWaterUseUtah(WRID = water_supply_id)[[1]] %>%
  slice(1:39) %>%
  pivot_longer(-c("Year", "Method of Measurement"), 
               names_to = "month", 
               values_to = "use_acre_feet") %>%
  mutate(ym = ym(paste0(Year, "-", month))) %>%
  filter(month != "Annual inAcre Feet") %>%
  dplyr::select(ym, use_acre_feet) %>%
  mutate(well = "Well 1")

# Well 2
water_use_2 <- getWaterUseUtah(WRID = water_supply_id)[[1]] %>%
  slice(51:nrow(.)) %>%
  dplyr::filter(!is.na(as.numeric(Year))) %>%
  pivot_longer(-c("Year", "Method of Measurement"), 
               names_to = "month", 
               values_to = "use_acre_feet") %>%
  mutate(ym = ym(paste0(Year, "-", month))) %>%
  filter(month != "Annual inAcre Feet") %>%
  dplyr::select(ym, use_acre_feet) %>%
  mutate(well = "Well 2")

# join water use data
both_wells <- water_use_1 %>%
  dplyr::bind_rows(water_use_2) %>%
  group_by(ym) %>%
  dplyr::summarize(use_acre_feet = sum(as.numeric(use_acre_feet), na.rm = TRUE))

# Here we combine the average monthly static water levels with the monthly water use:
well_munge <- well_monthly %>%
  left_join(both_wells, by = c("ym"))

# NPS tracks monthly total park visitors. Here we pull that information in for the park:
visitors <- 
  getUnitVisitation(units = "BRCA", startYear = 2000, endYear = 2023) %>%
  mutate(ym = ym(paste0(Year, "-", Month))) %>%
  dplyr::select(ym, RecreationVisitors)

well_data <- list("water_supply_id" = water_supply_id,
                  "well_munge" = well_munge,
                  "visitors" = visitors)

rm("visitors","well_munge","both_wells","water_use_1",
   "water_use_2","water_supply_id","well_data","well_monthly")


# ======================= Begin Analysis =======================================


# define conversion terms
km2_to_ft2 <-10763910.41671  # 1km2 = 10763910.41671 ft2
mm_to_in <- 0.0393701        # 1mm = .0393701 in
ft_to_in <- 12               # 1ft = 12 in 
day_to_s <- 86400            # 1day = 86400 seconds



# Select either watershed 1 (mammoth creek) or 2 (sevier)
i <- 2

site_no <- nldi_watershed$site_no[i]
site_name <- nldi_watershed$site_name[i]
watershed_area <- nldi_watershed$drain_sqkm[i] * km2_to_ft2
ws_centroid <-nldi_watershed[i,] %>% 
              sf::st_transform(4326) %>% 
              sf::st_centroid()


# NWIS discharge reported in CFS. Convert to CFD then divide by total drainage 
# area to get area averaged discharge in units in per day

discharge_daily <- nwis_stream_daily %>% 
  dplyr::select(c(Date, !!sym(site_no))) %>%
  dplyr::mutate(discharge_in_d = (!!sym(site_no) * ft_to_in * day_to_s) /
                  (watershed_area)) %>%
  dplyr::filter(!is.na(discharge_in_d))


#discharge_monthly <- nwis_stream_monthly %>% 
#  dplyr::select(c(ym, !!sym(site_no))) %>%
#  dplyr::mutate(discharge_in_d = (!!sym(site_no) * ft_to_in * day_to_s) /
#                  (watershed_area)) %>%
#  dplyr::filter(!is.na(discharge_in_d))


# Import WBM in mm/d. Convert to inch / day

# First, load daily and monthly data for park 
# Note, this was previously downloaded using getHistoricWBMGridMET() using
# an AOI that is much greater than what we need here. So, we're pulling only
# the points we want. This step processes a lot of data and takes a couple
# of minutes.

# accumswe represents the cumulative value for swe at any given time. To get
# melt, get the difference in SWE between timesteps. Any negative changes in 
# SWE represent melt. If melt + accumulation happens, this will be lost.
# Alternatively, just use SWE for precip


wbm_vars <- 
  list.files("data/misc/wbm_gridmet_hist_daily/temp/", 
             pattern = "runoff|rain|PET|accumswe",
             full.names = TRUE) %>%
  purrr::map(~data.table::fread(.)) %>% 
  dplyr::bind_rows()

wbm_xy_daily <- wbm_vars %>%
              dplyr::mutate(val = val/10) %>%
              tidyr::pivot_wider(names_from = "var", 
                     values_from = "val") %>%
                raster_puller(data = ., 
                              aoi = NULL, 
                              point = ws_centroid) %>%
                dplyr::mutate(runoff = mm_to_in*runoff, # Reported * 10 (dbl check)
                              date = as.Date(date),
                              melt = -1* (accumswe - lag(accumswe, 
                                                    default = accumswe[1])),
                              melt = ifelse(melt < 0, 0, melt),
                              effR = melt + rain) %>%
                dplyr::select(-ym)


#wbm_xy_monthly <- wbm_xy_daily %>%
#                  dplyr::mutate(y = year(date), m = month(date)) %>%
#                  dplyr::group_by(y,m) %>%
#                  dplyr::summarize(runoff = sum(runoff, na.rm. = TRUE),
#                                   melt = sum(melt, na.rm. = TRUE),
#                                   rain = sum(rain, na.rm. = TRUE),
#                                   PET = sum(PET, na.rm. = TRUE),
#                                   accumswe = sum(accumswe, na.rm. = TRUE),
#                                   effR =  sum(effR, na.rm. = TRUE),
#                                  .groups = "keep") %>%
#                  dplyr::mutate(ym = lubridate::make_date(year = y, month = m, day = 1)) %>%
#                  ungroup() %>% dplyr::select(-c(y,m))

  
# Import wbm for full watershed as mean

wbm_ws_daily <- wbm_vars %>%
                  dplyr::mutate(val = val/10) %>% # Reported * 10 (dbl check)
                  tidyr::pivot_wider(names_from = "var", 
                                    values_from = "val") %>%
                  raster_puller(data = ., 
                                aoi = nldi_watershed[i,1], 
                                point = NULL) %>%
                  dplyr::mutate(runoff = mm_to_in*runoff, 
                                  date = as.Date(date),
                                  melt = -1* (accumswe - lag(accumswe, 
                                                      default = accumswe[1])),
                                  melt = ifelse(melt < 0, 0, melt),
                                  effR = melt + rain) %>%
                  dplyr::group_by(date) %>%
                  dplyr::summarize(runoff = mean(runoff, na.rm. = TRUE),
                                   accumswe = mean(accumswe, na.rm. = TRUE),
                                   rain = mean(rain, na.rm. = TRUE),
                                   PET = mean(PET, na.rm. = TRUE),
                                   melt = mean(melt, na.rm. = TRUE),
                                   effR = mean(effR, na.rm. = TRUE),
                                   .groups = "keep") 
 
#wbm_ws_monthly <- wbm_ws_daily %>%
#                  dplyr::mutate(y = year(date), m = month(date)) %>%
#                  dplyr::group_by(y,m) %>%
#                  dplyr::summarize(runoff = sum(runoff, na.rm. = TRUE),
#                                   accumswe = sum(accumswe, na.rm. = TRUE),
#                                   rain = sum(rain, na.rm. = TRUE),
#                                   PET = sum(PET, na.rm. = TRUE),
#                                   melt = sum(melt, na.rm. = TRUE),
#                                   effR = sum(effR, na.rm. = TRUE),
#                                   .groups = "keep") %>%
#                  dplyr::mutate(ym = lubridate::make_date(year = y, month = m, day = 1)) %>%
#                  dplyr::ungroup() %>% 
#                    dplyr::select(-c(y,m))


tp <- "daily" # or "daily"
span <- "watershed" # or "centroid"


if (tp == "daily" & span == "centroid") { 
    wbm_nwis <- wbm_xy_daily
  } else if (tp == "daily" & span == "watershed") { 
    wbm_nwis <- wbm_ws_daily
  } else if (tp == "monthly" & span == "centroid") { 
    wbm_nwis <- wbm_xy_monthly
  } else if (tp == "monthly" & span == "watershed") {
    wbm_nwis <- wbm_ws_monthly
}


if (tp == "daily") {
  wbm_nwis <- wbm_nwis %>% 
    dplyr::left_join(., discharge_daily, by = c("date" = "Date")) %>%
    dplyr::filter(year(date) < max(year(discharge_daily$Date)))
  } else if (tp == "monthly") {
    wbm_nwis <- wbm_nwis %>%
  dplyr::left_join(.,discharge_monthly, by = c("ym" = "ym")) %>% 
    dplyr::filter(year(ym) < max(year(discharge_monthly$ym))) 
  }

# uncomment to save data
#saveRDS(wbm_nwis,"data/park/BRCA/Rainfall_runoff_modeling/mammoth_creek_daily_wbm_discharge.RDS")