NHD_navigate.R

#' Network navigation for upstream/downstream from a COMID of interest
#' @param inCOM  (list) list of input COMIDs 
#' @param nhdDF (sf data.frame) (data frame) valid data frame of NHD flowlines
#' @param withTrib (logical) flag for if the upstream navigation should include tributaries
#              or stick to mainstem level path
#
#' @return (list) list of COMIDs upstream of point
NetworkNav <- function(inCom, nhdDF, withTrib){
  
  if (missing(withTrib)){
    seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
  } else {
    seg <- nhdplusTools::get_UT(nhdDF, inCom)
  }
  return(seg)
}

#' Creates raw and dissolved segment layers with necessaary
#         upstream/downstream routing attribution
#'  @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#'  @param POIs  (sf data.frame) Existing POIs
#' 
#' @return (sf data.frame) data.frame of segments connecting POIs attributed
#'         with POI_ID for each upstream flowpath
segment_increment <- function(nhdDF, POIs){

  ptm<-proc.time()
  
  seg_POIs <- arrange(POIs, desc(LevelPathI), desc(Hydroseq)) %>%
    select(COMID, Hydroseq, LevelPathI) %>%
    group_by(LevelPathI) %>%
    # These next two levels arrange POIs hydrologically along the 
    #       level path in order to meet the requirements of the code below
    mutate(id = row_number(), 
           num = if(n() > 1) id[1L] + row_number()-1 else id) %>%
    ungroup()
  
  # Add an empty field for POI_Id population
  nhdDF <- mutate(nhdDF, POI_ID = 0)
  
  POI_ID_assign <- function(i, seg_POIs, nhd){
    ##########################################
    # Populates POI_ID per segment
    # 
    # Arguments:
    #   i : (integer) iterator
    #   seg_POIs : (data frame) POI data frame
    #   nhd : (dataframe) flowlines data frame (no geometry)
    #  
    # Returns:
    #   nhd_sub : (data frame) raw segments with POI_ID populated
    #             
    ##########################################
    library(dplyr)
    
    # If POI is most upstream POI on levelpath
    if (seg_POIs$num[i] == 1){
      # Assign POI_ID
      nhd_sub <- filter(nhd, Hydroseq >= seg_POIs$Hydroseq[i] & 
                          LevelPathI == seg_POIs$LevelPathI[i]) %>% 
        mutate(POI_ID =  seg_POIs$COMID[i])
      # or as you travel downstream on set of POIs below level path
    } else {
      # Assign POI_ID
      nhd_sub <- filter(nhd, LevelPathI == seg_POIs$LevelPathI[i] & 
                          Hydroseq >= seg_POIs$Hydroseq[i] & Hydroseq < seg_POIs$Hydroseq[i-1]) %>% 
        mutate(POI_ID = seg_POIs$COMID[i])
    }
    # return assigned flowlines
    return(select(nhd_sub, LevelPathI, Hydroseq, COMID, POI_ID) %>%
             filter(POI_ID != 0))
  }
  
  library(parallel)
  library(dplyr)
  clust <- makeCluster(4)
  POI_list <- parLapply(clust, c(1:nrow(seg_POIs)), POI_ID_assign, seg_POIs, st_drop_geometry(nhdDF))
  #POI_list <- lapply(c(1:nrow(seg_POIs)), POI_ID_assign, seg_POIs, st_drop_geometry(nhdDF))
  stopCluster(clust)
  
  inc_segs <- data.table::rbindlist(POI_list)
  
  print(proc.time()-ptm)
  return(inc_segs)
}

#' Creates finalized segments and routing
#'  @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#'  @param routing_fix  (sf data.frame) any additional routing fixes
#' 
#' @return (sf data.frame) data.frame of segments
segment_creation <- function(nhdDF, routing_fix = NULL){ 
  
  if(!"StartFlag" %in% names(nhdDF)) {
    nhdDF$StartFlag <- ifelse(nhdDF$Hydroseq %in% nhdDF$DnHydroseq, 0, 1)
  }
  
  in_segs <- filter(nhdDF, !is.na(POI_ID))
  
  # If there are routing fixes to account for if a POI with a DA of 0 is moved upsream or downstream
  if (is.data.frame(routing_fix)){
    routing_fix <- routing_fix %>%
      rename(COMID = oldPOI, new_COMID = COMID)
    
    # Above we generated the network using the initial set of POIs; here we crosswalk over the old COMIDs to the new
    nhd_fix <- nhdDF %>%
      left_join(routing_fix %>%
                   select(COMID, new_COMID), by = c("POI_ID" = "COMID")) %>%
      mutate(POI_ID = ifelse(is.na(new_COMID), POI_ID, new_COMID)) %>%
      filter(!POI_ID %in% routing_fix$COMID) %>%
      select(-new_COMID)
    
    in_segs <- filter(nhd_fix, !is.na(POI_ID))
  }
  
  # Dissolve flowlines to aggregated segments
  nsegments <- filter(in_segs, !is.na(POI_ID)) %>%
    group_by(POI_ID) %>%
    #arrange(desc(LevelPathI), desc(Hydroseq)) %>%
    summarize(TotalLength = sum(LENGTHKM),TotalDA = max(TotDASqKM), HW = max(StartFlag),
              do_union = FALSE) %>%
    #st_cast("MULTILINESTRING")  %>%
    inner_join(st_drop_geometry(filter(in_segs, minNet == 1)) %>%
                 select(COMID, Hydroseq, DnHydroseq), by = c("POI_ID" = "COMID"))
  
  # produce a short data frame for populating TO_POI for downstream segment
  to_from <- filter(st_drop_geometry(in_segs)) %>%
    left_join(filter(st_drop_geometry(nhdDF), !is.na(POI_ID)) %>% 
                select(COMID, Hydroseq, POI_ID), by = c("DnHydroseq" = "Hydroseq")) %>%
    select(COMID.x, Hydroseq, DnHydroseq, POI_ID.y) %>%
    rename(To_POI_ID = POI_ID.y) 
  
  # Add To_POI_ID to dissolved segments
  nsegments_fin <- nsegments %>% 
    left_join(select(to_from, COMID.x, To_POI_ID), by = c("POI_ID" = "COMID.x")) %>%
    select(POI_ID, TotalLength, TotalDA, HW, To_POI_ID) 
  
  return(list(diss_segs = nsegments_fin, raw_segs = in_segs))
}

#' Moves POI Upstream or downstream if it falls on COMID
#       of flowline with no corresponding catchment
#'  @param POIs_wgeom (sf data.frame) POIs
#'  @param nhdDF  (sf data.frame) valid data frame of NHD flowlines
#' 
#' @return (sf data.frame) data.frame of POIs with new COMID associated
DS_poiFix <- function(POIs_wgeom, nhd){
  nhd <- distinct(nhd)
  POIs <- st_drop_geometry(POIs_wgeom) %>%
    arrange(COMID) %>%
    filter(nexus == FALSE)

  # DF of downstream segment
  tocomDF <- select(st_drop_geometry(nhd), COMID, Hydroseq, TotDASqKM,
                    DnHydroseq, WBAREACOMI) %>%
    inner_join(select(st_drop_geometry(nhd), COMID_ds = COMID, Hydroseq, 
                      WBAREACOMI_down = WBAREACOMI, totda_ds = TotDASqKM), 
               by = c("DnHydroseq" = "Hydroseq")) %>%
    inner_join(select(st_drop_geometry(nhd), COMID_us = COMID, DnHydroseq,
                      WBAREACOMI_up = WBAREACOMI, totda_us = TotDASqKM),
               by = c("Hydroseq" = "DnHydroseq"))
  
  # Find segments with POIs where there is no corresponding catchment that are not terminal
  unCon_fl <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)# & Hydroseq != TerminalPa)
  unCon_POIs <- filter(POIs, COMID %in% unCon_fl$COMID)
  
  # Get specific fixes for waterbody inlets and outlets
  wbout <- filter(unCon_POIs, !is.na(Type_WBOut)) %>%
    inner_join(tocomDF, by = "COMID") %>%
    mutate(nonrefactor = ifelse(WBAREACOMI %in% WBAREACOMI_up, COMID_ds, 0),
           new_POI = COMID_us)

  wb_pois <- filter(unCon_POIs, !is.na(Type_WBIn)) %>%
    inner_join(tocomDF, by = "COMID") %>%
    mutate(nonrefactor = ifelse(WBAREACOMI < 0, COMID_ds, 0),
           new_POI = COMID_us) %>%
    rbind(wbout) %>%
    select(-c(nexus, Hydroseq, TotDASqKM, DnHydroseq, WBAREACOMI, COMID_ds,
              WBAREACOMI_down, totda_ds, COMID_us, WBAREACOMI_up, totda_us)) %>%
    rename(COMID = new_POI, oldPOI = COMID)

  # The rest can be resolved with drainage are ratio
  unCon_POIs <- filter(unCon_POIs, !COMID %in% wb_pois$oldPOI)

  poi_fix <- as.data.frame(do.call("rbind", lapply(unCon_POIs$COMID, movePOI_NA_DA, st_drop_geometry(nhd)))) %>%
    inner_join(POIs, by = c("oldPOI" = "COMID")) %>%
    inner_join(select(st_drop_geometry(nhd), COMID), by = c("oldPOI" = "COMID")) %>%
    select(-c(AreaSqKM, DnHydroseq, nexus, TotDASqKM)) %>%
    distinct() %>%
    bind_rows(wb_pois)
  
  # Fold in new POIs with existing POIs so all the "Type" attribution will carry over
  # using the minimum will ensure correct downstream hydrosequence gets carried over
  poi_orig <- filter(POIs, COMID %in% poi_fix$COMID) %>%
    bind_rows(poi_fix) %>%
    select(-oldPOI)

  list_df <- dplyr::group_by(poi_orig, COMID) |>
    group_split()
  
  compact <- function(l) {
    if(nrow(l) == 1) return(as.data.frame(l))
    lapply(names(l), \(x) {
      out <- unique(l[[x]][!is.na(l[[x]])])
      if(!length(out)) NA 
      else if(length(out) == 1) out 
      else {
        cat(paste("duplicate id for", unique(l$COMID),
                  "column", x, 
                  "values", paste(out, collapse = ", "), 
                  "using", out[1]), file = "POI-issues.log")
        out[1]
      }
    }) |> 
      setNames(names(l)) |>
      as.data.frame()
  }
  
  poi_merged <- bind_rows(lapply(list_df, compact))

  # # Combine POI information together for redundant pois  
  # poi_merged <- poi_orig %>% 
  #   select(-c(nexus, AreaSqKM, oldPOI, DnHydroseq, TotDASqKM)) %>%
  #   group_by(COMID) %>%
  #   summarise_each(funs(toString(na.omit(.)))) 
  # is.na(poi_merged) <- poi_merged == ""
  
  # Join new POI COMIDs and geometry with the old Type fields
  fin_POIs <- poi_merged %>%
    arrange(COMID) %>%
    bind_cols(get_node(filter(nhd, COMID %in% .$COMID) %>% arrange(COMID), position = "end")) %>%
    st_sf() %>%
    st_compatibalize(., POIs_wgeom)
  
  return (list(xWalk = poi_fix, new_POIs = fin_POIs))
}


#' Move POIs that fall on flowlines with no catchment upstream/downstream
#     to adjacent flowline with most similar total drainage area. Called from 
#     DS_poi_fix function above
#'  @param poi_fix (data.frame) POI data set of COMIDs to be changed
#'  @param nhdDF  (sf data.frame) valid data frame of NHD flowlines
#' 
#' @return (data frame) DF of POIs with new COMID associated
movePOI_NA_DA <- function(poi_fix, nhdDF){
  #print(poi_fix)
  nhdDF <- distinct(nhdDF)
  
  # Closest POI/US/DS
  up_segs <- unique(nhdplusTools::get_UM(nhdDF, poi_fix, sort=T)) 
  seg2fix <- filter(nhdDF, COMID == poi_fix) %>%
    distinct()
  
  # Sorted results and filter out all flowlines w/o catchments
  upstuff <- filter(nhdDF, COMID %in% unlist(up_segs)) %>% 
    arrange(Hydroseq) %>%
    filter(AreaSqKM > 0)
  
  down_segs <- unique(nhdplusTools::get_DM(nhdDF, poi_fix, sort=T))
  downstuff <- filter(nhdDF, COMID %in% unlist(down_segs)) %>% 
    arrange(Hydroseq)%>%
    filter(AreaSqKM > 0)
  
  # combine into one dataframe, select up/downstream seg with least change in total drainage area
  near_FL <- rbind(select(upstuff, COMID, TotDASqKM, AreaSqKM) %>% 
                    slice(1), 
                    select(downstuff, COMID, TotDASqKM, AreaSqKM) %>% 
                    slice(1))
  
  # If 1 or other adjacent flowlines are coupled with a catchment
  if (sum(near_FL$AreaSqKM) > 0){
    new_POI <- near_FL[(which.min(abs(seg2fix$TotDASqKM - near_FL$TotDASqKM))),] #near_FL$COMID
    new_POI$oldPOI <- poi_fix
    new_POI$DnHydroseq <-seg2fix$DnHydroseq
  } else {
    # Remove POI if not catchment associated with flowlines upstream or downstream
    print (poi_fix)
    print ("US and DS flowlines also have no catchment, removing POI")
    new_POI <- NA
  }
  return(new_POI)
}

#' Collaposes/merges POIs together based on drainage area ratio and data theme
#'  @param out_gpkg (geopackage) output geopackage to write intermediate results to
#'  @param POIs  (sf data.frame) Original  POIs
#'  @param Seg  (sf data.frame) dissolved segments
#'  @param Seg2  (sf data.frame) (data frame) raw segments/flowlines
#'  @param exp  (string) Type of thematic POI being moved round
#'  @param DA_diff  (float) threshold of drainage area difference to
#             consider when merging two POIs
#'  
#' @return (list of data.frames) 1 - New set of POIs
#          2/3 - correpsonding segments (both raw and dissolved)
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
  
  nexus <- filter(POIs, nexus == TRUE)
  
  POIs <- POIs %>% mutate_if(is.numeric, function(x) ifelse(is.infinite(x), NA, x)) %>%
    filter(nexus != 1)
  
  # downstream POIs and their drainage area ratios and incremental drainage areas
  segs_down <- Seg %>% 
    inner_join(select(st_drop_geometry(.), c(POI_ID, TotalDA, TotalLength)), 
               by = c("To_POI_ID" = "POI_ID")) %>%
    inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>%
    # Keep WBOut and In to preserve NHD Waterbody geometry, Keep Confluences and TE where they are
    filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_TE), all_vars(is.na(.))) %>%
    # Select POIs within the correct drainage area ratio and fit the size criteria
    mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.y - TotalDA.x) %>% 
    # If the drainage area ratio is within acceptable range (1.05)
    filter(DAR < (1 + DA_diff), POI_ID != To_POI_ID) %>% 
    # select and re-order
    select(POI_ID, HW, To_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID, 
           Type_WBIn, Type_WBOut, Type_Conf, Type_Term, Type_Elev, Type_Travel, nexus, DAR, IncDA)

  # upstream POIs and their drainage area ratios and incremental drainage areas
  up_POIs <- filter(POIs, is.na(Type_Conf))
  segs_up <- Seg %>% inner_join(select(filter(st_drop_geometry(.), POI_ID %in% up_POIs$COMID), 
                                       c(POI_ID, To_POI_ID, TotalDA, TotalLength)), 
                                by = c("POI_ID" = "To_POI_ID")) %>%
    inner_join(st_drop_geometry(POIs), by = c("POI_ID" = "COMID")) %>% 
    rename(From_POI_ID = POI_ID.y) %>%
    # Don't want to move outlets upstream, don't move confluences upstream
    filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_TE), all_vars(is.na(.))) %>%
    # Select POIs within the correct drainage area ratio and fit the size criteria
    mutate(DAR = TotalDA.y/TotalDA.x, IncDA = TotalDA.x - TotalDA.y) %>% 
    # If the drainage area ratio is within acceptable range (0.95)
    filter(DAR >= (1 - DA_diff), POI_ID != From_POI_ID) %>% 
    # select and re-order
    select(POI_ID, HW, From_POI_ID, Type_HUC12, Type_Gages, Type_TE, Type_NID, 
           Type_WBIn, Type_WBOut, Type_Conf, Type_Term, Type_Elev, Type_Travel, nexus, DAR, IncDA)
  
  # Filter by POI Type, this removes the POI theme 
  Types <- c("Type_Gages", "Type_HUC12", "Type_TE", "Type_Conf", "Type_NID", "Type_WBIn", "Type_WBOut",
             "Type_Term", "Type_Elev", "Type_Travel")
  Types <- Types[Types != exp]
  
  # Gages coupled with confluences, waterbodies, do not move these
  static_POIs <- POIs %>% 
    filter(!is.na(Type_Gages) & (!is.na(Type_Conf) | !is.na(Type_WBOut) | !is.na(Type_WBIn)))
  
  # compiled list of POIs to move up or down
  seg_choices <- filter(POIs, !COMID %in% static_POIs) %>%
    left_join(st_drop_geometry(segs_down) %>%
                select(POI_ID, To_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
    left_join(st_drop_geometry(segs_up) %>%
                select(POI_ID, From_POI_ID, DAR, IncDA), by = c("COMID" = "POI_ID")) %>%
    filter(!is.na(To_POI_ID) | !is.na(From_POI_ID)) %>%
    filter(IncDA.x < min_da_km | IncDA.y < min_da_km)
  
  # We only want to move POIs that are not coupled within another theme
  seg_sub <- seg_choices %>% 
    filter_at(Types, all_vars(is.na(.))) %>%
    # Don't want to move the static POIS either
    filter(!COMID %in% static_POIs$COMID) %>%
    select (COMID, To_POI_ID, DAR.x, IncDA.x, From_POI_ID, DAR.y, DAR.y, IncDA.y) %>%
    st_sf()
  
  # Commented for QAQCing 
  #st_write(seg_sub, out_gpkg, "Seg_choices", delete_layer = F, append = T)
  
  # Data frame to indicate whether POIs should be grouped upstream or downstream
  segsub_dir <- mutate(seg_sub, DirDA = ifelse(is.na(DAR.y), "Down", 
                          ifelse(is.na(DAR.x), "Up",
                                 ifelse(IncDA.x < IncDA.y, "Down", "Up"))),
           Move = ifelse(DirDA == "Down", To_POI_ID, From_POI_ID))
  
  # POIs whose moves correspond to each others locations
  Mult <- select(segsub_dir, -c(To_POI_ID, From_POI_ID)) %>% #select(segsub_dir, COMID, Move) %>% 
    inner_join(st_drop_geometry(.), by = c("Move" = "COMID"), suffix = c(".A", ".B")) 
  
  # Basicly just chose one of the pair in Mult.
  segsub_dir <- filter(segsub_dir, !COMID %in% Mult$Move)
  
  #1 - POIs that need to be moved downstream
  move_POI <- filter(POIs, COMID %in% segsub_dir$COMID)
  # POIs that are stationary and their infomration will be merged with upstream  POIs
  stationary_POI <- filter(POIs, !COMID %in% move_POI$COMID)
  
  # Final Set to be merged with that don't involve either sets above
  #FinalPOI <- POIs %>% filter(!COMID %in% SegSub_Dir$COMID & !COMID %in% SegSub_Dir$Move) %>%
  #  mutate(merged_COMID = NA)
  
  #2 - Join SegSub assignment to POIs to bring over POI attributes
  movedownPOI_withatt <- select(st_drop_geometry(segsub_dir), COMID, DirDA, Move) %>%
    left_join(st_drop_geometry(move_POI), by = "COMID")  
  
  # Join SegSub mod to downstream POIs, COMID.x is the final COMID
  merged_POIs <- stationary_POI %>%
    inner_join(movedownPOI_withatt,
               by = c("COMID" = "Move"), suffix = c(".x", ".y")) %>%
  #merged_POIs[mapply(is.infinite, merged_POIs)] <- NA
  #merged_POIs <- merged_POIs %>%
    #mutate_all(., list(~na_if(.,""))) %>%
    # Don't overwrite existing gages or HUC12s
    filter(is.na(!!as.symbol(paste0(exp, ".x")))) %>%
    # Bring over relevant attributes
    mutate(Type_HUC12 = ifelse(is.na(Type_HUC12.y), Type_HUC12.x, Type_HUC12.y)) %>%
    mutate(Type_Gages = ifelse(is.na(Type_Gages.x), Type_Gages.y, Type_Gages.x)) %>%
    mutate(Type_WBIn = ifelse(is.na(Type_WBIn.y), Type_WBIn.x, Type_WBIn.y)) %>% 
    mutate(Type_WBOut = ifelse(is.na(Type_WBOut.y), Type_WBOut.x, Type_WBOut.y)) %>%
    mutate(Type_TE = ifelse(is.na(Type_TE.y), Type_TE.x, Type_TE.y)) %>%
    mutate(Type_NID = ifelse(is.na(Type_NID.y), Type_NID.x, Type_NID.y)) %>%
    mutate(Type_Conf = ifelse(is.na(Type_Conf.y), Type_Conf.x, Type_Conf.y)) %>%
    mutate(Type_Term = ifelse(is.na(Type_Term.y), Type_Term.x, Type_Term.y)) %>%
    mutate(Type_Elev = ifelse(is.na(Type_Elev.y), Type_Elev.x, Type_Elev.y)) %>%
    mutate(Type_Travel = ifelse(is.na(Type_Travel.y), Type_Travel.x, Type_Travel.y)) %>%
    mutate(oldPOI = COMID.y, to_com = NA) %>% 
    select(COMID, Type_HUC12, Type_Gages, Type_WBIn, Type_WBOut, 
           Type_TE, Type_NID, Type_Conf, Type_Term, 
           Type_Elev, Type_Travel, DirDA, nexus = nexus.x, oldPOI = COMID.y) %>%
    group_by(COMID) %>%
    summarize_all(~paste(unique(na.omit(.)), collapse = ',')) %>%
    st_sf()
    
  print(colnames(POIs))
  print(colnames(merged_POIs))
  
  # Bring new POI data over to new data
  fin_POIs <- filter(POIs, !COMID %in% c(merged_POIs$oldPOI, merged_POIs$COMID)) %>% 
    rbind(merged_POIs %>% select(-c(DirDA, oldPOI))) %>%
    rbind(nexus)
    #inner_join(st_drop_geometry(merged_POIs) %>% select(COMID, Type_Gages_A), by = "COMID") %>%
    #mutate(Type_Gages = ifelse(!is.na(Type_Gages_A), Type_Gages_A, Type_Gages)) %>% select(-Type_Gages_A)
  
  changed_POIs <- POIs %>%
    inner_join(select(st_drop_geometry(merged_POIs), COMID, oldPOI)) %>%
    select(COMID, oldPOI, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, 
           Type_Term, Type_Elev, Type_Travel) %>%
    mutate(to_com = COMID) %>%
    st_compatibalize(., read_sf(out_gpkg, poi_moved_layer))
  st_write(changed_POIs, out_gpkg, poi_moved_layer, append = TRUE) # write_sf not appending?
  
  # Format POI moves to table
  xWalk <- select(st_drop_geometry(segsub_dir), Move, COMID) %>%
    filter(!is.na(Move), Move %in% merged_POIs$COMID) %>%
    rename(COMID = Move, oldPOI = COMID)
  st_write(xWalk, out_gpkg, "collapse_xWalk", append = TRUE) # write_sf not appending?

  newSegs <- segment_creation(mutate(Seg2, old_POI_ID = POI_ID), xWalk)

  newSegs$raw_segs <- distinct(newSegs$raw_segs, .keep_all = T)
  
  # Return POIs, segments, and raw segments
  return (list(allPOIs = fin_POIs, segs = newSegs$diss_segs, FL = newSegs$raw_segs))
}

#' Identifies and attributes structural POIs
#'  @param nhdDF (sf data.frame) valid data frame of NHD flowlines
#'  @param final_POIs  (sf data.frame) final POIs
#' 
#' @return (data.frame columns) 1/0 Columns indicating structural POIs and the structural network
structPOIsNet <- function(nhdDF, final_POIs){
  
  # subset to flowlines used for nsegments
  inc_segs <- filter(nhdDF, !is.na(POI_ID))
  vpu <- inc_segs$VPUID
  
  # Terminal outlets from the initial network
  termout <- filter(nhdDF, !Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% final_POIs$COMID)
  
  # POIs that are also terminal outlets
  out_POIs <- filter(nhdDF, COMID %in% final_POIs$COMID & TerminalFl == 1)
  
  # Confluence POIs
  conf_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_Conf)])
  
  # Waterbody outlet POIs
  wb_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[!is.na(final_POIs$Type_WBOut) | !is.na(final_POIs$Type_WBIn)])
  
  # Waterbody inlets POIs
  struc_flines <- termout %>% 
    bind_rows(out_POIs, conf_POIs, wb_POIs) %>% 
    arrange(COMID)
  
  struc_flines <- struc_flines[!st_is_empty(struc_flines), , drop = F] %>%
    st_cast("LINESTRING")
  
  struc_POIs <- get_node(struc_flines, position = "end") %>% 
    mutate(COMID = struc_flines$COMID, LevelPathI = struc_flines$LevelPathI) 
  
  # Add back in any levelpaths missing (shouldn't be much)
  miss_LP <- filter(inc_segs, COMID %in% final_POIs$COMID) %>% 
    filter(!LevelPathI %in% struc_POIs$LevelPathI)
  
  # Only bind if there are rows present
  if (nrow(miss_LP) > 0){
    lp_POIs <- c(miss_LP$LevelPathI, struc_POIs$LevelPathI)
  } else {
    lp_POIs <- struc_POIs$LevelPathI
  }
  
  # final stuctural POIs, order by COMID to match with strucPOIs
  final_struct <- filter(nhdDF, Hydroseq %in% lp_POIs) %>% 
    arrange(COMID)
  struct_POIs <- get_node(final_struct, position = "end") %>% 
    mutate(COMID = final_struct$COMID) 
  
  #  produce unique set of flowlines linking POIs
  final_net <- unique(unlist(lapply(unique(final_struct$COMID), NetworkNav, st_drop_geometry(nhdDF))))
  
  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
  structnet <- filter(nhdDF, COMID %in% final_net & grepl(paste0("^", vpu, ".*"), .data$VPUID)) 
  
  # Write out minimally-sufficient network and POIs as a list
  return(list(struc_POIs = struc_POIs, structnet = structnet))
}


#'  Splits a network based elevation change and size
#'  @param in_POI_ID (integer) POI_ID of aggregated flowline that needs to be
#'                   split based on elevation
#'  @param split_DF  (sf data.frame) flowlines attributed with POI_ID
#'  @param elev_diff (numeric) max elevation change threshold within a segment
#'  @param max_DA (numeric) minimum drainage area for resulting split
#' 
#' @return (sf data.frame) flowlines with new POI_IDs identified (elev_POI_ID)
split_elev <- function(in_POI_ID, split_DF, elev_diff, max_DA){
  
  # subset to a given POI_ID
  split_elev_DF <- filter(split_DF, POI_ID == in_POI_ID) %>%
    filter(AreaSqKM > 0) %>%
    mutate(inc_seg_area = c(TotDASqKM[1], (TotDASqKM - lag(TotDASqKM))[-1]))
  
  
  first_iter <- unique(split_elev_DF$iter)
  # Iterate through segs that need splitting
  for (i in c(first_iter:(max(split_DF$iter, na.rm = T) + 10))){
    #print(i)
    
    # first split
    elev_pois_init_iter <- split_elev_DF %>%
      # filter by iteration, comes into play when multiple splits 
      #        are required per single POI_ID
      filter(iter == i) %>%
      # Recalculate inc elev, length, and area based on the last split
      mutate(csum_length = cumsum(LENGTHKM),
             inc_elev = cumsum(MAXELEVSMO - MINELEVSMO),
             sum_area = cumsum(inc_seg_area)) %>%
      # Determine if split actually necessary on this iteration
      filter((max(MAXELEVSMO) - min(MINELEVSMO)) > (elev_diff), sum_area > max_DA) %>%
      # This is an important step, identify which row to split on based on incremental
      #      elevaton value.  Some flowlines exceed the elevation difference threshold
      # The 'split_index' values determiens which row to create a new POI on
      mutate(split_index = ifelse(nrow(filter(., inc_elev > (elev_diff))) == nrow(.), 
                                  1, nrow(filter(., inc_elev < (elev_diff))))) %>%
      # Take out if at last iteration to avoid duplicating existng POIs
      filter(!COMID == POI_ID)
    
    if(nrow(elev_pois_init_iter) == 0){
      #print("done")
      # Done iterating on possible splits
      return(split_elev_DF)}
    
    # Get the flowline POI
    elev_pois_init_pre <- elev_pois_init_iter %>%
      filter(row_number() == split_index) %>%
      mutate(new_POI_ID = COMID) %>%
      select(COMID, POI_ID, new_POI_ID, Hydroseq_cut = Hydroseq) 
    
    # Remaining rows downstream after split
    elev_pois_init_post <- elev_pois_init_iter %>%
      left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
      filter(Hydroseq < Hydroseq_cut) %>%
      ungroup() 
    
    # Test for if we need to split again
    if(nrow(elev_pois_init_post) > 0){
      #print("yope")
      # Re-cacl the iter based on results above
      split_elev_DF <- split_elev_DF %>%
        left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
        mutate(elev_POI_ID = ifelse((is.na(elev_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, elev_POI_ID),
               iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), iter + 1, orig_iter)) %>%
        select(-c(new_POI_ID, Hydroseq_cut)) %>%
        ungroup()
      
    } else {
      split_elev_DF <- split_elev_DF %>%
        left_join(select(elev_pois_init_pre, -COMID), by = "POI_ID") %>%
        mutate(elev_POI_ID = ifelse((is.na(elev_POI_ID) & Hydroseq >= Hydroseq_cut), new_POI_ID, elev_POI_ID),
               iter = ifelse((!is.na(new_POI_ID) & Hydroseq < Hydroseq_cut), 0, orig_iter)) %>%
        select(-c(new_POI_ID, Hydroseq_cut)) %>%
        ungroup()
      
    }
  }
  return(split_elev_DF)
}

#'  Retrieves waterbodies from NHDArea layer or NHD Hi-Res for ResOpsUs 
#'  locations if absent from NHD waterbody
#'  @param nhd (sf data.frame) Fowlines for given VPU
#'  @param WBs  (sf data.frame) waterbodiess for discretization within VPU
#'  @param data_paths (list) data paths to data layers
#' 
#'  @return (list) wbs - sf data frame for NHDArea and HR waterbodies
#'                 wb_table - table of flowlines and outlet info for each 
#'                            feature in wb
HR_Area_coms <- function(nhd, WBs, data_paths, crs){
  
  # Pull out rows for VPU that are NHDArea
  nhd_area_resops <- WBs %>%
    filter(resops_flowlcomid %in% nhd$COMID, source == "NHDAREA")
  
  # Pull out rows for VPU that are NHD HR
  nhd_hr_wb_resops <-  WBs %>%
    filter(member_comid != 65000300139130) %>% # R09, in Canada
    filter(resops_flowlcomid %in% nhd$COMID, source == "HR ID AVAILABLE") %>%
    st_drop_geometry()
  
  # Get reachcodes for waterbody outlets, so we have an idea of which
  #     NHD HR 4-digit geodatabase we may need to retrieve
  RC <- filter(nhd, COMID %in% nhd_hr_wb_resops$resops_flowlcomid)$REACHCODE
  
  # If no NHDArea or HR waterbodies needed return NULL
  if (nrow(nhd_area_resops) == 0 & nrow(nhd_hr_wb_resops) == 0){
    return(NA)
  }
  
  # If NHDArea feature needed retrieve from National GDB
  if (nrow(nhd_area_resops) > 0){
    nhd_area_vpu <- read_sf(data_paths$nhdplus_gdb, "NHDArea") %>%
      filter(COMID %in% nhd_area_resops$member_comid) %>%
      mutate(source = "NHDv2Area")

    wb <- st_transform(nhd_area_vpu, crs)
  }
  
  # If NHDHR feature needed
  if (nrow(nhd_hr_wb_resops) > 0){
    # HUC04 we need to download
    huc04 <- substr(RC, 1, 4)

    # Download NHD HR HUC04 if we dont' have it, other wise load and
    # Bind NHDHR waterbodies into one layer
    hr_wb <- do.call("rbind", lapply(unique(huc04), function(x){
      print(x)
      vpu <- unique(substr(x, 1, 2))

      if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
                                paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")))){
        download_nhdplushr(data_paths$nhdplus_dir, unique(huc04))
      }
      
      gdb <- list.files(file.path(data_paths$nhdplus_dir, substr(vpu, 1, 2)), 
                        pattern = paste0("_", x, "_.+gdb"), full.names = TRUE)
    
      # Format to similar to NHDArea/Waterbody layers
      read_sf(gdb, "NHDWaterbody")})) #%>%
    
    # convert to lower case and designate new shape field
    names(hr_wb) <- tolower(names(hr_wb)) 
    st_geometry(hr_wb) <- "shape"
    
    hr_wb <- filter(hr_wb, permanent_identifier %in% WBs$member_comid) %>%
      rename(GNIS_NAME = gnis_name, GNIS_ID = gnis_id, 
             Permanent_Identifier = permanent_identifier, AREASQKM = areasqkm, 
             FTYPE = ftype, FCODE = fcode, FDATE = fdate, 
             REACHCODE = reachcode) %>%
      #select(-c(permanent_identifier, visibilityfilter, vpuid)) %>%
      st_zm(.) %>%
      st_as_sf() %>%
      mutate(source = "NHDHR",
             wb_id = ifelse(!is.na(GNIS_ID), GNIS_ID, Permanent_Identifier)) %>%
      st_transform(crs)

      # Bind or create new object
      if(exists("wb")){
        hr_wb <- st_compatibalize(hr_wb, wb)

        wb <- data.table::rbindlist(list(wb, hr_wb), fill = TRUE) %>%
          st_as_sf()
      } else {
        wb <- hr_wb
      }

  }
    
  # # get the outlt rows from the table
  # resops_outlet <- read.csv("data/reservoir_data/ISTARF-CONUS.csv") %>%
  #    filter(WBAREACOMI %in% wb$COMID | HR_NHDPLUSID %in% wb$COMID)
  
  # Clear out columns not needed
  WBs_sub <- WBs %>%
    select(-any_of(c("COMID", "GNIS_ID")))
    
  WBs_fin <- st_drop_geometry(WBs_sub) %>%
    mutate(member_comid = member_comid) %>%
    inner_join(select(wb, Permanent_Identifier, GNIS_ID, GNIS_NAME2 = GNIS_NAME), 
               by = c("member_comid" = "Permanent_Identifier")) %>%
    mutate(wb_id = ifelse(GNIS_ID %in% c(NA, " "), member_comid, GNIS_ID), 
           GNIS_NAME = ifelse(is.na(GNIS_NAME), GNIS_NAME2, GNIS_NAME)) %>%
    select(-c(GNIS_ID, GNIS_NAME2)) %>%
    st_as_sf() 
  
  # Create table of all flowlines that intersect the waterbody
  nhd_wb <- st_intersects(st_transform(nhd, st_crs(wb)), wb) 
  comid <- nhd[lengths(nhd_wb) > 0,]$COMID
  nhd_wb_all <- nhd_wb[lengths(nhd_wb) > 0] %>%
    purrr::set_names(comid) %>%
    stack() %>%
    # Ind is the default name for the set_names
    rename(comid = ind, nhd_wb = values) %>%
    mutate(wb_comid = wb[nhd_wb,]$Permanent_Identifier,
           outlet = ifelse(comid %in% WBs_fin$resops_flowlcomid, "outlet", NA),
           comid = as.integer(as.character(comid))) %>%
    #left_join(select(wb_table, DAM_ID, DAM_NAME, resops_FlowLcomid), 
    #          by = c("comid" = "resops_FlowLcomid")) %>%
    left_join(select(st_drop_geometry(wb), wb_id, Permanent_Identifier, GNIS_ID, 
                     source), 
              by = c("wb_comid" = "Permanent_Identifier"))
  
  return(list(nhd_wb_table = nhd_wb_all, wb = WBs_fin))
}

#'  Creates wb inlet and outlet events for splitting in hyRefactor
#'          for waterbodies derived from NHDArea and NHDHR waterbodies
#'  @param WBs  (sf data.frame) return from HR_Area_coms
#'  @param nhd_wb (sf data.frame) flowlines that intersect waterbodies
#'  @param type (character) whether to derive inlet or outlet points
#' 
#'  @return (sf data.frame) dataframe of WB inlet and outlet points to split
#'  
WB_event <- function(WBs, nhd_wb, type){
  # split into features and table
  WBs_table <- WBs$nhd_wb_table
  WBs_layer <- WBs$wb

  if (type == "outlet"){
    # get the outlet comid from the ResOps Table
    outlet_fl <- filter(nhd_wb, COMID %in% filter(WBs_table, outlet == "outlet")$comid)
    
    # Get the downstream flowline for continuity
    ds_fl <- filter(nhd_wb, Hydroseq %in% outlet_fl$DnHydroseq,
                    LevelPathI %in% outlet_fl$LevelPathI) %>%
      rbind(outlet_fl) %>%
      arrange(desc(Hydroseq)) %>%
      group_by(LevelPathI) %>%
      # union together 
      summarize(do_union = TRUE) %>%
      st_cast("LINESTRING")
    
    WBs_HR <- filter(WBs_layer, source == "HR ID AVAILABLE")
    
    # For NHD HR waterbody outlets its a bit more complicated
    if (nrow(WBs_HR) > 0){
      # Get the flowlines intersecting the HR waterbody and find one with the
      #     max DA
      outlet_wb_int <- nhd_wb[lengths(st_intersects(st_transform(nhd_wb, st_crs(WBs_HR)),
                                                                 WBs_HR)) > 0,] %>%
        group_by(wb_id) %>%
        filter(TotDASqKM == max(TotDASqKM)) %>%
        ungroup()

      # get the ds flo with the same levepath (JIC)
      ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_wb_int$Hydroseq,
                      LevelPathI %in% outlet_wb_int$LevelPathI)

      outlet_fl <- rbind(outlet_wb_int, ds_fl)

      # Cast flowlines within NHDHR waterbody to point
      WB_FL_pnts <- outlet_wb_int %>%
        st_cast("POINT") %>%
        group_by(wb_id) %>%
        mutate(pnt_id = row_number()) %>%
        ungroup()

      outlet_pnts <- do.call("rbind", lapply(unique(WB_FL_pnts$wb_id), function(x){
        fl <- filter(WB_FL_pnts, wb_id == x)
        wb <- filter(WBs_HR, wb_id == x)
        
        # Determine which points intersect waterbody
        WB_outlet_pnts <- fl[lengths(st_intersects(fl, wb)) > 0,] %>%
          st_drop_geometry() %>%
          group_by(wb_id) %>%
          mutate(within_wb_id = row_number()) %>%
          filter(within_wb_id >= max(within_wb_id)) %>%
          ungroup() %>%
          select(wb_id, orig_pnt_id = pnt_id, within_wb_id)
        
        # Deriv new linestring by concating points from most upstream point
        #       within waterbody to downstream so we can split at FL/waterbody
        #       nexus
        outlet_FL <- fl %>%
          inner_join(WB_outlet_pnts, by = "wb_id") %>%
          select(wb_id, pnt_id, orig_pnt_id, within_wb_id) %>%
          filter(pnt_id >= orig_pnt_id) %>%
          group_by(wb_id) %>%
          summarize(do_union = F) %>%
          st_cast("LINESTRING") %>%
          filter(wb_id %in% wb$wb_id)
        
        outlet_pnt <- sf::st_intersection(outlet_FL, WBs_HR) %>%
          st_cast("MULTIPOINT") %>%
          st_cast("POINT") %>% # was point
          group_by(wb_id) %>%
          filter(row_number() == max(row_number(), na.rm = T)) %>%
          ungroup() %>%
          mutate(id = row_number()) %>%
          filter(wb_id == wb_id.1)
      }))
      
      # Derive the events
      if(exists("wb_events")){
        hr_events <- get_flowline_index(st_transform(nhd_wb, st_crs(outlet_pnts)),
                                                     outlet_pnts) %>%
          inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id), by = "COMID") %>%
          filter(wb_id %in% WBs_HR$wb_id) %>%
          mutate(event_type = type) %>%
          cbind(select(outlet_pnts, geom)) %>%
          st_as_sf() %>%
          st_transform(st_crs(wb_events)) %>%
          st_compatibalize(., wb_events)

        wb_events <- rbind(wb_events, hr_events) %>%
          distinct()
      } else {
        wb_events <- do.call("rbind", lapply(unique(outlet_pnts$wb_id), function(x){
          outlet_pnt <- filter(outlet_pnts, wb_id == x)
          get_flowline_index(nhd_wb, outlet_pnt)})) %>%
          left_join(distinct(st_drop_geometry(outlet_pnts), resops_flowlcomid, wb_id), 
                    by = c("COMID" = "resops_flowlcomid")) %>%
          mutate(event_type = type) %>%
          cbind(select(outlet_pnts, geom)) %>%
          st_as_sf()
      }

    }
    
  # For inlet points its alot easier for both NHDARea and NHDHR
  } else {
    start_pts <- get_node(nhd_wb, position = "start") %>%
      cbind(st_drop_geometry(nhd_wb))
    
    inlet_FL <- nhd_wb[lengths(st_intersects(
      st_transform(start_pts, st_crs(WBs_layer)), WBs_layer)) == 0,] %>%
      rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq, 
                   LevelPathI %in% .$LevelPathI)) %>%
      arrange(desc(Hydroseq)) %>%
      unique()
    
    inlet_ls <- inlet_FL %>%
      group_by(LevelPathI) %>%
      st_cast("POINT") %>%
      summarize(do_union = F) %>%
      st_cast("LINESTRING") %>%
      ungroup()
    
    inlet_pnts_int <- sf::st_intersection(
      st_transform(inlet_ls, st_crs(WBs_layer)), WBs_layer) #%>%
      #st_cast("LINESTRING")
      
    int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
    if(nrow(int_point) > 0){
      inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
        st_cast("LINESTRING")
      inlet_pnts_int <- rbind(int_point, inlet_fl)
    }
    
    inlet_pnts <- inlet_pnts_int %>%
      group_by(LevelPathI) %>%
      st_cast("POINT") %>%
      mutate(id = row_number()) %>%
      filter(id == min(id), !is.na(wb_id)) %>%
      ungroup() 

    wb_events <- get_flowline_index(st_transform(nhd_wb, st_crs(inlet_pnts)),
                                    inlet_pnts) %>%
      inner_join(select(st_drop_geometry(nhd_wb), COMID, wb_id, LevelPathI), by = "COMID") %>%
      mutate(event_type = type) %>%
      inner_join(select(inlet_pnts, LevelPathI), by = "LevelPathI") %>%
      #group_by(COMID, LEVELPATHI) %>%
      mutate(nexus = T) %>%
      rename(POI_identifier = wb_id) %>%
      select(-c(id, offset)) %>%
      st_as_sf()
  }
  return(wb_events)
} 


#'  Creates POIs for gages using refactor criteria
#'  @param temp_POIs  (sf data.frame) current data frame of POIs
#'  @param gages_info (sf data.frame) VPU specific streamgages from 01_gage_selection
#'  @param nhd (sf data.frame) flowline data.frame 
#'  @param combine_meters (integer) refactor threshold (m) for if if two adjacent fl should be combined
#'  @reach_meas_thresh (integer) threshold added or substracted from reach_meas to determine if a gage will split fl
#' 
#'  @return (sf data.frame) dataframe of gage POIs
#'  
gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_meas_thresh){
  # Create streamgage POIs
  gage_POIs <- POI_creation(select(st_drop_geometry(gages_info), COMID, site_no), nhd, "Gages") %>%
    st_compatibalize(., tmp_POIs)
  
  # Avoid these for refactoring
  avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)
  
  # Create events for streamgages for flowline splitting
  events <- gages_info %>%
    # bring over NHD information
    inner_join(select(st_drop_geometry(gage_POIs), Type_Gages), by = c("site_no" = "Type_Gages")) %>%
    # Get over the NHD attributes
    #inner_join(select(st_drop_geometry(nhd), AreaSqKM, REACHCODE, LENGTHKM, COMID, FromMeas, ToMeas), 
    #           by = "COMID") %>%
    # Apply reach measure thresholds
    filter(reach_meas - FromMeas > reach_meas_thresh & AreaSqKM > 2 & 
             ToMeas - reach_meas > reach_meas_thresh & LENGTHKM > (combine_meters / m_per_km)) %>%
    select(COMID, REACHCODE = reachcode, REACH_meas = reach_meas, Type_Gages = site_no) %>%
    filter(!COMID %in% avoid$COMID) %>%
    mutate(event_type = "streamgage") %>%
    st_compatibalize(., tmp_POIs) %>%
    mutate(nexus = TRUE)
  
  # Reset gage flag if even is created
  gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
    addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE) 
  
  if(nrow(events) > 0){
    tmp_POIs <- data.table::rbindlist(list(gage_POIs_nonevent, 
                                           select(events, COMID, Type_Gages, nexus)), fill = TRUE) %>%
      mutate(nexus = ifelse(is.na(nexus), FALSE, nexus)) %>%
      st_as_sf()
  } else {
    events <- NA
    tmp_POIs <- mutate(tmp_POIs, nexus = FALSE)
  }

  
  return(list(events = events, tmp_POIs = tmp_POIs))
}

#'  Creates Waterbody POIs, calls a few other functions
#'  @param nhd (sf data.frame) flowline data.frame 
#'  @param WBs_VPU (sf data.frame) 
#'  @param data_paths (list) data paths to data layers
#'  @param crs (integer) CRS to use (epsg code) 
#' 
#'  @return (sf data.frame) dataframe of waterbody outlet POIs
#'  
wbout_POI_creaton <- function(nhd, wb_table, data_paths, tmp_POIs, crs){
  
  wb_components_table <- st_drop_geometry(wb_table) %>%
    dplyr::mutate(member_comid = strsplit(member_comid, ",")) %>%
    tidyr::unnest(cols = member_comid) %>%
    mutate(member_comid = as.integer(member_comid)) %>%
    distinct() 
  
  # R17, some wierd stuff going on in the waterbodies
  if("947070203" %in% wb_components_table$member_comid){
    wb_components_table <- wb_components_table %>%
      filter(!wb_id %in% c(1503847, 1513298))
  }

  # Bring wb_id to NHD data frame
  if(("wb_id") %in% colnames(nhd)){
    nhd <- select(nhd, -wb_id)
  }

  nhd <- nhd %>%
    left_join(distinct(st_drop_geometry(wb_components_table), wb_id, member_comid) %>%
                mutate(member_comid = as.numeric(member_comid)),
              by = c("WBAREACOMI" = "member_comid"))


  # Resops that are either NHDArea or MHD HR source and are missing in the
  #       nhd attributiong table
  resops_wb_other_sources <- filter(wb_table, source ==
                                      "HR ID AVAILABLE" | (source == "NHDAREA" &
                                      !wb_id %in% nhd$wb_id) &
                                      accounted == 0)

  # Resops that are not accounted for and NHDv2 waterbody or NHDArea sourced
  resops_main <- filter(wb_table, !is.na(resops_flowlcomid) &
                          !wb_id %in% resops_wb_other_sources$wb_id,
                        accounted == 0)

  # Remaining WB that don't have outlet COMID defined and aren't in ResOPsUS
  WB_outlet <- filter(wb_table, is.na(resops_flowlcomid),
                      !wb_id %in% tmp_POIs$Type_WBOut)

  # Size-based waterbody outlet POIs ot in NHDArea, HR, or
  #        manually defined with ResOPS
  wbout_COMIDs <- nhd %>%
    filter(wb_id %in% WB_outlet$wb_id) %>%
    group_by(wb_id) %>%
    slice(which.min(Hydroseq)) %>%
    ungroup() %>%
    switchDiv(., nhd) %>%
    select(COMID, wb_id) %>%
    mutate(wbtype = "size") %>%
    st_drop_geometry() %>%
    distinct()

  # ResOpsUS defined outlets
  resops_main_outlets <- select(st_drop_geometry(resops_main), COMID = resops_flowlcomid,
                                wb_id) %>%