utils.R

needs_layer <- function(db, layer) {
  if(file.exists(db)) {
    layers <- st_layers(db)
    if(layer %in% layers$name)
      return(FALSE)
  }
  TRUE
}

layer_exists <- function(db, layer) {
  if(file.exists(db)) {
    layers <- st_layers(db)
    if(layer %in% layers$name)
      return(TRUE)
  }
  FALSE
}

# Subset NHD RDS by hydrologic region
VPU_Subset <- function(nhdPath, VPU){
  # Read in COMIDs of outlets
  RegOutlets <- jsonlite::read_json(file.path("cache", "RegOutlets.json"))
  print (paste0("subsetting for VPU ", VPU))
  
  if (!VPU %in% names(RegOutlets)){
    
    # For regions that terminate to the NHDPlus domain
    # 1,2,3,4,8,9,12,15,17,18
    # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
    nhd <- readRDS(nhdPath) %>% 
      filter(grepl(paste0("^", VPU, ".*"), .data$VPUID)) %>% filter(FTYPE != "Coastline")
    
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID)
    
  } else if (VPU %in% c("05", "06", "07", "10U", "14")){
    
    # Subset VPUs that are single outlet
    outlet <- RegOutlets[VPU]
    #print (outlet)
    
    nhd <- readRDS(nhdPath)
    
    #nhd_US <- get_UT(nhd, outlet)
    
    nhd <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, outlet)))
     
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID, VPUID == VPU)
    
  } else if (VPU == "10L"){

    # Subset by flowlines in Region 10
    nhd <- readRDS(nhdPath) %>% 
      filter(grepl(paste0("^",substr(VPU,1,2),".*"), .data$VPUID)) # take out this VPUID
    
    # R10U upstream
    nhd_US_10U <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, RegOutlets["10U"])))
    
    # R10L upstream
    nhd_US_10L <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, RegOutlets["10L"])))
    
    nhd_10L <- nhd_US_10L %>% filter(!COMID %in% nhd_US_10U$COMID)
    
    keep <- prepare_nhdplus(nhd_10L, 0, 0, 0, FALSE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID)
    
  } else {
    
    outlets <- unlist(RegOutlets[VPU])
    
    nhd <- readRDS(nhdPath) 
    
    #outlets <- c(14320629, 941140164, 15334434)
    
    for (out in outlets){
      #print (out)
      nhd_US <- get_UT(nhd, out)
      keep <- prepare_nhdplus(nhd %>% filter(COMID %in% unlist(nhd_US)), 0, 0, 0, FALSE)
      #print (dim(keep))
      ifelse( exists("final") , final <- rbind(final, keep), final <- keep)
    }
    
    nhd <- nhd %>% filter(COMID %in% final$COMID)
  }
  return (nhd)
}

return_NonDendirtic <- function(VPU, nhd, nhdPath){
  # Returns non-dendritic flowlines for assignment of coastal/non-dendritic catchments/flowlines
  nhdAll <- readRDS(nhdPath) %>%
    filter(grepl(paste0("^",VPU,".*"), .data$VPUID))
    nhd_ND <- nhdAll %>% filter(!COMID %in% nhd$COMID) %>% filter(TerminalFl == 1)
    
  return(nhd_ND)
  }

Merge_hydReg <- function(feat, out_gpkg){
  if(needs_layer(out_gpkg, feat)) {
    feat_01 <- read_sf("cache/GF_01.gpkg", paste0(feat,"_01")) 
    feat_02 <- read_sf("cache/GF_02.gpkg", paste0(feat,"_02"))
    feat_03 <- read_sf("cache/GF_03.gpkg", paste0(feat,"_03"))
    feat_04 <- read_sf("cache/GF_04.gpkg", paste0(feat,"_04"))
    feat_05 <- read_sf("cache/GF_05.gpkg", paste0(feat,"_05"))
    feat_06 <- read_sf("cache/GF_06.gpkg", paste0(feat,"_06"))
    feat_07 <- read_sf("cache/GF_07.gpkg", paste0(feat,"_07"))
    feat_08 <- read_sf("cache/GF_08.gpkg", paste0(feat,"_08"))
    feat_09 <- read_sf("cache/GF_09.gpkg", paste0(feat,"_09"))
    feat_10U <- read_sf("cache/GF_10U.gpkg", paste0(feat,"_10U"))
    feat_10L <- read_sf("cache/GF_10L.gpkg", paste0(feat,"_10L"))
    feat_11 <- read_sf("cache/GF_11.gpkg", paste0(feat,"_11"))
    feat_12 <- read_sf("cache/GF_12.gpkg", paste0(feat,"_12"))
    feat_13 <- read_sf("cache/GF_13.gpkg", paste0(feat,"_13"))
    feat_14 <- read_sf("cache/GF_14.gpkg", paste0(feat,"_14"))
    feat_15 <- read_sf("cache/GF_15.gpkg", paste0(feat,"_15"))
    feat_16 <- read_sf("cache/GF_16.gpkg", paste0(feat,"_16"))
    feat_17 <- read_sf("cache/GF_17.gpkg", paste0(feat,"_17"))
    feat_18 <- read_sf("cache/GF_18.gpkg", paste0(feat,"_18"))

    allfeat <- list(feat_01, feat_02, feat_03, feat_04, feat_05, feat_06, feat_07, feat_08, feat_09,
                    feat_10L, feat_10U, feat_11, feat_12, feat_13, feat_14, feat_15, feat_16, feat_17,
                    feat_18)
                  
  
    all_sf<-do.call(rbind, allfeat)
    write_sf(all_sf, out_gpkg, feat)
  }
}

siteAtt <- function(hydReg){
  print ("gage")
  out_gpkg <- file.path("cache",paste0("GF_", hydReg, ".gpkg"))
  
  data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
  
  # Layers
  allPOIs <- read_sf(out_gpkg, paste0("POIS_", hydReg)) 
  gagesIII_pois <- read_sf(out_gpkg, paste0("GagesIII_", hydReg)) # GAGESIII POIs
  TE_pois <- read_sf(out_gpkg, paste0("TE_", hydReg)) # Thermoelectric POIs
  NID_pois <- read_sf(out_gpkg, paste0("NID_", hydReg)) # NID POIs
  nsegment_raw <- read_sf(out_gpkg, paste0("nsegment_raw_", hydReg)) # Minimally-sufficient network attributed with POI_ID
  nsegment <- read_sf(out_gpkg, paste0("Nsegment_", hydReg)) # Minimally-sufficient network dissolved by POI_ID
  WBs_VPU <- read_sf(out_gpkg, paste0("WB_", hydReg)) %>% st_drop_geometry() %>% mutate(COMID = as.integer(COMID))
  WBout_POIs <- read_sf(out_gpkg, paste0("WBout_", hydReg)) %>% st_drop_geometry()
  WBin_POIs <- read_sf(out_gpkg, paste0("WBin_", hydReg)) %>% st_drop_geometry()
  
  #*************************************************************************
  # Attribute gagesIII POIs with adjusted measure and drainage area information
  gagesIII <- read_sf(data_paths$gagesiii_points_path, "GagesIII_gages") %>%
    filter(grepl(paste0("^", hydReg, ".*"), .data$NHDPlusReg))

  # Load GAGESIII POIs from gpkg at cache dir
  gagesIIIPOIs <- gagesIII_pois %>% st_drop_geometry(.) %>% select(Type_GagesIII, COMID) %>%
    inner_join(gagesIII, by = "COMID") %>% select(Type_GagesIII, COMID, COMID_meas)

  # Derive segment measure attributes
  nsegment_gages <- st_drop_geometry(nsegment_raw) %>% left_join(gagesIIIPOIs, by = "COMID") %>%
    select(Hydroseq, Type_GagesIII, COMID_meas, LENGTHKM, POI_ID, TotDASqKM) %>%
    mutate(COMID_meas = 1 - replace_na(COMID_meas/100,0), Seg_meas = COMID_meas * LENGTHKM) %>% group_by(POI_ID) %>%
    summarize(HS = min(Hydroseq), Length = sum(LENGTHKM), Percent = sum(Seg_meas)/sum(LENGTHKM),
              TotDASqKM = max(TotDASqKM), Type_GagesIII =  first(na.omit(Type_GagesIII))) %>% filter(!is.na(Type_GagesIII))

  # Calcualte estimated proporational drainage area
  gagesIIIPOI_DA <- gagesIIIPOIs %>% left_join(nsegment_raw, by = "COMID") %>%
    select(Type_GagesIII, COMID, COMID_meas, AreaSqKM, TotDASqKM) %>%
    mutate(Inc_DA = (1-(COMID_meas / 100)) * AreaSqKM) %>%
    mutate(Total_DA = TotDASqKM - (AreaSqKM - Inc_DA))

  # PUt everything together
  gagesIII_stats <- gagesIIIPOI_DA %>% inner_join(nsegment_gages, by = "Type_GagesIII") %>%
    select(Type_GagesIII, COMID, Length, Percent, TotDASqKM.y) %>%
    rename(Seg_Length = Length, Seg_Measure = Percent, Adj_NHD_DA = TotDASqKM.y) %>%
    mutate(VPU = hydReg)

  # Compile data into GagesIII folder
  if (file.exists("data/GAGESIII_gages/gagesIII_MDA.rds")){
    # open file and bind to it
    curRDS_gages <- readRDS("data/GAGESIII_gages/gagesIII_MDA.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(gagesIII_stats)
    saveRDS(curRDS_gages, "data/GAGESIII_gages/gagesIII_MDA.rds")
  } else {
    saveRDS(gagesIII_stats, "data/GAGESIII_gages/gagesIII_MDA.rds")
  }
  #************************************************************************

  #************************************************************************
  # Read in TE shapefile
  print ("TE")
  TE_fac <- st_read(file.path(data_paths$TE_points_path, "/2015_TE_Model_Estimates_lat.long_COMIDs.shp")) %>%
    filter(COMID %in% TE_pois$COMID)

  # Cast TE as points and project to USGS albers
  TE_fac_alb <- TE_fac %>% st_cast('POINT') %>% st_transform(., 6703)

  # Subset segs based on TE and project to USGS albs
  TE_seg_pnts <- nsegment %>% filter(POI_ID %in% TE_fac_alb$COMID) %>% st_transform(., 6703) %>% group_by(POI_ID) %>%
    st_line_merge(.) %>% st_cast("POINT")

  #for (POI_ID in TE_fac_alb$COMID){
  TE_fac_alb$TE_measure <- unlist(lapply(1:nrow(TE_fac_alb), function(x){
    COM <- TE_fac_alb$COMID[x]
    fac <- TE_fac_alb %>% filter(COMID == COM)
    seg <- TE_seg_pnts %>% filter(POI_ID == COM)

    # Get index of nearest vertex on segment and rebuild two line IDs, then get length of each
    nearVert <- st_nearest_feature(fac, seg)
    # Pull out the line upstream and downstream of snapping point
    Upline <- seg[1:nearVert,] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)
    Downline <- seg[nearVert:nrow(seg),] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)

    # Return the measure
    return (Upline / (Upline + Downline))
  }))

  # Derive proportional total (probably need proporational incremental ratio too)
  TE_data <- st_drop_geometry(TE_fac_alb) %>%
    left_join(st_drop_geometry(nsegment_raw) %>% select(COMID, AreaSqKM, TotDASqKM), by = "COMID") %>%
    mutate(Inc_DA = TE_measure * AreaSqKM, Total_DA = TotDASqKM - (AreaSqKM - Inc_DA)) %>%
    select(EIA_PLANT_, LATITUDE, LONGITUDE, COMID, TE_measure, Inc_DA, Total_DA)

  TE_data_upstream <- TE_data %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("COMID" = "To_POI_ID")) %>%
    group_by(COMID) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " "))

  TE_data_final <- TE_data_upstream %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("COMID" = "POI_ID")) %>%
    group_by(COMID) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " ")) %>%
    rename(Downstream_POI = To_POI_ID) %>% select(-POI_ID) %>% distinct() %>% mutate(VPU = hydReg)

  if (file.exists("data/TE_points/TE_adj.rds")){
    # open file and bind to it
    curRDS_TE <- readRDS("data/TE_points/TE_adj.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(TE_data_final)
    saveRDS(curRDS_TE, "data/TE_points/TE_adj.rds")
  } else {
    saveRDS(TE_data_final, "data/TE_points/TE_adj.rds")
  }
  #************************************************************************
  
  #************************************************************************
  # Read in NID shapefile
  print ("NID")
  NID_fac <- read.csv(file.path(data_paths$NID_points_path,"/NID_attributes_20170612.txt"), stringsAsFactors = F) %>%
    filter(FlowLcomid %in% NID_pois$COMID)

  NID_snap_alb <- readRDS(file.path(data_paths$NID_points_path,"/NAWQA_NID_snap.rds")) %>%
    filter(Source_FeatureID %in% NID_fac$Source_FeatureID) %>% st_transform(., 6703) %>%
    inner_join(NID_fac, by = "Source_FeatureID")

  # Subset segs based on TE and project to USGS albs
  NID_seg_pnts <- nsegment %>% filter(POI_ID %in% unique(NID_fac$FlowLcomid)) %>% st_transform(., 6703) %>% group_by(POI_ID) %>%
    st_line_merge(.) %>% st_cast("POINT")

  #options(warn = -1)
  NID_snap_alb$NID_measure <- unlist(lapply(1:nrow(NID_snap_alb), function(x){
    fac <- NID_snap_alb[x,]
    COM <- fac$FlowLcomid
    seg <- NID_seg_pnts %>% filter(POI_ID == COM)

    # Get index of nearest vertex on segment and rebuild two line IDs, then get length of each
    nearVert <- st_nearest_feature(fac, seg)
    # Pull out the line upstream and downstream of snapping point
    Upline <- seg[1:nearVert,] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)
    Downline <- seg[nearVert:nrow(seg),] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)

    # Return the measure
    return (Upline / (Upline + Downline))
  }))
  #options(warn = oldw)

  # Derive proportional total (probably need proporational incremental ratio too)
  NID_data <- st_drop_geometry(NID_snap_alb) %>%
    left_join(st_drop_geometry(nsegment_raw) %>% select(COMID, AreaSqKM, TotDASqKM), by = c("FlowLcomid" = "COMID")) %>%
    mutate(Inc_DA = NID_measure * AreaSqKM, Total_DA = TotDASqKM - (AreaSqKM - Inc_DA)) %>%
    select(Source_FeatureID, VPU.x, Outlet, FlowLcomid, WBAreaSQKM, NIDSASQKM, NLCDSASQKM, MAX_SA, MultiMain, Catchcomid, FalconeDate,
           YEAR_COMPLETED, MAX_DISCHARGE, NID_STORAGE, NID_measure, Inc_DA, Total_DA)

  NID_data_upstream <- NID_data %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("FlowLcomid" = "To_POI_ID")) %>%
    group_by(FlowLcomid) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " "))

  NID_data_final <- NID_data_upstream %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("FlowLcomid" = "POI_ID")) %>%
    group_by(FlowLcomid) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " ")) %>%
    rename(Downstream_POI = To_POI_ID) %>% select(-POI_ID) %>% distinct() %>% mutate(VPU = hydReg)

  if (file.exists("data/NID_points/NID_adj.rds")){
    # open file and bind to it
    curRDS_NID <- readRDS("data/NID_points/NID_adj.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(NID_data_final)
    saveRDS(curRDS_NID, "data/NID_points/NID_adj.rds")
  } else {
    saveRDS(NID_data_final, "data/NID_points/NID_adj.rds")
  }
  
  #****************************************************************************
  # Get information on all inflow/outflows in nsegment network
  # index the  WBAREACOMI to In WB POIs
  print("waterbodies")
  
  WBinPOIs_tab <- WBin_POIs %>% select(COMID) %>%
    left_join(st_drop_geometry(nsegment_raw), by = "COMID") %>% select(COMID, DnHydroseq)%>%
    left_join(st_drop_geometry(nsegment_raw), by = c("DnHydroseq" = "Hydroseq")) %>%
    mutate(COMID = COMID.x) %>% select(COMID, WBAREACOMI) 
  
  # Noticed there are some "NHDArea" features coded into the WBARECOMI field (prob mistakes)
  # Should we QAQC in previous step and remove these?
  wb_pol_intab <- WBs_VPU %>% inner_join(WBinPOIs_tab, by = c("COMID" = "WBAREACOMI")) %>%
    select(COMID, COMID.y) %>% rename(seg_COMID = COMID.y) %>%
    group_by(COMID) %>% summarize(num_poi_inlets = n(), Upstream_POIs = paste0(unique(seg_COMID), collapse = " "))
  
  # index the WBAREACOMI to out WB POIs
  # 3 WBoutPOIs have no WBareaCOMID
  WBoutPOIs_tab <- WBout_POIs %>% select(COMID) %>% left_join(st_drop_geometry(nsegment_raw), by = "COMID") %>%
    select(COMID, WBAREACOMI) 
  
  wb_pol_outtab <- WBs_VPU %>% inner_join(WBoutPOIs_tab, by = c("COMID" = "WBAREACOMI")) %>%
    select(COMID, COMID.y)  %>% rename(seg_COMID = COMID.y) %>% 
    group_by(COMID) %>% summarize(num_poi_outlets = n(), outlet = paste0(unique(seg_COMID), collapse = " "))
  
  # WBs that have multiple outlets  
  multiOutWBs <- wb_pol_outtab %>% filter(num_poi_outlets > 1)
  
  wb_pol_tab <- wb_pol_outtab %>% left_join(wb_pol_intab, by = "COMID")
  
  # write to csv file
  wbmore_tab <- WBs_VPU %>% select(-c(FDATE, RESOLUTION, Shape_Length, Shape_Area, GNIS_ID, REACHCODE, FTYPE, FCODE, 
                                      PurpCode, PurpDesc, MeanDCode, ONOFFNET, LakeArea)) %>%
    inner_join(wb_pol_tab, by = "COMID")
  #write.csv(wbmore_tab, paste0("cache/SiteInfo/WBPOIs_r", hydReg, "_tab.csv"))
  
  #****************************************************************************
  # Create table of spatial assocation of waterbodies with gages, NID, and TE
  wbpoi_with_gage <- allPOIs %>% filter((!is.na(Type_WBIn) | !is.na(Type_WBOut) & !is.na(Type_GagesIII)))
  wbpoi_with_NID <- allPOIs %>% filter((!is.na(Type_WBIn) | !is.na(Type_WBOut) & !is.na(Type_NID)))
  wbpoi_with_TE <- allPOIs %>% filter((!is.na(Type_WBIn) | !is.na(Type_WBOut)) & !is.na(Type_TE))
  
  # determine gages by waterbody (using segs WBAREACOMI); only wb (not nhdrea)
  wb_gage <- gagesIII_pois %>% left_join(st_drop_geometry(nsegment_raw), by = "COMID") %>%
    select(COMID, WBAREACOMI, Type_GagesIII) %>% filter(WBAREACOMI > 0, WBAREACOMI %in% WBs_VPU$COMID) 
  
  # Waterbody inlets that are gages
  wb_ingage_tab <- st_drop_geometry(gagesIII_pois) %>% filter(Type_WBIn == 1) %>%
    left_join(WBinPOIs_tab, by = "COMID") %>% select(COMID, Type_GagesIII, WBAREACOMI) %>%
    rename(WBin_Gage = Type_GagesIII)
  
  # find if gage at upstream segment of a WBin POI
  wb_upgage_tab <-  st_drop_geometry(gagesIII_pois) %>% select(COMID, Type_GagesIII) %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("COMID" = "POI_ID")) %>%
    inner_join(WBinPOIs_tab, by = c("To_POI_ID" = "COMID")) %>% group_by(WBAREACOMI) %>% 
    summarize(num_us_gages = n(), Upstream_gages = paste0(unique(Type_GagesIII), collapse = ", "))
  
  # WB outlets 
  wb_outgage_tab <- gagesIII_pois %>% filter(Type_WBOut == 1) %>% 
    inner_join(WBoutPOIs_tab, by = "COMID") %>% select(COMID, WBAREACOMI, Type_GagesIII) %>%
    rename(WBout_Gage = Type_GagesIII)
  
  # find if gage on downstream segment of a WBout POI
  wb_dngage_tab <- WBout_POIs %>% select(COMID) %>% inner_join(WBoutPOIs_tab, by = "COMID") %>%
    left_join(nsegment %>% select(POI_ID, To_POI_ID), by = c("COMID" = "POI_ID")) %>%
    inner_join(gagesIII_pois, by = c("To_POI_ID" = "COMID")) %>% rename(Downstream_gage = Type_GagesIII)
  
  # Read in NID sites for the VPU (and get unique COMIDS)there may be more than one per COMID
  NID_fac <- read.csv(paste0(data_paths$NID_points_path,"/NID_attributes_20170612.txt"), stringsAsFactors = F) %>%
    select(-c(ReachCode, Measure, Source_FeatureID, VPU, SNAPTYPE, comment, MultiNIDs, NHDnet_Err, ONNHDPLUS,
              MultiMain, Catchcomid)) %>% filter(FlowLcomid %in% NID_POIs$COMID) %>% 
    left_join(nsegment_raw %>% select(COMID, WBAREACOMI), by = c("FlowLcomid" = "COMID"))
  
  # gen list of WB NID, put on a table, there may be more than one per WB
  wb_NID <- NID_fac %>% inner_join(WBs_VPU, by = c("WBAREACOMI" = "COMID")) %>% group_by(WBAREACOMI) %>% 
    summarize(Outlet = max(Outlet), WBAreaSQKM = max(WBAreaSQKM, na.rm = T), NIDSASQKM = max(NIDSASQKM, na.rm = T),
              NLCDSASQKM = max(NLCDSASQKM, na.rm = T), MAX_SA = max(MAX_SA, na.rm = T), FalconeDate = max(FalconeDate, na.rm = T),
              NLCDASADIST = max(NLCDSADIST), NIDID = paste0(unique(NIDID), collapse = " "), 
              YEAR_COMPLETED = max(YEAR_COMPLETED, na.rm = T), MULT_YEARS = paste0(unique(YEAR_COMPLETED), collapse = " "),
              MAX_DISCHARGE = max(MAX_DISCHARGE), NID_STORAGE = max(NID_STORAGE, na.rm = T)) %>% 
    mutate_all(funs(stringr::str_replace_all(., "-Inf", "NA"))) %>% mutate(WBAREACOMI = as.integer(WBAREACOMI))
  
  # Do we really need to read the facilities in again
  TE_fac <- read_sf(data_paths$TE_points_path, "2015_TE_Model_Estimates_lat.long_COMIDs") %>% 
    mutate(COMID=as.integer(COMID)) %>% filter(COMID > 0) %>%
    inner_join(., select(st_drop_geometry(nhd), c(COMID, WBAREACOMI, VPUID)), by = "COMID") %>%
    filter(grepl(paste0("^", hydReg, ".*"), .data$VPUID)) %>% switchDiv(., nhd)
  
  # gen list of WB NID, put on a table, there may be more than one per WB
  wb_TE <- TE_fac %>% inner_join(WBs_VPU, by = c("WBAREACOMI" = "COMID")) %>% 
    select(EIA_PLANT_, COMID, WBAREACOMI) 
  
  wb_finTable <- wbmore_tab %>% left_join(wb_NID, by = c("COMID" = "WBAREACOMI")) %>% mutate(VPU = hydReg)
  
  # create WB table: WB_r"xx_info
  wb_POI_tab <- WBs_VPU %>% rename(WBAREACOMI=COMID) %>%
    select(WBAREACOMI) %>%
    left_join(wb_outgage_tab %>% select(WBAREACOMI, WBout_Gage)) %>%
    left_join(wb_dngage_tab %>% select(WBAREACOMI, Downstream_gage)) %>%
    left_join(wb_ingage_tab %>% select(WBAREACOMI, WBin_Gage)) %>% #fix to be able to handle multiple gages
    left_join(wb_upgage_tab %>% select(WBAREACOMI, Upstream_gages)) %>%
    left_join(wb_NID) %>% select(WBAREACOMI, Outlet, NIDID) %>%
    mutate(VPU = hydReg)
  
  if (file.exists("data/NHDPlusNationalData/wbod_stats.rds")){
    # open waterbody file and bind to it
    curRDS_wbod <- readRDS("data/NHDPlusNationalData/wbod_stats.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(wb_finTable)
    saveRDS(curRDS_wbod, "data/NHDPlusNationalData/wbod_stats.rds")
    
    # open POI file and bind to it
    curRDS_POIwbod <- readRDS("data/NHDPlusNationalData/wbod_POIs.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(wb_POI_tab)
    
  } else {
    saveRDS(wb_finTable, "data/NHDPlusNationalData/wbod_stats.rds")
    saveRDS(wb_POI_tab, "data/NHDPlusNationalData/wbod_POIs.rds")
  }
}

# make sf1 compatible with sf2
st_compatibalize <- function(sf1, sf2) {
  sf1 <- st_transform(sf1, st_crs(sf2))
  
  g <- attr(sf1, "sf_column")
  gp <- attr(sf2, "sf_column")
  names(sf1)[names(sf1) == g] <- gp
  attr(sf1, "sf_column") <- gp
  sf1
}