diff --git a/workspace/06-2_aggregate_cats.Rmd b/workspace/06-2_aggregate_cats.Rmd index e02546d8485251086ce4c55dbc8fdab8321a43bf..0e239b8c574188e898b9a0895eab3d28d6d47af5 100644 --- a/workspace/06-2_aggregate_cats.Rmd +++ b/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()