Merge branch 'main' into 'main'

navigation edge cases

See merge request !158

workspace/02_NHD_navigate.Rmd

@@ -343,8 +343,10 @@ mapview(filter(tmp_POIs, Type_NID != ""), layer.name = "NID POIs", col.regions =
 #  Derive or load Waterbody POIs ----------------------
 if(all(is.na(tmp_POIs$Type_WBIn))) {
 
-  wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer)
-  tmp_POIs <- wb_layers$POIs
+  wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
+  
+  wb_layers <- wbin_POIcreation(nhd, WBs_VPU, data_paths, crs, split_layer, wbout_COMIDs)
+  tmp_POIs <- addType(wb_layers$POIs, tmp_POIs, "WBIn")
   
   if(!all(is.na(wb_layers$events))) {
     wb_inlet_events <- select(wb_layers$events, -c(id, offset, Hydroseq, ToMeas)) %>%

workspace/R/NHD_navigate.R

@@ -713,6 +713,10 @@ split_tt <- function(in_POI_ID, split_DF, tt_diff, max_DA){
     filter(AreaSqKM > 0) %>%
     mutate(inc_seg_area = c(TotDASqKM[1], (TotDASqKM - lag(TotDASqKM))[-1]))
   
+  if(nrow(split_tt_DF) == 0) {
+    return(NULL)
+  }
+  
   first_iter <- unique(split_tt_DF$iter)
   # Iterate through segs that need splitting
   for (i in c(first_iter:max(split_DF$iter, na.rm = T))){
@@ -1005,15 +1009,22 @@ WB_event <- function(WBs, nhd_wb, type){
       summarize(do_union = F) %>%
       st_cast("LINESTRING")
     
-    inlet_pnts <- sf::st_intersection(inlet_ls, WBs_layer) %>%
-      #group_by("LINESTRING") %>%  
-      st_cast("POINT") %>%
-      st_as_sf()
-      #mutate(row_id = row_number()) %>%
-      #group_by(LevelPathI) %>%
-      #filter(row_id == min(row_id)) %>%
-      #ungroup() 
-    
+    inlet_pnts <- sf::st_intersection(inlet_ls, WBs_layer)
+
+    # the intersection can return linestring features.    
+    if(sf::st_geometry_type(inlet_pnts, by_geometry = FALSE) == "LINESTRING") {
+      inlet_pnts <- inlet_pnts %>%
+        st_cast("POINT") %>%
+        st_as_sf() %>%
+        group_by(COMID) %>%
+        filter(row_number() == 1) %>%
+        ungroup() 
+    } else {
+      inlet_pnts <- inlet_pnts %>%
+        st_cast("POINT") %>%
+        st_as_sf()
+    }
+
     wb_events <- get_flowline_index(nhd_wb, inlet_pnts) %>%
       inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI, LevelPathI), by = "COMID") %>%
       mutate(event_type = type) %>%
@@ -1122,7 +1133,8 @@ wbout_POI_creaton <- function(nhd, WBs_VPU, data_paths, crs){
   
 }
 
-wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, split_layer){
+wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, 
+                             split_layer, wbout_COMIDs){
   # Create waterbody inlet POIs
   wbin_COMIDs <- filter(nhd, WB == 0, 
                         DnHydroseq %in% filter(nhd, WB == 1)$Hydroseq) %>%
@@ -1166,10 +1178,8 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, split_layer){
     }
   }
   
-  wbout_COMIDs <- filter(tmp_POIs, !is.na(Type_WBOut))
-  tmp_POIs <- POI_creation(filter(st_drop_geometry(wbin_COMIDs), !COMID %in% wbout_COMIDs$COMID), 
-                           nhd, "WBIn") %>%
-    addType(., tmp_POIs, "WBIn")
+  wbin_POIs <- POI_creation(filter(st_drop_geometry(wbin_COMIDs), !COMID %in% wbout_COMIDs$COMID), 
+                           nhd, "WBIn")
   
   # Get NHDArea and HR waterbodies
   WBs_VPU_areaHR <- HR_Area_coms(nhd, WBs_VPU, data_paths, crs)
@@ -1193,7 +1203,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, split_layer){
     if (nrow(nhd_inlet) > 0){
       # Get inlet events and bind with outlets
       wb_inlet_events <- WB_event(WBs_VPU_areaHR, nhd_inlet, "inlet") %>%
-        st_compatibalize(tmp_POIs) %>%
+        st_compatibalize(wbin_POIs) %>%
         mutate(nexus = TRUE) %>%
         inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq, ToMeas), by = "COMID")
       
@@ -1210,27 +1220,30 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs, split_layer){
           inner_join(select(st_drop_geometry(wbin_fl_upstream), usCOMID = COMID, toCOMID), by = c("COMID" = "toCOMID")) %>%
           mutate(dsCOMID = COMID, COMID = usCOMID)
         
-        tmp_POIs <- POI_creation(select(st_drop_geometry(wb_inlet_POIs), COMID, Type_WBIn = WBAREACOMI), 
-                                 nhd, "WBIn") %>%
-          addType(., tmp_POIs, "WBIn")
-        
-        wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
+        if(nrow(wb_inlet_POIs) > 0) {
+          wbin_POIs <- bind_rows(POI_creation(select(st_drop_geometry(wb_inlet_POIs), COMID, Type_WBIn = WBAREACOMI), 
+                                              nhd, "WBIn"), wbin_POIs)
+          
+          wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
+        }
       }
 
       if(nrow(wb_inlet_events) > 0){
         
-        tmp_POIs <- data.table::rbindlist(list(tmp_POIs,
-                                               select(wb_inlet_events, COMID, Type_WBIn = WBAREACOMI, nexus)), fill = TRUE) %>%
+        wbin_POIs <- data.table::rbindlist(list(wbin_POIs, 
+                                                select(wb_inlet_events, COMID, 
+                                                       Type_WBIn = WBAREACOMI, nexus)), 
+                                           fill = TRUE) %>%
           st_as_sf()
       }
       
-      return(list(POIs = tmp_POIs, events = wb_inlet_events))
+      return(list(POIs = wbin_POIs, events = wb_inlet_events))
     } else {
       print("no waterbody inlets events")
-      return(list(POIs = tmp_POIs, events = NA))
+      return(list(POIs = NA, events = NA))
     }
   } else {
-    return(list(POIs = tmp_POIs, events = NA))
+    return(list(POIs = NA, events = NA))
   }
   
 }