clean up for #78

workspace/06-2_aggregate_cats.Rmd

@@ -51,54 +51,50 @@ reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs)
 # refactored has the original COMIDs and other attributes
 refactored <- st_transform(read_sf(out_refac_gpkg, refactored_layer), crs) 
 
-# Summarize total drainage area for reconciled flowlines
-reconciled_DA <- select(reconciled, ID, toID) %>%
-  left_join(select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
-  mutate(ID = as.integer(ID), toID = as.integer(toID),
-         area = ifelse(is.na(area), 0, area)) %>%
-  mutate(rec_DA = calculate_total_drainage_area(.)) %>%
-  left_join(select(lookup, reconciled_ID, NHDPlusV2_COMID), by = c("ID" = "reconciled_ID")) %>%
-  mutate(NHDPlusV2_COMID = as.integer(NHDPlusV2_COMID))
-  
 # Create unified POI/outlet layer
 if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
   
-  # subset reconciled flowlines by minimum drainage area for aggcats
-  reconciled_sub <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) 
+  # paths that will for sure be "in network" after aggregate
+  dominant_paths <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)
   
   # Get the end of the reconciled flowlines and bind with data frame
-  reconciled_end <- get_node(filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) , "end") %>%
-    cbind(st_drop_geometry(filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm) ))
+  reconciled_end <- get_node(dominant_paths, "end") %>%
+    cbind(st_drop_geometry(dominant_paths))
   
   # If there are events on split flowlines
   if(nrow(read_sf(out_refac_gpkg, split_layer)) > 0){
     events <- read_sf(out_refac_gpkg, split_layer) 
-    
-    # Get rows of lookup for split events
-    lookup_events <- filter(lookup, NHDPlusV2_COMID %in% events$COMID)
-    reconciled_events <- filter(reconciled_end, ID %in% lookup_events$reconciled_ID)
 
     # Outlets that are not associated with flowlines with events
-    poi_outlets <- filter(outlets_all, !Type_Gages %in% events$site_no) %>%
+    poi_outlets_mapped <- filter(outlets_all, !Type_Gages %in% events$site_no) %>%
       # Reset Type_Gages attrbute to NA if event exists
       distinct(.keep_all = TRUE) %>%
-      filter(TotDASqKM > aggregate_da_thresh_sqkm)
-    
-    # Map non-event outlets ID to reconciled ID 
-    poi_outlets_mapped <- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)) %>%
+      filter(TotDASqKM > aggregate_da_thresh_sqkm) %>%
+      st_drop_geometry() %>%
+      select(COMID, type) %>%
+      map_outlet_ids(dominant_paths) %>%
       mutate(COMID = as.integer(COMID)) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
       st_as_sf() 
-
+    
+    # Summarize total drainage area for reconciled flowlines
+    reconciled_DA <- left_join(select(st_drop_geometry(reconciled), ID, toID),
+                               select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
+      mutate(ID = as.integer(ID), toID = as.integer(toID), area = ifelse(is.na(area), 0, area)) %>%
+      mutate(rec_DA = calculate_total_drainage_area(.)) %>%
+      left_join(select(lookup, reconciled_ID, NHDPlusV2_COMID), by = c("ID" = "reconciled_ID")) %>%
+      mutate(NHDPlusV2_COMID = as.integer(NHDPlusV2_COMID)) %>%
+      select(ID, rec_DA, NHDPlusV2_COMID)
+    
     # Pois which share same COMID as events, reset Type_Gages to NA 
     poi_events_mapped <- filter(outlets_all, COMID %in% events$COMID) %>%
       mutate(Type_Gages = NA) %>%
       distinct() %>%
       # Subset to POIs that have other POI flags than Type_Gage
-      filter_at(vars(-c(COMID, TotDASqKM, type, geom)), any_vars(!is.na(.)), TotDASqKM > 1) %>%
-      inner_join(select(st_drop_geometry(reconciled_DA), ID, rec_DA, NHDPlusV2_COMID), 
+      filter_at(vars(-c(COMID, TotDASqKM, type, geom)), any_vars(!is.na(.)), 
+                TotDASqKM > aggregate_da_thresh_sqkm) %>%
+      inner_join(reconciled_DA, 
                  by = c("COMID" = "NHDPlusV2_COMID")) %>% 
       group_by(COMID) %>% 
       # Max DA associated with each orig_COMID will be the proper reconciled ID value
@@ -110,7 +106,7 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
     events_mapped <- select(st_drop_geometry(events), COMID, Type_Gages = site_no) %>%
       #IDs that events are associated with
       mutate(ID = reconciled_end[st_nearest_feature(st_transform(events, st_crs(reconciled_end)), 
-                                                 reconciled_end),]$ID) %>%
+                                                    reconciled_end),]$ID) %>%
       mutate(type = "outlet", COMID = as.character(COMID)) %>%
       # Replace event geometry (off-flowline) with POI geometry (on-flowline)
       inner_join(reconciled_end, by = "ID") %>%
@@ -145,7 +141,7 @@ if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
     
     # Map non-event outlets ID to reconciled ID 
     outlets_POI<- map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), 
-                                     filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm)) %>%
+                                 dominant_paths) %>%
       inner_join(select(outlets_all, -type), by = "COMID") %>%
       select(-TotDASqKM) %>%
       st_as_sf()