Major mods see MR for details

workspace/06-2_aggregate_cats.Rmd

@@ -24,6 +24,9 @@ source("R/utils.R")
 source("R/config.R")
 source("R/NHD_navigate.R")
 source("R/hyRefactor_funs.R")
+
+old_path <- Sys.getenv("PATH") 
+Sys.setenv(PATH = paste(old_path, "C:\\Users\\abock\\AppData\\Roaming\\npm", sep = ";"))
 ```
 
 Derive set of outlets to aggregate catchments for. There are three sets of 
@@ -37,39 +40,49 @@ nhd <- read_sf(out_refac_gpkg, nhd_flowline)
 
 # Load outlets
 outlets_POI <- read_sf(out_refac_gpkg, outlets_layer) %>%
-  distinct()
+  distinct() %>%
+  filter(hy_id %in% nhd$COMID)
+
+POIs_data <- read_sf(nav_gpkg, poi_data_table)
 
 # reconciled have no ID applied in refactoring flowlines
 reconciled <- st_transform(read_sf(out_refac_gpkg, reconciled_layer), crs) 
 divides <- st_transform(read_sf(out_refac_gpkg, divide_layer), crs) 
+
+lookup <- read_sf(out_refac_gpkg, lookup_table_refactor) %>%
+  inner_join(select(st_drop_geometry(nhd), COMID, TotDASqKM), by = c("NHDPlusV2_COMID" = "COMID"))
+
+highest_DA <- lookup %>%
+  group_by(reconciled_ID) %>%
+  filter(TotDASqKM == max(TotDASqKM))
 ```
 
 ```{r agg_catchments I}
 if(needs_layer(out_agg_gpkg, agg_cats_layer)){
  
   # https://github.com/r-spatial/sf/issues/1094#issuecomment-988933580 for why this is here.
-  sf::st_geometry(divides[sf::st_is_empty(divides),]) <- 
+  sf::st_geometry(divides[sf::st_is_empty(divides),]) <-
     sf::st_cast(sf::st_geometry(divides[sf::st_is_empty(divides),]), "MULTIPOLYGON")
-  
-  #divides <- hyRefactor::clean_geometry(divides)  
+
+  #divides <- hyRefactor::clean_geometry(divides)
 
   # POI outlets
   outlets <- select(st_drop_geometry(outlets_POI), ID = reconciled_ID, type)
 
   # Identify reconciled flowpaths that are terminal and above DA threshold
-  reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID)) 
-  
+  reconciled_DA <- filter(reconciled, TotDASqKM > aggregate_da_thresh_sqkm, is.na(toID))
+
   # Ones not in existing pois
   missing_outlets <- filter(reconciled_DA, !ID %in% outlets$ID) %>%
     pull(ID)
-  
+
   # bind missing outlets, re-assign mis-attributed terminals
   if(length(missing_outlets) > 0){
     outlets <- bind_rows(outlets,
                          data.frame(ID = missing_outlets,
                                     type = rep("terminal", length(missing_outlets))))
     write_sf(filter(reconciled, ID %in% missing_outlets), out_agg_gpkg, "missing_outlets")
-    
+
     outlets <- outlets %>%
       left_join(select(st_drop_geometry(reconciled), ID, toID), by = "ID") %>%
       mutate(type = ifelse(type == "terminal" & !is.na(toID), "outlet", type)) %>%
@@ -83,10 +96,10 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
                                    da_thresh = aggregate_da_thresh_sqkm,
                                    only_larger = TRUE,
                                    post_mortem_file = cache_split)
-    
+
   agg_cats$cat_sets$set <- sapply(agg_cats$cat_sets$set, paste, collapse = ",")
   agg_cats$fline_sets$set <- sapply(agg_cats$fline_sets$set, paste, collapse = ",")
-  
+
   write_sf(agg_cats$cat_sets, out_agg_gpkg, agg_cats_layer)
   write_sf(agg_cats$fline_sets, out_agg_gpkg, agg_fline_layer)
   
@@ -96,18 +109,34 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
     st_as_sf()
   
   # Build final POI set
-  final_POIs <- st_drop_geometry(agg_cats$fline_set) %>%
+  POIs_att <- select(rec_outlets, ID, toID) %>% 
     # Bring over POI information
-    left_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
-    # Bring over physical end node XY
-    left_join(select(rec_outlets, ID), by = "ID") %>%
-    # Mark additional POI/outlets created during aggregate cats
-    mutate(Type_Con = ifelse(!ID %in% outlets$ID, 1, 0)) %>%
+    inner_join(st_drop_geometry(outlets_POI), by = c("ID" = "reconciled_ID")) %>%
+    select(-c(type, member_COMID))
+  
+  POIs_missing <- filter(select(rec_outlets, -toID), !ID %in% POIs_att$ID) %>%
+    inner_join(st_drop_geometry(agg_cats$fline_sets), by = "ID") %>%
+    arrange(ID) %>%
+    mutate(poi_id = max(POIs_att$poi_id) + row_number()) %>%
+    select(ID, toID, poi_id, event_id = event_identifier, Type_Term = orig_levelpathID) %>%
+    inner_join(select(highest_DA, NHDPlusV2_COMID, reconciled_ID), by = c("ID" = "reconciled_ID")) %>%
+    select(ID, toID, hy_id = NHDPlusV2_COMID, poi_id, event_id, Type_Term) %>%
+    mutate(Type_Term = ifelse(is.na(toID), Type_Term, NA),
+           Type_Con = ifelse(!is.na(toID), 1, NA)) %>%
+    st_as_sf()
+  
+  final_POIs <- data.table::rbindlist(list(POIs_att, select(POIs_missing, 
+                                                            -c(Type_Term, Type_Con))), fill = TRUE) %>%
     st_as_sf()
     
+  
+  unaggregated_outlets <- filter(outlets_POI, !reconciled_ID %in% final_POIs$ID)
+  double_outlets <- final_POIs %>% group_by(ID) %>% filter(n() > 1)
+  
   # Write out
   write_sf(final_POIs, out_agg_gpkg, mapped_outlets_layer)
-
+  write_sf(unaggregated_outlets, out_agg_gpkg, "unmapped_outlets")
+  write_sf(double_outlets, out_agg_gpkg, "double_outlets")
 } else {
   agg_cats <- list(cat_sets = read_sf(out_agg_gpkg, agg_cats_layer),
                    fline_sets = read_sf(out_agg_gpkg, agg_fline_layer))
@@ -117,6 +146,25 @@ if(needs_layer(out_agg_gpkg, agg_cats_layer)){
 
 ```
 
+```{r Long form POIs}
+POIs <- final_POIs %>%
+  arrange(ID) 
+
+final_POI_geom <- POIs %>%
+  select(ID) %>%
+  cbind(st_coordinates(.)) %>%
+  group_by(ID, X, Y) %>%
+  mutate(geom_id = cur_group_id()) %>%
+  ungroup()
+
+final_POIs_table <- POIs %>%
+  inner_join(select(st_drop_geometry(final_POI_geom), -X, -Y), by = "ID") %>%
+  distinct()
+
+write_sf(final_POIs_table, out_agg_gpkg, poi_data_table)
+```
+
+
 ```{r lookup table}
 refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID) %>%
   dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%