From cf39c1434c57165204735c993eecd892fc69b6dc Mon Sep 17 00:00:00 2001
From: mikejohnson51 <mikecp11@gmail.com>
Date: Fri, 15 Nov 2024 14:42:59 -0700
Subject: [PATCH] add flowline accounting

---
 NAMESPACE                      |   2 +-
 R/aggregate_along_mainstems.R  | 325 ++++++++++++++++++---------------
 R/aggregate_to_distribution.R  | 199 +++++++++++---------
 R/make_hf_gpkg_from_other.R    |   3 -
 R/modify_collapse_headwaters.R |  22 ++-
 5 files changed, 308 insertions(+), 243 deletions(-)

diff --git a/NAMESPACE b/NAMESPACE
index e118998e..7e381a90 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -32,11 +32,11 @@ export(drop_extra_features)
 export(flowpaths_to_linestrings)
 export(flush_prefix)
 export(get_minimal_network)
-export(gpkg)
 export(hl_to_outlet)
 export(hyaggregate_log)
 export(layer_exists)
 export(make_hf_gpkg_from_refactor)
+export(make_hf_gpkg_from_reference)
 export(make_hf_gpkg_from_uniform_aggregate)
 export(map_outlet_ids)
 export(middle_massage)
diff --git a/R/aggregate_along_mainstems.R b/R/aggregate_along_mainstems.R
index 871b09b1..557ce7dd 100644
--- a/R/aggregate_along_mainstems.R
+++ b/R/aggregate_along_mainstems.R
@@ -1,18 +1,18 @@
-st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))
+st_erase = function(x, y)
+  st_difference(x, st_union(st_combine(y)))
 
-add_network_type = function(network_list, verbose = TRUE){
-  
-  network_list$flowpaths = network_list$flowpaths %>% 
-    mutate(has_divide = id %in% network_list$catchments$id) %>% 
+add_network_type = function(network_list, verbose = TRUE) {
+  network_list$flowpaths = network_list$flowpaths %>%
+    mutate(has_divide = id %in% network_list$catchments$id) %>%
     filter(!duplicated(.))
   
-  network_list$catchments = network_list$catchments %>% 
-    mutate(has_flowpath = id %in% network_list$flowpaths$id) %>% 
+  network_list$catchments = network_list$catchments %>%
+    mutate(has_flowpath = id %in% network_list$flowpaths$id) %>%
     filter(!duplicated(.))
   
   
   network_list
-}  
+}
 
 #' Aggregate along network mainstems
 #' @description Given a set of ideal catchment sizes, plus the
@@ -36,30 +36,31 @@ aggregate_along_mainstems = function(network_list,
                                      min_length_km,
                                      verbose = TRUE,
                                      cache_file = NULL) {
-  
   hyaggregate_log("INFO", "\n---  Aggregate Along Mainstem ---", verbose)
-  hyaggregate_log("INFO", glue("ideal_size_sqkm --> {ideal_size_sqkm}"), verbose)
-  hyaggregate_log("INFO", glue("min_length_km --> {min_length_km}"),     verbose)
-  hyaggregate_log("INFO", glue("min_area_sqkm --> {min_area_sqkm}"),     verbose)
-  
-  tmp = network_list$flowpaths %>% 
-    st_drop_geometry() %>% 
-    as.data.frame() |> 
-    select(id = toid, hl_un = poi_id) %>% 
-    distinct() |> 
-    filter(!is.na(hl_un)) %>% 
-    group_by(id) %>% 
-    slice(1) %>% 
-    ungroup() 
-  
-  fline = network_list$flowpaths %>% 
-    st_drop_geometry() %>% 
-    mutate(hl_dn = poi_id) %>% 
-    left_join(tmp, by = "id") %>% 
-    distinct() 
+  hyaggregate_log("INFO",
+                  glue("ideal_size_sqkm --> {ideal_size_sqkm}"),
+                  verbose)
+  hyaggregate_log("INFO", glue("min_length_km --> {min_length_km}"), verbose)
+  hyaggregate_log("INFO", glue("min_area_sqkm --> {min_area_sqkm}"), verbose)
+  
+  tmp = network_list$flowpaths %>%
+    st_drop_geometry() %>%
+    as.data.frame() |>
+    select(id = toid, hl_un = poi_id) %>%
+    distinct() |>
+    filter(!is.na(hl_un)) %>%
+    group_by(id) %>%
+    slice(1) %>%
+    ungroup()
+  
+  fline = network_list$flowpaths %>%
+    st_drop_geometry() %>%
+    mutate(hl_dn = poi_id) %>%
+    left_join(tmp, by = "id") %>%
+    distinct()
   
   index_table = fline %>%
-    group_by(levelpathid) |> 
+    group_by(levelpathid) |>
     mutate(
       ind = cs_group(
         areasqkm,
@@ -68,41 +69,42 @@ aggregate_along_mainstems = function(network_list,
         hl_un,
         ideal_size_sqkm,
         min_area_sqkm,
-        min_length_km)
+        min_length_km
+      )
     ) %>%
     ungroup()   %>%
     group_by(levelpathid, ind) %>%
     mutate(set = cur_group_id(), n = n()) %>%
     ungroup() %>%
-    select(set, id, toid, levelpathid,
-           hydroseq, member_comid,
-           poi_id, n)
-  
-  index_table |>
-    group_by(id) |>
-    mutate(n = n()) |>
-    filter(n > 1)
-  
-
+    select(set, id, toid, levelpathid, hydroseq, member_comid, poi_id, n)
   
   v = aggregate_sets(network_list, index_table)
   
   v = add_network_type(v, verbose = verbose)
-
-  hyaggregate_log("SUCCESS",
-                  glue("Merged to idealized catchment size of {ideal_size_sqkm} sqkm: {nrow(network_list$flowpaths) - nrow(v$flowpaths)} features removed"),
-                  verbose)
   
-  if(!is.null(cache_file)) {
-    
+  v$lookup = 
+    index_table |> 
+    group_by(set) |> 
+    mutate(flowline_id = paste(id, collapse = ","),
+           member_comid = paste(member_comid, collapse = ",")) |> 
+    ungroup() |> 
+    select(id = set, flowline_id, member_comid ) |> 
+    distinct()
+  
+  hyaggregate_log(
+    "SUCCESS",
+    glue(
+      "Merged to idealized catchment size of {ideal_size_sqkm} sqkm: {nrow(network_list$flowpaths) - nrow(v$flowpaths)} features removed"
+    ),
+    verbose
+  )
+  
+  if (!is.null(cache_file)) {
     tmp = list()
     tmp$aggregate_along_mainstems_catchment = v$catchments
     tmp$aggregate_along_mainstems_flowpath = v$flowpaths
-
-    write_hydrofabric(tmp,
-                      cache_file,
-                      verbose, 
-                      enforce_dm = FALSE)
+    
+    write_hydrofabric(tmp, cache_file, verbose, enforce_dm = FALSE)
     
     rm(tmp)
   }
@@ -124,21 +126,28 @@ aggregate_along_mainstems = function(network_list,
 #' @param lmin a threshold, or target, cumulative size
 #' @return a vector of length(areas) containing grouping indexes
 #' @export
-#' 
+#'
 
-cs_group <- function(areas, lengths, exclude_dn, exclude_un, ideal_size_sqkm, amin, lmin) {
-  
+cs_group <- function(areas,
+                     lengths,
+                     exclude_dn,
+                     exclude_un,
+                     ideal_size_sqkm,
+                     amin,
+                     lmin) {
   areas[is.na(areas)] = 0
   lengths[is.na(lengths)] = 0
   
-  if(length(areas) == 1){ return(1) }
+  if (length(areas) == 1) {
+    return(1)
+  }
   
   break_index = which(!is.na(exclude_dn))
   break_index2 = which(!is.na(exclude_un)) - 1
   
   break_index = sort(c(break_index, break_index2))
   
-  if(length(break_index) != 0){
+  if (length(break_index) != 0) {
     sub_areas = splitAt(areas, break_index)
     sub_lengths = splitAt(lengths, break_index)
   } else {
@@ -146,21 +155,39 @@ cs_group <- function(areas, lengths, exclude_dn, exclude_un, ideal_size_sqkm, am
     sub_lengths = list(lengths)
   }
   
-  if(all(lengths(sub_areas) != lengths(sub_lengths))){
+  if (all(lengths(sub_areas) != lengths(sub_lengths))) {
     stop("Yuck~")
   }
   
   o1 = lapply(sub_areas, assign_id, athres = ideal_size_sqkm)
-
-  o2 = lapply(1:length(sub_areas),   function(i) { pinch_sides(   x = sub_areas[[i]],   ind = o1[[i]], thres = amin) })
-
-  o3 = lapply(1:length(sub_lengths), function(i) { pinch_sides(   x = sub_lengths[[i]], ind = o2[[i]], thres = lmin) })
-
-  o4 = lapply(1:length(sub_areas),   function(i) { middle_massage(x = sub_areas[[i]],   ind = o3[[i]], thres = amin) })
-
-  o5 = lapply(1:length(sub_lengths), function(i) { middle_massage(x = sub_lengths[[i]], ind = o4[[i]], thres = lmin) })
   
-  for(i in 1:length(o5)){ o5[[i]] = o5[[i]] + 1e9*i }
+  o2 = lapply(1:length(sub_areas), function(i) {
+    pinch_sides(x = sub_areas[[i]],
+                ind = o1[[i]],
+                thres = amin)
+  })
+  
+  o3 = lapply(1:length(sub_lengths), function(i) {
+    pinch_sides(x = sub_lengths[[i]],
+                ind = o2[[i]],
+                thres = lmin)
+  })
+  
+  o4 = lapply(1:length(sub_areas), function(i) {
+    middle_massage(x = sub_areas[[i]],
+                   ind = o3[[i]],
+                   thres = amin)
+  })
+  
+  o5 = lapply(1:length(sub_lengths), function(i) {
+    middle_massage(x = sub_lengths[[i]],
+                   ind = o4[[i]],
+                   thres = lmin)
+  })
+  
+  for (i in 1:length(o5)) {
+    o5[[i]] = o5[[i]] + 1e9 * i
+  }
   
   unlist(o5)
   
@@ -175,19 +202,20 @@ cs_group <- function(areas, lengths, exclude_dn, exclude_un, ideal_size_sqkm, am
 #' @importFrom sf st_set_geometry st_geometry
 #' @export
 
-pinch_sides = function(x, ind, thres){
-  
+pinch_sides = function(x, ind, thres) {
   tmp_areas = unlist(lapply(split(x, ind), sum))
   
-  if(length(tmp_areas) == 1){ return(ind) }
+  if (length(tmp_areas) == 1) {
+    return(ind)
+  }
   
   n = as.numeric(names(tmp_areas))
   
-  if(tmp_areas[1] < thres){
+  if (tmp_areas[1] < thres) {
     names(tmp_areas)[1] = names(tmp_areas[2])
   }
   
-  if(tmp_areas[length(tmp_areas)] < thres){
+  if (tmp_areas[length(tmp_areas)] < thres) {
     names(tmp_areas)[length(tmp_areas)] = names(tmp_areas[length(tmp_areas) - 1])
   }
   
@@ -206,18 +234,19 @@ pinch_sides = function(x, ind, thres){
 #' @return a vector of length(x) containing grouping indexes
 #' @export
 
-middle_massage = function(x, index_values, threshold){
-  
+middle_massage = function(x, index_values, threshold) {
   tmp_areas = unlist(lapply(split(x, index_values), sum))
   
-  if(length(tmp_areas) == 1){ return(index_values) }
+  if (length(tmp_areas) == 1) {
+    return(index_values)
+  }
   
   n = as.numeric(names(tmp_areas))
   
-  if(any(tmp_areas < threshold)){
+  if (any(tmp_areas < threshold)) {
     tmp = which(tmp_areas < threshold)
     
-    for(j in 1:length(tmp)){
+    for (j in 1:length(tmp)) {
       base = as.numeric(tmp[j])
       edges = c(base - 1, base + 1)
       becomes = names(which.min(tmp_areas[edges]))
@@ -244,37 +273,39 @@ middle_massage = function(x, index_values, threshold){
 #' @export
 
 agg_length_area   <- function(l, a, lthres, athres) {
-  
   ids = 1:length(l)
   
-  if(length(ids) != 1){
-    
-    if(!is.null(lthres)){
-      for (i in 1:(length(l)-1)) {
+  if (length(ids) != 1) {
+    if (!is.null(lthres)) {
+      for (i in 1:(length(l) - 1)) {
         if (l[i] < lthres) {
-          ids[(i+1):length(l)] = ids[(i+1):length(l)] - 1
-          l[i+1] = l[i] + l[i+1]
-          l[i]   = l[i+1]
-          a[i+1] = a[i] + a[i+1]
-          a[i] =   a[i+1]
+          ids[(i + 1):length(l)] = ids[(i + 1):length(l)] - 1
+          l[i + 1] = l[i] + l[i + 1]
+          l[i]   = l[i + 1]
+          a[i + 1] = a[i] + a[i + 1]
+          a[i] =   a[i + 1]
         }
       }
     }
     
-    if(!is.null(athres)){
-      for (i in 1:(length(a)-1)) {
+    if (!is.null(athres)) {
+      for (i in 1:(length(a) - 1)) {
         if (a[i] < athres) {
-          ids[(i+1):length(a)] = ids[(i+1):length(a)] - 1
-          a[i+1] = a[i] + a[i+1]
-          a[i] =   a[i+1]
+          ids[(i + 1):length(a)] = ids[(i + 1):length(a)] - 1
+          a[i + 1] = a[i] + a[i + 1]
+          a[i] =   a[i + 1]
         }
       }
     }
     
-    if(is.null(athres)){ athres = 0 }
-    if(is.null(lthres)){ lthres = 0 }
+    if (is.null(athres)) {
+      athres = 0
+    }
+    if (is.null(lthres)) {
+      lthres = 0
+    }
     
-    if(a[length(a)] < athres | l[length(l)] < lthres){
+    if (a[length(a)] < athres | l[length(l)] < lthres) {
       ids[length(ids)] = pmax(1, ids[length(ids)] - 1)
     }
   }
@@ -291,7 +322,7 @@ agg_length_area   <- function(l, a, lthres, athres) {
 #' @noRd
 
 splitAt <- function(x, pos) {
-  unname(split(x, cumsum(seq_along(x) %in% (pos + 1)) ))
+  unname(split(x, cumsum(seq_along(x) %in% (pos + 1))))
 }
 
 #' Index a Vector by Cumulative Sum
@@ -302,8 +333,7 @@ splitAt <- function(x, pos) {
 #' @return a vector of length(a)
 #' @export
 
-assign_id = function(x, athres){
-  
+assign_id = function(x, athres) {
   cumsum <- 0
   group  <- 1
   result <- numeric()
@@ -331,21 +361,30 @@ assign_id = function(x, athres){
 #' @importFrom nhdplusTools get_sorted
 
 aggregate_sets = function(network_list, index_table) {
-  
   set_topo = index_table %>%
     group_by(set) %>%
-    mutate(member_comid  = paste(member_comid, collapse = ","),
-           poi_id  = paste(poi_id[!is.na(poi_id)], collapse = ","),
-           poi_id  = ifelse(poi_id == "", NA, poi_id)) %>%
+    mutate(
+      member_comid  = paste(member_comid, collapse = ","),
+      poi_id  = paste(poi_id[!is.na(poi_id)], collapse = ","),
+      poi_id  = ifelse(poi_id == "", NA, poi_id)
+    ) %>%
     arrange(hydroseq) %>%
-    select(set, id, toid, levelpathid, poi_id,
-           hydroseq, member_comid) %>%
+    select(set, id, toid, levelpathid, poi_id, hydroseq, member_comid) %>%
     ungroup()
   
-  set_topo_fin = left_join(select(set_topo, set, id = toid, hydroseq,
-                                  levelpathid, poi_id, member_comid),
-                           select(set_topo, toset = set, id),
-                           by = "id") %>%
+  set_topo_fin = left_join(
+    select(
+      set_topo,
+      set,
+      id = toid,
+      hydroseq,
+      levelpathid,
+      poi_id,
+      member_comid
+    ),
+    select(set_topo, toset = set, id),
+    by = "id"
+  ) %>%
     group_by(set) %>%
     mutate(toset = ifelse(is.na(toset), 0, toset)) %>%
     filter(set != toset) %>%
@@ -361,22 +400,22 @@ aggregate_sets = function(network_list, index_table) {
       st_as_sf() %>%
       select(set) %>%
       rename_geometry("geometry")
-  }, error = function(e){
+  }, error = function(e) {
     NULL
-  }) 
-
+  })
+  
   flowpaths_out  = tryCatch({
     filter(index_table, n > 1) %>%
-    inner_join(network_list$flowpaths, by = "id") %>%
-    st_as_sf() %>%
-    select(set) %>%
-    union_linestrings('set') %>%
-    rename_geometry("geometry")
-  }, error = function(e){
+      inner_join(network_list$flowpaths, by = "id") %>%
+      st_as_sf() %>%
+      select(set) %>%
+      union_linestrings('set') %>%
+      rename_geometry("geometry")
+  }, error = function(e) {
     NULL
   })
   
-  flowpaths_out = flowpaths_out |> 
+  flowpaths_out = flowpaths_out |>
     bind_rows(single_flowpaths) %>%
     left_join(set_topo_fin, by = "set") %>%
     rename(id = set, toid = toset)
@@ -388,9 +427,9 @@ aggregate_sets = function(network_list, index_table) {
       st_as_sf() %>%
       select(set) %>%
       rename_geometry("geometry")
-  }, error = function(e){
+  }, error = function(e) {
     NULL
-  }) 
+  })
   
   catchments_out  = tryCatch({
     filter(index_table, n != 1) %>%
@@ -398,51 +437,49 @@ aggregate_sets = function(network_list, index_table) {
       st_as_sf() %>%
       select(set) %>%
       union_polygons('set') %>%
-      mutate(areasqkm = add_areasqkm(.)) 
-  }, error = function(e){
+      mutate(areasqkm = add_areasqkm(.))
+  }, error = function(e) {
     NULL
   })
   
-  catchments_out = catchments_out %>% 
+  catchments_out = catchments_out %>%
     bind_rows(single_catchments) %>%
     left_join(set_topo_fin, by = "set") %>%
-    select(id = set, toid = toset) 
+    select(id = set, toid = toset)
   
   ####
-
-  if(!is.null(catchments_out)){
+  
+  if (!is.null(catchments_out)) {
     mps = suppressWarnings({
-      catchments_out %>% 
-        st_cast("MULTIPOLYGON") %>% 
-        st_cast("POLYGON") %>% 
-        fast_validity_check() %>% 
-        add_count(id) 
+      catchments_out %>%
+        st_cast("MULTIPOLYGON") %>%
+        st_cast("POLYGON") %>%
+        fast_validity_check() %>%
+        add_count(id)
     })
     
-    if(nrow(mps) > nrow(catchments_out)){
-      fixers = filter(mps, n > 1) %>% 
-        mutate(areasqkm = add_areasqkm(.),
-               tmpID = 1:n()) %>% 
-        group_by(id) 
+    if (nrow(mps) > nrow(catchments_out)) {
+      fixers = filter(mps, n > 1) %>%
+        mutate(areasqkm = add_areasqkm(.), tmpID = 1:n()) %>%
+        group_by(id)
       
-      keep = slice_max(fixers, areasqkm) %>% 
-        ungroup() 
+      keep = slice_max(fixers, areasqkm) %>%
+        ungroup()
       
-      to_merge = filter(fixers, !tmpID %in% keep$tmpID) %>% 
+      to_merge = filter(fixers, !tmpID %in% keep$tmpID) %>%
         ungroup()
       
-      good_to_go = filter(mps, n == 1) %>% 
-        bind_rows(select(keep, id, toid)) %>% 
+      good_to_go = filter(mps, n == 1) %>%
+        bind_rows(select(keep, id, toid)) %>%
         fast_validity_check()
       
       catchments_out = suppressWarnings({
-        terra::combineGeoms(vect(good_to_go), vect(to_merge)) %>% 
-          st_as_sf() %>% 
+        terra::combineGeoms(vect(good_to_go), vect(to_merge)) %>%
+          st_as_sf() %>%
           st_cast("POLYGON")
       })
-    } 
+    }
   }
   
   prepare_network(network_list = list(flowpaths = flowpaths_out, catchments = catchments_out))
 }
-
diff --git a/R/aggregate_to_distribution.R b/R/aggregate_to_distribution.R
index 9278e362..c0d100da 100644
--- a/R/aggregate_to_distribution.R
+++ b/R/aggregate_to_distribution.R
@@ -41,7 +41,6 @@ aggregate_to_distribution = function(gpkg = NULL,
                                      overwrite = FALSE,
                                      cache = FALSE,
                                      verbose = TRUE) {
-  
   if (cache &  is.null(outfile)) {
     stop("cache cannot be written if outfile is NULL")
   }
@@ -74,36 +73,35 @@ aggregate_to_distribution = function(gpkg = NULL,
   network_list <- read_hydrofabric(gpkg,
                                    catchments = divide,
                                    flowpaths = flowpath,
-                                   crs = crs) |> 
-    prepare_network() |> 
+                                   crs = crs) |>
+    prepare_network() |>
     add_network_type(verbose = FALSE)
   
-  if(!"member_comid" %in% names(network_list$flowpaths)){
+  if (!"member_comid" %in% names(network_list$flowpaths)) {
     network_list$flowpaths$member_comid = network_list$flowpaths$id
   }
   
   # Add outlets
   if (!is.null(pois)) {
-    
     names(pois) = tolower(names(pois))
     
     bad = filter(pois, !id %in% network_list$flowpaths$id)
     
     outflows = pois %>%
-      st_drop_geometry() %>% 
+      st_drop_geometry() %>%
       select(poi_id, id)
-
-    network_list$flowpaths  = left_join(mutate(network_list$flowpaths, poi_id = NULL), 
-                                        st_drop_geometry(outflows), 
+    
+    network_list$flowpaths  = left_join(mutate(network_list$flowpaths, poi_id = NULL),
+                                        st_drop_geometry(outflows),
                                         by = 'id')
     
     unique(pois$poi_id) |> length()
     unique(network_list$flowpaths$poi_id) |> sort()
     
-    if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))){
-      warning("All POIs are not indexable: reference", call.  = FALSE )
+    if (length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))) {
+      warning("All POIs are not indexable: reference", call.  = FALSE)
     }
-   
+    
     
   } else {
     network_list$flowpaths$poi_id    = NA
@@ -116,15 +114,12 @@ aggregate_to_distribution = function(gpkg = NULL,
     tmp$base_catchments = network_list$catchments
     tmp$base_flowpaths = network_list$flowpaths
     
-    write_hydrofabric(tmp,
-                      cache_file,
-                      verbose, 
-                      enforce_dm = FALSE)
+    write_hydrofabric(tmp, cache_file, verbose, enforce_dm = FALSE)
     
     rm(tmp)
   }
   
-  main_agg =  
+  main_agg =
     aggregate_along_mainstems(
       network_list,
       ideal_size_sqkm,
@@ -132,10 +127,10 @@ aggregate_to_distribution = function(gpkg = NULL,
       min_length_km,
       verbose = verbose,
       cache_file = cache_file
-  ) 
+    )
   
-  if(length(unique(pois$poi_id)) != length(unique(na.omit(main_agg$flowpaths$poi_id)))){
-    warning("All POIs are not indexable: mainstem agg",call.  = FALSE )
+  if (length(unique(pois$poi_id)) != length(unique(na.omit(main_agg$flowpaths$poi_id)))) {
+    warning("All POIs are not indexable: mainstem agg", call.  = FALSE)
   }
   
   collapse_agg = collapse_headwaters2(
@@ -143,104 +138,128 @@ aggregate_to_distribution = function(gpkg = NULL,
     min_area_sqkm,
     min_length_km,
     verbose = verbose,
-    cache_file = cache_file)
+    cache_file = cache_file
+  )
   
-  collapse_agg$catchments = clean_geometry(collapse_agg$catchments, 
-                                           ID = "id", 
-                                           keep = NULL)
+  collapse_agg$catchments = clean_geometry(collapse_agg$catchments, ID = "id", keep = NULL)
   
-  if(length(unique(pois$poi_id)) != length(unique(na.omit(collapse_agg$flowpaths$poi_id)))){
-    warning("All POIs are not indexable: collapse", call.  = FALSE )
+  if (length(unique(pois$poi_id)) != length(unique(na.omit(collapse_agg$flowpaths$poi_id)))) {
+    warning("All POIs are not indexable: collapse", call.  = FALSE)
   }
   
   network_list = collapse_agg
   rm(collapse_agg)
   rm(main_agg)
-
-if(!is.null(pois)){
-  network_list$pois =  network_list$flowpaths %>% 
-    st_drop_geometry() %>%
-    select(id, poi_id) %>% 
-    mutate(poi_id = as.integer(poi_id)) |> 
-    filter(!is.na(poi_id)) %>%
-    left_join(select(pois, poi_id), by = "poi_id")
-} 
+  
+  if (!is.null(pois)) {
+    network_list$pois =  network_list$flowpaths %>%
+      st_drop_geometry() %>%
+      select(id, poi_id) %>%
+      mutate(poi_id = as.integer(poi_id)) |>
+      filter(!is.na(poi_id)) %>%
+      left_join(select(pois, poi_id), by = "poi_id")
+  }
   
   network_list$flowpaths = suppressWarnings({
-    hydroloom::add_streamorder(network_list$flowpaths) 
-  }) 
-  
-  network_list$flowpaths = 
-    select(network_list$flowpaths, 
-           id, toid, 
-           mainstem = levelpathid, 
-           order = stream_order, 
-           member_comid, 
-           poi_id, 
-           hydroseq, 
-           lengthkm, 
-           areasqkm, 
-           tot_drainage_areasqkm = tot_drainage_area, 
-           has_divide) %>% 
+    hydroloom::add_streamorder(network_list$flowpaths)
+  })
+  
+  network_list$flowpaths =
+    select(
+      network_list$flowpaths,
+      id,
+      toid,
+      flowline_id,
+      mainstem = levelpathid,
+      order = stream_order,
+      member_comid,
+      poi_id,
+      hydroseq,
+      lengthkm,
+      areasqkm,
+      tot_drainage_areasqkm = tot_drainage_area,
+      has_divide
+    ) %>%
     mutate(divide_id = ifelse(id %in% network_list$catchments$id, id, NA))
   
-  if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))){
-    warning("All POIs are not indexable: final",call.  = FALSE )
+  if (length(unique(pois$poi_id)) != length(unique(na.omit(network_list$flowpaths$poi_id)))) {
+    warning("All POIs are not indexable: final", call.  = FALSE)
   }
   
-  topo = st_drop_geometry(network_list$flowpaths) %>% 
+  topo = st_drop_geometry(network_list$flowpaths) %>%
     select(divide_id, toid)
   
-  network_list$divides = select(network_list$catchments, id,  areasqkm) %>% 
-    mutate(divide_id = id, has_flowline = TRUE, ds_id = NA, type = "network") %>% 
+  network_list$divides = select(network_list$catchments, id, areasqkm) %>%
+    mutate(
+      divide_id = id,
+      has_flowline = TRUE,
+      ds_id = NA,
+      type = "network"
+    ) %>%
     left_join(topo, by = "divide_id")
   
   network_list$catchments = NULL
-
+  
   network_list$network  = st_drop_geometry(network_list$flowpaths) %>%
-   select(
-     id,
-     toid,
-     member  = member_comid,
-     divide_id,
-     any_of('poi_id'),
-     mainstem,
-     hydroseq,
-     order,
-     lengthkm, areasqkm,
-     tot_drainage_areasqkm) %>%
-   separate_longer_delim(col = 'member', delim = ",") %>%
-   mutate(hf_id_part = sapply( strsplit(member, "[.]"), FUN = function(x){ x[2] }),
-          hf_id_part = ifelse(is.na(hf_id_part), 1L, floor(as.numeric((hf_id_part)))),
-          hf_id = sapply( strsplit(member, "[.]"), FUN = function(x){ as.numeric(x[1]) }),
-          member = NULL,
-          hf_source = "NHDPlusV2"
-   ) %>%
-   left_join(st_drop_geometry(select(network_list$divides, divide_id, type, ds_id)), by = "divide_id") |> 
+    select(
+      id,
+      toid,
+      member  = member_comid,
+      flowline_id,
+      divide_id,
+      any_of('poi_id'),
+      mainstem,
+      hydroseq,
+      order,
+      lengthkm,
+      areasqkm,
+      tot_drainage_areasqkm
+    ) %>%
+    separate_longer_delim(col = 'member', delim = ",") %>%
+    mutate(
+      hf_id_part = sapply(
+        strsplit(member, "[.]"),
+        FUN = function(x) {
+          x[2]
+        }
+      ),
+      hf_id_part = ifelse(is.na(hf_id_part), 1L, floor(as.numeric((
+        hf_id_part
+      )))),
+      hf_id = sapply(
+        strsplit(member, "[.]"),
+        FUN = function(x) {
+          as.numeric(x[1])
+        }
+      ),
+      member = NULL,
+      hf_source = "NHDPlusV2"
+    ) %>%
+    left_join(st_drop_geometry(select(
+      network_list$divides, divide_id, type, ds_id
+    )), by = "divide_id") |>
     distinct()
   
-  if(length(unique(pois$poi_id)) != length(unique(na.omit(network_list$network$poi_id)))){
-    warning("All POIs are not indexable: network", call.  = FALSE )
+  if (length(unique(pois$poi_id)) != length(unique(na.omit(network_list$network$poi_id)))) {
+    warning("All POIs are not indexable: network", call.  = FALSE)
   }
   
-  if(!is.null(vpu)){
+  if (!is.null(vpu)) {
     network_list$network$vpu = vpu
   } else {
     network_list$network$vpu = NA
   }
-
- 
- if(!all(st_geometry_type(network_list$divides) == "POLYGON")){
-   warning("MULTIPOLYGONS FOUND VPU: ", vpu)
- }
-
-  if (!is.null(outfile)) {
   
-    outfile = write_hydrofabric(
-      network_list,
-      outfile,
-      verbose = verbose, 
-      enforce_dm = FALSE)
+  
+  if (!all(st_geometry_type(network_list$divides) == "POLYGON")) {
+    warning("MULTIPOLYGONS FOUND VPU: ", vpu)
+  }
+  
+  if (!is.null(outfile)) {
+    outfile = write_hydrofabric(network_list,
+                                outfile,
+                                verbose = verbose,
+                                enforce_dm = FALSE)
     
     return(outfile)
     
diff --git a/R/make_hf_gpkg_from_other.R b/R/make_hf_gpkg_from_other.R
index 54b79de7..ef87b63a 100644
--- a/R/make_hf_gpkg_from_other.R
+++ b/R/make_hf_gpkg_from_other.R
@@ -93,9 +93,6 @@ make_hf_gpkg_from_uniform_aggregate = function(gpkg){
 #' @return file.path
 #' @export
 
-
-gpkg = '/Volumes/Transcend/ngen/CONUS-hydrofabric/01_reference/reference_01.gpkg'
-
 make_hf_gpkg_from_reference = function(gpkg){
   
   if (is.null(gpkg)) {
diff --git a/R/modify_collapse_headwaters.R b/R/modify_collapse_headwaters.R
index c98aae0f..2ee014da 100644
--- a/R/modify_collapse_headwaters.R
+++ b/R/modify_collapse_headwaters.R
@@ -36,12 +36,16 @@ build_headwater_collapse = function(network_list,
   df$hl1 = network_list$flowpaths$poi_id[match(df$id, network_list$flowpaths$id)]
   df$hl2 = network_list$flowpaths$poi_id[match(df$becomes, network_list$flowpaths$id)]
   
+  df$fl1 = network_list$flowpaths$flowline_id[match(df$id, network_list$flowpaths$id)]
+  df$fl2 = network_list$flowpaths$flowline_id[match(df$becomes, network_list$flowpaths$id)]
+  
   tmp = suppressWarnings({
     group_by(df, becomes) %>%
       mutate(member_comid = paste0(mC2[1], "," ,paste(mC1, collapse = ",")),
-             poi_id = as.numeric(paste(unique(na.omit(c(hl1, hl2))), collapse = ","))) %>%
+             poi_id = as.numeric(paste(unique(na.omit(c(hl1, hl2))), collapse = ",")),
+             fl_id = paste0(fl2[1], "," ,paste(fl1, collapse = ","))) %>%
       ungroup() %>%
-      select(-mC1, -mC2, -hl1, -hl2)
+      select(-mC1, -mC2, -hl1, -hl2, -fl1, -fl2)
   })
   
   
@@ -59,9 +63,14 @@ collapse_headwaters2 = function(network_list,
   
   hyaggregate_log("INFO", "\n--- Collapse Network Inward ---\n", verbose)
   
+  network_list$flowpaths =  network_list$flowpaths |> 
+    left_join(select(network_list$lookup, id, flowline_id), by = 'id')
+  
   start <- nrow(network_list$flowpaths)
   
-  mapping_table <- build_headwater_collapse(network_list, min_area_sqkm, min_length_km)
+  mapping_table <- build_headwater_collapse(network_list, 
+                                            min_area_sqkm, 
+                                            min_length_km)
   
   count = 0
   
@@ -74,6 +83,7 @@ collapse_headwaters2 = function(network_list,
     ind = match(mapping_table$becomes, network_list$flowpaths$id)
 
     network_list$flowpaths$member_comid[ind] = mapping_table$member_comid
+    network_list$flowpaths$flowline_id[ind] = mapping_table$fl_id
     
     fl = filter(network_list$flowpaths, !id %in% mapping_table$id)
     
@@ -90,7 +100,9 @@ collapse_headwaters2 = function(network_list,
     
     network_list  = prepare_network(network_list = list(flowpaths = fl, catchments = cat))
     
-    mapping_table = build_headwater_collapse(network_list, min_area_sqkm, min_length_km)
+    mapping_table = build_headwater_collapse(network_list, 
+                                             min_area_sqkm, 
+                                             min_length_km)
   }
   
   network_list = add_network_type(network_list, verbose)
@@ -100,7 +112,7 @@ collapse_headwaters2 = function(network_list,
   if (!is.null(cache_file)) {
     tmp = list()
     tmp$collapse_headwaters_catchments = network_list$catchments
-    tmp$collapse_headwaters_flowpaths = network_list$flowpaths
+    tmp$collapse_headwaters_flowpaths  = network_list$flowpaths
     
     write_hydrofabric(tmp,
                       cache_file,