NHD_navigate.R

# # NHDPlus V2 Region
# if(!exists("VPU") && !nchar(VPU <- Sys.getenv("HYDREG")) > 1) 
#   VPU <- "01"
# 
# # Source Data
# crs <- 5070
# data_paths <- data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
# huc12_xWalk <-  file.path(data_paths$nhdplus_dir, "CrosswalkTable_NHDplus_HU12.csv")
# HUC12 <- file.path(data_paths$nhdplus_dir, "HUC12.rds")
# 
# # Defined during NHD_Navigate.Rmd
# nav_gpkg <- file.path("cache", paste0("GF_", VPU,".gpkg"))
# nhdflowline <- "NHDFlowline" # NHD flowlines per VPU
# temp_POIs <- paste0("tmp_POIs_", VPU) # Running POI list, precedes "POIs_VPU"
# WBs <-  paste0("WB_", VPU) # Waterbodies within VPU
# pois_all <- paste0("POIs_", VPU) # All POIs binded together
# poi_moved <- paste0("POIs_mv_", VPU) # POIs moved from original COMID assignment
# poi_xWalk <- paste0("poi_xWalk_", VPU) # POIs that changed COMIDS during the navigate part of the workflow
# n_segments <- paste0("nsegment_", VPU) # Minimally-sufficient network dissolved by POI_ID
# pois_merge <- paste0("merPOIs_", VPU) # All POIs binded together
# 
# # Making an "events" script to address "events" information
# #     i.e. gages, NID, TE Plants, waterbodies
# gage_gpkg <- "cache/Gage_info.gpkg"
# gage_cand <- "Potential_Gages"
# VPU_gage_Table <- file.path("cache", paste0("R", VPU, "_Gages.csv"))
# CONUS_gage_Table <- "data/gages_MDA.rds"
# 
# # Defined during NonDend.Rmd
# fin_gpkg <- file.path("cache", paste0(VPU, "_final.gpkg"))
# lookup_miss <- paste0("lookup_missing_", VPU)
# missing_cats <- paste0("miss_cats_", VPU)
# protoHRUs <- paste0("poi_cats_", VPU)
# cat_borders <- paste0("perimeters_", VPU)


NetworkNav <- function(inCom, nhdDF, withTrib){
  ##########################################
  # Network navigation for upstream/downstream from a COMID of interest
  # 
  # Arguments:
  #   inCOM : (list) list of input COMIDs that are 'dangles'
  #   nhdDF : (data frame) valid data frame of NHD flowlines
  #   withTrib : (logical) flag for if the upstream navigation should include tributaries
  #              or stick to mainstem level path
  #
  # Returns:
  #   mydata : (list) list of COMIDs to connect dangle to network
  ##########################################
  #print (inCom)
  if (missing(withTrib)){
    seg <- nhdplusTools::get_UM(nhdDF, inCom, include = TRUE)
  } else {
    seg <- nhdplusTools::get_UT(nhdDF, inCom)
  }
  return(seg)
}


NetworkConnection <- function(incom, nhd){
  ##########################################
  # Connects dangles in the network that are not
  #          terminal outlets
  # 
  # Arguments:
  #   inCOM : (list) list of input COMIDs that are 'dangles'
  #   nhdDF : (data frame) valid data frame of NHD flowlines
  #
  # Returns:
  #   mydata : (list) list of COMIDs to connect dangle to network
  ##########################################
  # data frame of connections that need to be made
  upnet_DF <- filter(nhd, COMID %in% incom) %>%
    filter(!DnHydroseq %in% Hydroseq)
  
  # while the number of dangles is greater than 0
  while (length(upnet_DF$COMID) > 0){
    # create item for number of dangles
    count <- dim(upnet_DF)[1]
    print (dim(upnet_DF))
    # find out which level paths are downstream of dangling huc12 POIs
    DSLP <- upnet_DF %>% pull(DnLevelPat)#[upnet_DF$COMID %in% incom]
    # Get the COMID of the hydroseq with level path value
    # the lowest downstream flowline within the levelpath
    inCom2 <- nhd$COMID[nhd$Hydroseq %in% DSLP]
    # Run the upstream navigation code
    upNet <- unique(unlist(lapply(inCom2, NetworkNav, nhd)))
    # Append result to existing segment list
    incom <- append(incom, upNet)
    # Get the same variable as above
    upnet_DF <- filter(nhd, COMID %in% incom, !DnHydroseq %in% Hydroseq)
    # Get the count
    count2 <- dim(upnet_DF)[1]
    # if the count has remained the same we are done and return the flowline list
    if (count == count2){
      return (incom)
    }
  }
  # Not sure this other return is needed
  return(incom)
}


switchDiv <- function(insegs, nhd){
  ##########################################
  # Swith divergence path for POIs from minor (2) to major (1)
  # reduces the number of short, spurious segments and POIS as locations such as waterbody outlets
  # 
  # Arguments:
  #   inSegs : (list) list of input COMIDs that are 'dangles'
  #   nhd : (data frame) valid data frame of NHD flowlines
  #
  # Returns:
  #   mydata : (list) list of COMIDs for major diversions
  ##########################################
  div_segs <- filter(nhd, COMID %in% insegs$COMID)
  if (2 %in% div_segs$Divergence){
    print ("Switching divergence to other fork")
    # Look Upstream
    upstream <- st_drop_geometry(nhd) %>% 
      inner_join(st_drop_geometry(div_segs) %>% 
                   filter(Divergence == 2) %>% 
                   select(COMID, Hydroseq), 
                 by = c("DnMinorHyd" = "Hydroseq"))
    
    # From Upstream Segment switch POI to the downstream major/main path
    insegs_maj <- st_drop_geometry(nhd) %>% 
      inner_join(upstream %>% select(COMID.y, DnHydroseq), 
                 by = c("Hydroseq" = "DnHydroseq")) %>% 
      select(COMID, COMID.y)
    
    insegs <- insegs %>% 
      left_join(insegs_maj, by = c("COMID" = "COMID.y")) %>%
      mutate(COMID = ifelse(!is.na(COMID.y), COMID.y, COMID)) %>% select(-COMID.y)

  } else {
    print ("no divergences present in POI set")
  }
  return(insegs)
}

POI_creation<-function(srcData, nhd, IDfield){
  ##########################################
  # Creates POIs for a given data theme
  # 
  # Arguments:
  #   srcData : (data frame) DF with two columns: 
  #             1 - COMID
  #             2 - Unique ID value for POI (Streamgage ID, etc.)
  #   nhd : (data frame) valid data frame of NHD flowlines
  #   IDfield : (character) Name of 'Type' field to be modified in POI 
  #
  # Returns:
  #   mydata : (simple features) POIs for the specific data theme
  ##########################################
  
  # Give generic ID to Identity field
  colnames(srcData) <- c("COMID", "ID")
  
  sub_segs <- filter(nhd, COMID %in% srcData$COMID) 
  #merge(st_drop_geometry(sub_segs %>% select(COMID))
  POIs <- sub_segs %>% 
    get_node(., position = "end") %>% 
    mutate(COMID = sub_segs$COMID) %>%
    mutate(Type_HUC12 = NA, Type_WBIn = NA, Type_WBOut = NA, Type_Gages = NA, Type_TE = NA, Type_NID = NA, Type_Conf = NA) %>%
    inner_join(srcData %>% select(COMID, ID), by = "COMID") %>% 
    mutate(!!(paste0("Type_", IDfield)) := ID) %>% 
    select(COMID, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, geometry)

  #if("geom" %in% colnames(Seg_ends)) Seg_ends <- Seg_ends %>% rename(geometry=geom)

  return(POIs)
}

addType<-function(new_POIs, POIs, IDfield, bind){
  ##########################################
  # Checks if existing POIs are co-located with new POI set
  #        Adds the type attribute for co-located POIs of multiple themes
  # 
  # Arguments:
  #   POIs : (data frame) Existing POIs
  #   newPOIs: (data frame) newly-derived POIs of a given data theme
  #   IDfield : (character) Name of 'Type' field to be modified in POI 
  #
  # Returns:
  #   mydata : (data frame) Existing POIs with Type field modified for
  #            duplicate POIs for given data theme
  ##########################################
  new_POIs <- st_compatibalize(new_POIs, POIs) 
  
  POIs_exist <- POIs %>% 
    rename(ID = !!(paste0("Type_", IDfield)))
  
  POIs_newAtt <- filter(new_POIs, COMID %in% POIs$COMID) %>%
    rename(ID2 = !!(paste0("Type_", IDfield)))

  POIs_fin <- POIs_exist %>% left_join(st_drop_geometry(POIs_newAtt) %>% 
                                 select(COMID, ID2), by = "COMID", all.x = TRUE) %>%
    mutate(ID = ifelse(!is.na(ID2), ID2, ID)) %>% 
    rename(!!(paste0("Type_", IDfield)) := ID) %>% 
    select(-ID2) 
 
  # Bind unless indicated not to
  if(missing(bind)){
    POIs_fin <- rbind(POIs_fin, filter(new_POIs, !COMID %in% POIs_newAtt$COMID))
  }
  
  return(POIs_fin)
}

segment_increment <- function(nhd, POIs){
  ##########################################
  # Creates raw and dissolved segment layers with necessaary
  #         upstream/downstream routing attribution
  # 
  # Arguments:
  #   nhd : (data frame) minimally-sufficient flowlines that connect POIs
  #          and meet requirements of NHDPlusTools
  #   POIs : (data frame) Raw POI data frame
  #   Field : (character) ID field of unique Segment identifier
  #  
  # Returns:
  #   segList : (data frame) raw segments
  #             
  ##########################################
  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
  nhd <- mutate(nhd, 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(nhd))
  #POI_list <- lapply(c(1:nrow(seg_POIs)), POI_ID_assign, seg_POIs, st_drop_geometry(nhd))
  stopCluster(clust)
  
  inc_segs <- data.table::rbindlist(POI_list)
  
  print(proc.time()-ptm)
  return(inc_segs)
}


segment_creation <- function(nhd, routing_fix){
  ##########################################
  # Creates raw and dissolved segment layers with necessaary
  #         upstream/downstream routing attribution
  # 
  # Arguments:
  #   inSegs : (data frame) segments with incremental IDs (POI_ID populated)
  #   nhd : (data frame) full nhd data frame
  #   routing_fix : (logical) routing fixes if available
  #  
  # Returns:
  #   segList : (data frame) raw segments
  #             
  ##########################################
  
  if(!"StartFlag" %in% names(nhd)) {
    nhd$StartFlag <- ifelse(nhd$COMID %in% nhd$tocomid, 0, 1)
  }
  
  in_segs <- filter(nhd, !is.na(POI_ID))
  
  routing_fix <- routing_fix %>%
    rename(COMID = oldPOI, new_COMID = COMID)
  
  # If there are routing fixes to account for if a POI with a DA of 0 is moved upsream or downstream
  if (missing(routing_fix) == FALSE){
    # Above we generated the network using the initial set of POIs; here we crosswalk over the old COMIDs to the new
    nhd_fix <- nhd %>%
      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(nhd), !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))
}

DS_poiFix <- function(POIs_wgeom, nhd){
  ##########################################
  # Moves POI Upstream or downstream if it falls on COMID
  #       of flowline with no corresponding catchment
  #
  # Arguments:
  #   POIs : (data frame) POI data set
  #   nhd : (data frame) valid data frame of NHD flowlines
  #
  # Returns:
  #   repPOIs_unique : (data frame) DF of POIs with new COMID associated
  ##########################################
  POIs <- st_drop_geometry(POIs_wgeom) %>%
    arrange(COMID)

  # DF of downstream segment
  tocomDF <- select(nhd, COMID, DnHydroseq) %>%
    inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq"))

  # Find segments with POIs where there is no corresponding catchment
  unCon_POIs <- filter(nhd, COMID %in% POIs$COMID, AreaSqKM == 0)

  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"))
  
  # 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_fix_unique <- filter(POIs, COMID %in% poi_fix$COMID) %>% 
    bind_rows(poi_fix) %>% 
    group_by(COMID) %>% 
    filter(n() > 1) %>%
    summarise_each(funs(min(unique(.[!is.na(.)]), na.rm = T))) %>%
    mutate_all(~na_if(., -Inf))

  # Join new POI COMIDs and geometry with the old Type fields
  new_POIs <- bind_rows(filter(poi_fix, !COMID %in% poi_fix_unique$COMID), poi_fix_unique) %>%
    arrange(COMID) %>%
    bind_cols(get_node(filter(nhd, COMID %in% .$COMID) %>% arrange(COMID), position = "end")) %>%
    mutate(to_com = NA) %>%
    st_sf() %>%
    st_compatibalize(., POIs_wgeom)

  # If the DS replacement flowlines with Inc Area > 0 equals the number of replacement comIDs 
  replacement_tocoms <- st_drop_geometry(new_POIs) %>%
    inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq), by = c("DnHydroseq" = "Hydroseq")) %>%
    mutate(to_com = COMID.y) %>%
    select(oldPOI, COMID.x, COMID.y, AreaSqKM)
  
  # Add downstream COMIDS to the newCOMID
  fin_POIs <- new_POIs %>% 
    left_join(replacement_tocoms, by = "oldPOI") %>% 
    mutate(to_com = ifelse(!is.na(COMID.y), COMID.y, NA)) %>% 
    select(-c(TotDASqKM, AreaSqKM.x, DnHydroseq, COMID.x, COMID.y, AreaSqKM.y)) 
  
  return (list(xWalk = select(poi_fix, COMID, oldPOI), new_POIs = fin_POIs))
}


movePOI_NA_DA <- function(poi_fix, nhd){
  ##########################################
  # 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
  # 
  # Arguments:
  #   poi_fix : (data frame) POI data set of COMIDs to be changed
  #   nhd : (data frame) valid data frame of NHD flowlines
  #
  # Returns:
  #   newPOI : (data frame) DF of POIs with new COMID associated
  ##########################################
  
  # Closest POI/US/DS
  up_segs <- nhdplusTools::get_UM(nhd, poi_fix, sort=T)
  seg2fix <- filter(nhd, COMID == poi_fix)
  
  # Sorted results and filter out all flowlines w/o catchments
  upstuff <- filter(nhd, COMID %in% unlist(up_segs)) %>% 
    arrange(factor(COMID, levels = up_segs)) %>%
    filter(AreaSqKM > 0)
  
  down_segs <- nhdplusTools::get_DM(nhd, poi_fix, sort=T)
  downstuff <- filter(nhd, COMID %in% unlist(down_segs)) %>% 
    arrange(factor(COMID, levels = down_segs)) %>%
    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)
}

#out_gages <- POI_move_down(nav_gpkg, final_POIs, nsegments_fin, filter(nhd_Final, !is.na(POI_ID)), "Type_Gages", .05)
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
  ##########################################
  # Moves down/up POIs based on data theme
  # 
  # Arguments:
  #   out_gpkg : (geopackage) Output geopackage
  #   POIs: (data frame) Original  POIs
  #   Seg : (data frame) dissolved segments
  #   Seg2 : (data frame) raw segments
  #   exp : (string) Type of thematic POI being moved round
  #   DA_diff : (float) threshold of drainage area difference to
  #             consider when merging two POIs
  #
  # Returns:
  #   list : 1 - New set of POIs
  #          2/3 - correpsonding segments (both raw and dissolved)
  ##########################################
  
  POIs <- POIs %>% mutate_if(is.numeric, function(x) ifelse(is.infinite(x), NA, x))
  
  # 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, DAR, IncDA)

  # downstream 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, DAR, IncDA)
  
  # compiled list of POIs to move up or down
  seg_choices <- 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)

  # Filter by POI Type
  Types <- c("Type_Gages", "Type_TE", "Type_Conf", "Type_NID", "Type_WBIn", "Type_WBOut")
  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)))
  
  seg_sub <- seg_choices %>% 
    filter_at(Types, all_vars(is.na(.))) %>%
    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, COMID, Move) %>% 
    inner_join(st_drop_geometry(.) %>% 
                 select(COMID, Move), 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")) %>%
    # Bring over relevant attributes
    mutate(Type_HUC12 = ifelse(is.na(Type_HUC12.y), Type_HUC12.x, Type_HUC12.y)) %>%
    mutate(Type_Gages_A = ifelse(is.na(Type_Gages.x) & !is.na(Type_Gages.y), Type_Gages.y, Type_Gages.x)) %>%
    # Gages_B is incase there are two gages being merged together, not writing out for now
    mutate(Type_Gages_B = ifelse(is.na(Type_Gages.y) & !is.na(Type_Gages.x), Type_Gages.y, NA)) %>%
    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(oldPOI = COMID.y, to_com = NA) %>% 
    st_sf()
  
  # Bring new POI data over to new data
  fin_POIs <- filter(POIs, !COMID %in% merged_POIs$oldPOI) %>% 
    left_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, to_com)) %>%
    select(COMID, oldPOI, Type_HUC12, Type_Gages, Type_TE, Type_NID, Type_WBIn, Type_WBOut, Type_Conf, to_com) %>%
    st_compatibalize(., read_sf(out_gpkg, poi_moved))
  st_write(changed_POIs, out_gpkg, poi_moved, append = TRUE) # write_sf not appending?
  
  # Format POI moves to table
  xWalk <- select(st_drop_geometry(segsub_dir), Move, COMID) %>%
    filter(!is.na(Move)) %>%
    rename(COMID = Move, oldPOI = COMID)
  st_write(xWalk, out_gpkg, poi_xWalk, append = TRUE) # write_sf not appending?

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

  # Return POIs, segments, and raw segments
  return (list(allPOIs = fin_POIs, segs = newSegs$diss_segs, FL = newSegs$raw_segs))
}

writePOIs <- function(POIs, out_gpkg, Type){
  ##########################################
  # Write out final POI datasets with information
  # 
  # Arguments:
  #   POIs : (data frame) POI data set 
  #   out_gkpg : (geopackage) Geopackage where final POI layers written
  #   Type : (character) Type of POI being written; default is write features for all types
  #
  # Returns:
  #   finPOIs : (data frame) DF of final POIs 
  ##########################################
  
  print ("Writing out final POIs")
  # If type is missing write out all flowlines
  if (missing(Type)){
    print ("Writing out all POI Types")
    lyrs <- st_layers(out_gpkg)
    # get subcategory of POIs
    POI_names <- lyrs$name[!is.na(lyrs$geomtype) & lyrs$geomtype== "Point"]
    print (POI_names)
    for (poi in POI_names){
      print (poi)
      poi_type <- unlist(strsplit(poi, "_"))[1]
      sub_POIs <- POIs %>%
        filter(!is.na(!!as.name(paste0("Type_", poi_type))))
      
      write_sf(sub_POIs, out_gpkg, poi)
    }
  } else {
    sub_POIs <- POIs %>%
      filter(!is.na(!!as.name(paste0("Type_", Type))))
    write_sf(sub_POIs, out_gpkg, poi)
  }
}


structPOIsNet <- function(nhd, final_POIs){
  ##########################################
  # Produce the structural POIs
  # 
  # Arguments:
  #   ncombined : (data frame) final Segments 
  #   nhd : (data frame) valid data frame of NHD flowlines
  #   finalPOIs :  (data frame) final POIs
  #   out_gkpg : (geopackage) Geopackage where final POI layers written
  #
  # Returns:
  #   writes Structural POIs and segments to geopackage 
  ##########################################
  
  # subset to flowlines used for nsegments
  inc_segs <- filter(nhd, !is.na(POI_ID))
  
  # Terminal outlets from the initial network
  termout <- filter(nhd, !Hydroseq %in% DnHydroseq & TerminalFl == 1 & !COMID %in% final_POIs$COMID)
  
  # POIs that are also terminal outlets
  out_POIs <- filter(nhd, COMID %in% final_POIs$COMID & TerminalFl == 1)
  
  # Confluence POIs
  conf_POIs <- filter(inc_segs, COMID %in% final_POIs$COMID[final_POIs$Type_Conf == 1])
  
  # 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_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(nhd, 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(nhd))))
  
  # Subset NHDPlusV2 flowlines to navigation results and write to shapefile
  structnet <- filter(nhd, 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))
}

outlets_close <- function(nhd, cats){
  ##########################################
  # Find little outlets and determine closest valid catchment with HRU ID
  # 
  # Arguments:
  #   nhd : (data frame) valid data frame of NHD flowlines
  #   cats : (data frame) valid data frame of NHD catchments
  #
  # Returns:
  #   missing_outlets : (data frame) DF of outlets and closest proto-HRU
  ##########################################
  
  # Networks that terminate and are smaller than the threshold, previously defined 
  little_outlets <- filter(nhd, TerminalFl == 1,
                           TotDASqKM <= min_da_km &
                             COMID %in% filter(cats, is.na(POI_ID))$FEATUREID) %>%
    arrange(COMID)
  
  # convert flowlines to end points
  xy <- get_node(little_outlets) %>%
    st_transform(26904)
  
  # All cats except those with same COMId as little_outlets
  nonout_cats <- filter(cats, !FEATUREID %in% little_outlets$COMID) %>%
    mutate(ID = row_number())
  
  # cast to multlinestirng sf feature for proximity analysis
  cats_lines <- filter(cats, FEATUREID %in% little_outlets$COMID) %>%
    group_by(FEATUREID) %>%
    summarize(do_union = FALSE) %>%
    st_cast("MULTILINESTRING") %>%
    arrange(FEATUREID)
  
  # find and measure distance to closest HRUs
  missing_outlets <- bind_cols(st_drop_geometry(little_outlets) %>% select(COMID), xy) %>% 
    st_sf() %>%
    # Distance from little outlet to catchment boundary
    mutate(dist2edge = round(as.numeric(st_distance(., cats_lines, by_element = TRUE)), 2)) %>%
    # Nearest catchment that is not its own
    mutate(nearest_cat = st_nearest_feature(., nonout_cats)) %>%
    # Distance to nearest HRU
    mutate(dist2_HRU = round(as.numeric(st_distance(., nonout_cats[nearest_cat,], by_element = TRUE)), 2)) %>% 
    #st_drop_geometry() %>%
    inner_join(select(st_drop_geometry(cats), FEATUREID, HUC_12), by = c("COMID" = "FEATUREID")) %>%
    inner_join(select(st_drop_geometry(nonout_cats), FEATUREID, HUC_12, POI_ID, ID), by = c("nearest_cat" = "ID")) %>%
    rename(HUC_12_outlet = HUC_12.x, FEATUREID_nrHRU = FEATUREID, HUC_12_nrHRU = HUC_12.y) %>%
    filter(!dist2_HRU > dist2edge)
  
  return(missing_outlets)
}

perims <- function(cats){
  ##########################################
  # Determines length of shared perimeters with adjacent hRUS
  # 
  # Arguments:
  #   cats : (data frame) valid data frame of NHD catchments
  #
  # Returns:
  #   perimeters : (data frame) DF of shared perimeter lengths for cats
  ##########################################
  
  # Perimeters and catchments that touch
  Touching_List <- st_touches(cats)
  perimeters <- lwgeom::st_perimeter(cats)
  
  # Get perimeter lengths and IDs of neightbors of catchments
  neighbors <- all.length.list <- lapply(1:length(Touching_List), function(from) {
    lines <- st_intersection(cats[from,], cats[Touching_List[[from]],])
    lines <- st_cast(lines) # In case of multiple geometries
    l_lines <- st_length(lines)
    res <- data.frame(origin = from,
                      perimeter = as.vector(perimeters[from]),
                      touching = Touching_List[[from]],
                      t.length = as.vector(l_lines),
                      t.pc = as.vector(100*l_lines/perimeters[from]))
    res})
  
  perimeters <- do.call("rbind", neighbors)
  
  return(perimeters)
}


assignable_cats <- function(perimeters, cats, missing_outlets){
  ##########################################
  # Assigns cats with no POI_ID a POI_ID based on proxmity or 
  #         other spatial relationships where determined
  # 
  # Arguments:
  #   perimeters : (data frame) shared perimeters data frame 
  #   cats : (data frame) valid data frame of NHD catchments
  #   missing outlets : (data frame) outlets of catchments w/o POI_ID
  #
  # Returns:
  #   cats : (data frame) DF of cats with populated POI_ID
  ##########################################
  
  # limit the result to cats that have not been assigned POI_ID
  all.length.df <- perimeters %>%
    inner_join(select(st_drop_geometry(cats), FEATUREID, POI_ID, HUC_12, ID), 
               by = c("origin" = "ID")) %>%
    left_join(select(st_drop_geometry(cats), FEATUREID, POI_ID, HUC_12, ID), by = c("touching" = "ID")) %>%
    rename(FEATUREID = FEATUREID.x, shared_length = t.length, HUC_12_noPOI = HUC_12.x, 
           FEATUREID_bdHRU = FEATUREID.y, HUC_12_bdHRU = HUC_12.y, POI_ID_bdHRU = POI_ID.y) %>%
    select(origin, FEATUREID, shared_length, perimeter, HUC_12_noPOI, FEATUREID_bdHRU, HUC_12_bdHRU, POI_ID_bdHRU, POI_ID.x)
  
  # cats without POI surrounded by POI cats
  surrounded <- filter(all.length.df, is.na(POI_ID.x)) %>%
    # origin is the CAT from which the neighboring relationships are determined
    group_by(FEATUREID) %>%
    mutate(n_cats = n()) %>%
    filter(perimeter == sum(shared_length), !is.na(POI_ID_bdHRU), n_distinct(POI_ID_bdHRU) == 1) %>%
    select(FEATUREID, POI_ID_bdHRU) %>%
    distinct()
  
  # populate cats surrounded by HRUs with HRUs POI_ID
  cats <- cats %>%
    left_join(surrounded, by = "FEATUREID") %>%
    mutate(POI_ID = ifelse(!is.na(POI_ID_bdHRU), POI_ID_bdHRU,  POI_ID))
  
  # catchments taht are partially surrounded
  part_surrounded <- filter(all.length.df, !FEATUREID %in% surrounded$FEATUREID) %>%
    # origin is the CAT from which the neighboring relationships are determined
    group_by(FEATUREID, POI_ID_bdHRU) %>%
    summarize(perim_POI = sum(shared_length)/perimeter) %>%
    distinct() %>%
    # Join the catchment to missing outlets if the same COMID for another
    #      source of information
    inner_join(missing_outlets, by = c("FEATUREID" = "COMID")) %>% 
    # if the distance to the HRU with the POI_ID populated is greater than the
    #   distance from the outlet to the edge, disregard
    filter(!dist2_HRU > dist2edge, POI_ID_bdHRU == POI_ID) %>%
    st_as_sf()
  
  # for each surrounded catchment
  for (COMID in part_surrounded$FEATUREID){
    # get the upstream flowline/cat contributiong if it exists
    up <- get_UT(nhd, COMID)
    
    nhdsub <- filter(nhd, COMID %in% up) %>%
      select(COMID) %>%
      mutate(new_POI_ID = filter(part_surrounded, FEATUREID %in% up)$POI_ID)
    
    cats <- cats %>%
      left_join(st_drop_geometry(nhdsub), by = c("FEATUREID" = "COMID")) %>%
      mutate(POI_ID = ifelse(FEATUREID %in% up, nhdsub$new_POI_ID, POI_ID)) %>%
      select(-new_POI_ID)
  }
  
  return(cats) 
}