Some minor mods, and better POI collapse

workspace/R/NHD_navigate.R

@@ -98,7 +98,7 @@ segment_creation <- function(nhdDF, routing_fix){
   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 (missing(routing_fix) == FALSE){
+  if (!is.na(routing_fix)){
     routing_fix <- routing_fix %>%
       rename(COMID = oldPOI, new_COMID = COMID)
     
@@ -207,19 +207,22 @@ DS_poiFix <- function(POIs_wgeom, nhd){
 #' 
 #' @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 <- nhdplusTools::get_UM(nhdDF, poi_fix, sort=T)
-  seg2fix <- filter(nhdDF, COMID == poi_fix)
+  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(factor(COMID, levels = up_segs)) %>%
+    arrange(Hydroseq) %>%
     filter(AreaSqKM > 0)
   
-  down_segs <- nhdplusTools::get_DM(nhdDF, poi_fix, sort=T)
+  down_segs <- unique(nhdplusTools::get_DM(nhdDF, poi_fix, sort=T))
   downstuff <- filter(nhdDF, COMID %in% unlist(down_segs)) %>% 
-    arrange(factor(COMID, levels = down_segs)) %>%
+    arrange(Hydroseq)%>%
     filter(AreaSqKM > 0)
   
   # combine into one dataframe, select up/downstream seg with least change in total drainage area
@@ -966,11 +969,16 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
   gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
     addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE) 
   
-  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()
-  
+  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))
 }
@@ -985,7 +993,7 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
 #'  
 wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
   # Create waterbody outlet POIs for waterbodies that are in NHDv2 waterbody set
-  wbout_COMIDs <- filter(nhd, dend == 1 & WB == 1) %>%
+  wbout_COMIDs <- nhd %>% #filter(nhd, dend == 1 & WB == 1)
     group_by(WBAREACOMI) %>%
     slice(which.min(Hydroseq)) %>%
     switchDiv(., nhd) %>%
@@ -1074,7 +1082,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
     filter(minNet == 1) 
   
   # Headwater Waterbodies that may need the network extended through the inlet
-  need_wbin <- st_drop_geometry(filter(WBs_VPU, source == "NHDv2WB")) %>%
+  need_wbin <- st_drop_geometry(filter(WBs_VPU, source %in% c("ref_WB", "NHDv2WB"))) %>%
     dplyr::select(COMID)%>%
     dplyr::filter(!COMID %in% wbin_COMIDs$WBAREACOMI)
   
@@ -1151,7 +1159,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
           mutate(dsCOMID = COMID, COMID = usCOMID)
         
         if(nrow(wb_inlet_POIs) > 0) {
-          wbin_POIs <- bind_rows(wbin_POIs, POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI), 
+          wbin_POIs <- bind_rows(wbin_POIs, POI_creation2(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI), 
                                               nhd, "WBIn"))
           
           wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
@@ -1208,7 +1216,7 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
       unique()
     
     if(nrow(gage_reach) == 0){
-      print("Wilton")
+      print("no gage reaches")
     }
     
     if(nrow(gage_event) == 0){
@@ -1333,11 +1341,9 @@ wb_inlet_collapse <- function(tmp_POIs, nhd, events){
 #'  @param events (sf data.frame) waterbody inlet events
 #' 
 #'  @return (sf data.frame) dataframe of wb inlet POIs collapse
-#' 
-#' wb_poi_collapse <- function(tmp_POIs, nhd, events){
 wb_poi_collapse <- function(tmp_POIs, nhd, events){
   gage_dist_node <- function(x, wb_ds_ds, gage_add, events){
-    print (x)
+    print (x) #6116850
     wb_out_fl <- filter(wb_ds_ds, COMID == x)
     gage_ds <- filter(wb_ds_ds, Hydroseq %in% wb_out_fl$Hydroseq |
                         Hydroseq %in% wb_out_fl$DnHydroseq) 
@@ -1353,15 +1359,14 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
       unique()
     
     if(nrow(gage_reach) == 0){
-      print("Wilton")
+      print("no gages")
     }
     
     if(nrow(gage_event) == 0){
-      #print("akermayun")
-      return("Akermayun")
+      return("no events")
     } else if(gage_event$COMID != wb_out_fl$COMID) {
       gage_reach <- gage_reach %>%
-        filter(REACHCODE == gage_event$reachcode) %>%
+        filter(REACHCODE == unique(gage_event$reachcode)) %>%
         mutate(gage_dist = ifelse(gage_event$nexus == TRUE,
                                   total_length * (1 - (gage_event$reach_meas/100)),
                                   total_length)) %>%
@@ -1369,10 +1374,10 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
                wbout_comid = x)
     } else if(gage_event$COMID == wb_out_fl$COMID){
       if(nrow(wb_event) >0){
-        wb_out_meas <- wb_event$REACH_meas
+        wb_out_meas <- min(wb_event$REACH_meas)
         wb_RC <- wb_event$REACHCODE
       } else {
-        wb_out_meas <- wb_out_fl$FromMeas
+        wb_out_meas <- min(wb_out_fl$FromMeas)
         wb_RC <- wb_out_fl$REACHCODE
       }
       
@@ -1383,7 +1388,7 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
       
       # wb info
       wb_reach <- gage_reach %>%
-        filter(REACHCODE == wb_RC) %>%
+        filter(REACHCODE == unique(wb_RC)) %>%
         mutate(wb_dist = total_length * (1 - (wb_out_meas /100)))
       
       gage_reach <- gage_reach %>%
@@ -1393,14 +1398,11 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
     }
   }
   
-  #events <- read_sf(temp_gpkg, split_layer) %>%
-  #  rbind(st_compatibalize(wb_,.))
   
   # Previously identified streamgages within Gage_Selection.Rmd
   streamgages_VPU <- gages %>%
     rename(COMID = comid) %>%
     filter(COMID %in% nhd$COMID) %>%
-    #st_drop_geometry() %>%
     switchDiv(., nhd) 
   
   # get waterbody outlets
@@ -1421,10 +1423,16 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
   gage_add <- filter(streamgages_VPU, site_no %in% gage_wb$Type_Gages) %>%
     select(COMID, reachcode, reach_meas, site_no) %>%
     inner_join(select(st_drop_geometry(gage_wb), site_no = Type_Gages, nexus), 
-               by = "site_no")
+               by = "site_no") %>%
+    distinct()
   
   output <- lapply(wb_out$COMID, 
                    gage_dist_node, wb_ds_ds, gage_add, events)
+  output_length <- output[lengths(output) > 1]
+  
+  if(length(output_length) == 0){
+    return(list(POIs = tmp_POIs, events_ret = NA))
+  }
   
   output_full <- do.call("rbind", output[lengths(output) > 1]) %>%
     filter(gage_dist < 1) %>%
@@ -1438,12 +1446,11 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
   
   gage_POI <- filter(tmp_POIs, COMID %in% output_full$gage_comid) %>%
     select(COMID, Type_HUC12_ds = Type_HUC12, Type_Gages_ds = Type_Gages, 
-           Type_TE_ds = Type_TE, Type_Term_ds = Type_Term, nexus) %>%
+           Type_TE_ds = Type_TE, nexus) %>%
     st_drop_geometry() %>%
     group_by(COMID) %>%
     summarise(Type_HUC12_ds = last(na.omit(Type_HUC12_ds)), 
               Type_Gages_ds = last(na.omit(Type_Gages_ds)),
-              Type_Term_ds = last(na.omit(Type_Term_ds)),
               Type_TE_ds = last(na.omit(Type_TE_ds)),
               nexus = last(na.omit(nexus)))
   
@@ -1452,9 +1459,8 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
     inner_join(select(gage_POI, -nexus), by = c("gage_comid" = "COMID")) %>%
     mutate(Type_HUC12 = ifelse(!is.na(Type_HUC12_ds), Type_HUC12_ds, Type_HUC12),
            Type_Gages = ifelse(!is.na(Type_Gages_ds), Type_Gages_ds, Type_Gages),
-           Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE),
-           Type_Term = ifelse(!is.na(Type_Term_ds), Type_Term_ds, Type_Term)) %>%
-    select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds, Type_Term_ds))
+           Type_TE = ifelse(!is.na(Type_TE_ds), Type_TE_ds, Type_TE)) %>%
+    select(-c(Type_HUC12_ds, Type_Gages_ds, Type_TE_ds))
   
   tmp_POIs_fin <- filter(tmp_POIs, !COMID %in% c(WB_POI$COMID, WB_POI$gage_comid)) %>%
     rbind(select(WB_POI, -gage_comid))
@@ -1466,4 +1472,135 @@ wb_poi_collapse <- function(tmp_POIs, nhd, events){
   
   
   return(list(POIs = tmp_POIs_fin, events_ret = events))
+}
+
+#'  Creates Waterbody POIs, calls a few other functions
+#'  @param (pois) sf data frame of POIs
+#'  @param move_category (character) POI data theme to move
+#'  @param DAR (numeric) drainage area threshold to move within
+#'  @param dist (numeric) maximum river distance between two points to move within
+#'  @param keep_category (character) POI data themes to keep static
+#' 
+#'  @return (sf data.frame, table) dataframe of pois, table of points that have moved
+poi_move <- function(pois, move_category, DAR, dist, keep_category) {
+  # filter out features with identical geometry
+  
+  # Add row_number
+  pois_edit <- pois %>%
+    mutate(nexus = ifelse(is.na(nexus), 0, nexus))
+  
+  # Don't consider points already moved
+  if("moved" %in% colnames(pois_edit)){
+    pois_tomove <- filter(pois_edit, is.na(moved)) # change from poi_edit
+    pois_moved_pre <- filter(pois_edit, !is.na(moved))}
+  
+  # If 'keep' category included
+  if(!missing(keep_category)){
+    poi2move <- filter(pois_tomove, !is.na(.data[[move_category]]) & nexus == FALSE) %>%
+      filter(if_all(!!as.symbol(keep_category), function(x) is.na(x))) %>%
+      # Never move these
+      filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.)))
+    
+    pois2keep <- filter(pois_tomove, !id %in% poi2move$id) 
+    #is.na(.data[[move_category]]) & nexus == FALSE) #%>%
+    #filter(if_all(!!as.symbol(keep_category), function(x) is.na(x)))
+  } else {
+    # POIs to move
+    poi2move <- pois_tomove %>%
+      filter_at(vars(Type_WBOut, Type_WBIn, Type_Conf, Type_Term), all_vars(is.na(.))) %>%
+      filter(nexus == 0) %>%
+      filter(!is.na(.data[[move_category]]))
+    
+    pois2keep <- filter(pois_tomove, !id %in% poi2move$id)
+  }
+  
+  # Get relevant NHD data
+  nhd_poi1 <- filter(nhd, COMID %in% pois2keep$COMID)
+  nhd_poi2 <- filter(nhd, COMID %in% poi2move$COMID)
+  # Ensure they are on same level path
+  nhd_poi2 <- filter(nhd_poi2, LevelPathI %in% nhd_poi1$LevelPathI)
+  
+  # Join NHD data
+  pois2keep_nhd <- pois2keep %>% 
+    inner_join(select(st_drop_geometry(nhd_poi1), COMID, LevelPathI, Hydroseq,
+                      DA_keep = TotDASqKM, Pathlength_keep = Pathlength), by = "COMID") %>%
+    rename(COMID_keep = COMID)
+  
+  # Join NHD data
+  pois2move_nhd <- select(poi2move, COMID, !!as.symbol(move_category), id_move = id) %>% 
+    inner_join(select(st_drop_geometry(nhd_poi2), COMID, LevelPathI, Hydroseq, TotDASqKM, Pathlength), 
+               by = "COMID")
+  
+  # Candidates to move
+  pois2move_cand <-pois2move_nhd %>%
+    inner_join(select(st_drop_geometry(pois2keep_nhd), COMID_keep, DA_keep, LevelPathI,
+                      Pathlength_keep, id_keep = id, nexus), 
+               by = "LevelPathI") %>%
+    mutate(river_dist = abs(Pathlength - Pathlength_keep), DAR_poi = abs(DA_keep/TotDASqKM),
+           move_dir = ifelse(Pathlength < Pathlength_keep, "Up", "Down")) %>%
+    group_by(id_move, move_dir) %>%
+    ungroup() %>%
+    filter((river_dist < dist) & (DAR_poi > (1 - DAR)) & (DAR_poi < (1 + DAR))) %>%
+    select(!!as.symbol(move_category), id_move, COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
+    st_drop_geometry()
+  
+  move_distinct <- pois2move_cand %>%
+    group_by(id_keep) %>%
+    filter(row_number() == 1) %>%
+    ungroup() %>%
+    distinct(id_move, COMID_move = COMID, id_keep, COMID_keep, river_dist, DAR_poi, move_dir, nexus) %>%
+    group_by(id_move) %>%
+    slice(which.min(abs(1 - DAR_poi))) 
+  
+  if(nrow(move_distinct) == 0){
+    print("no POIs to move")
+    return(pois)
+  }
+  
+  pois2_move <- filter(st_drop_geometry(pois_tomove), id %in% move_distinct$id_move) %>%
+    select_if(~sum(!is.na(.)) > 0) %>%
+    select(-c(COMID, nexus)) %>%
+    inner_join(select(move_distinct, id_move, id_keep), by = c("id" = "id_move"))
+  
+  move_fields <- colnames(select(pois2_move, -c(id, id_keep)))
+  
+  if(length(move_fields) == 1){
+    pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
+      inner_join(select(pois2_move, id_move = id, id_keep, 
+                        new_val = !!as.symbol(move_category)), by = c("id" = "id_keep")) %>%
+      mutate(moved := ifelse(is.na(!!as.symbol(move_category)),
+                             id_move, moved),
+             !!as.symbol(move_category) := ifelse(is.na(!!as.symbol(move_category)),
+                                                  new_val, !!as.symbol(move_category)))
+    
+    moved_points <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
+      mutate(moved_value = move_category)
+  } else {
+    for (field in move_fields){
+      pois2_keep <- filter(pois_tomove, id %in% pois2_move$id_keep, !id %in% pois2_move$id) %>%
+        inner_join(select(pois2_move, id_move = id, id_keep, new_val = !!as.symbol(field)), 
+                   by = c("id" = "id_keep")) %>%
+        mutate(moved := ifelse(is.na(!!as.symbol(field)),
+                               id_move, moved),
+               !!as.symbol(field) := ifelse(is.na(!!as.symbol(field)),
+                                            new_val, !!as.symbol(field)))
+      
+      pois_moved <- filter(pois2_keep, !is.na(new_val), !is.na(moved)) %>%
+        mutate(moved_value = field)
+      
+      if(!exists("moved_points")){
+        moved_points <- pois_moved
+      } else {
+        moved_points <- rbind(moved_points, pois_moved)
+      }
+    }
+  }
+  
+  
+  pois_final <- data.table::rbindlist(list(filter(pois_edit, !id %in% c(moved_points$id_move, pois2_keep$id)),
+                                           select(pois2_keep, -c(new_val, id_move, new_val))), fill = TRUE) %>%
+    st_as_sf()
+  
+  return(list(final_pois = pois_final, moved_points = moved_points))
+  
 }
\ No newline at end of file