Renumbered

workspace/03_hyRefactor_AK.Rmd

+---
+title: "hyRefactor_AK.Rmd"
+output: html_document
+editor_options:
+  chunk_output_type: console
+---
+
+Project: GFv2.0  
+Script purpose: Run hyRefactor workflow on flowline network for merrit flowlines/catchments for AK
+Date: 2/29/2019  
+Author: [Dave Blodgett](dblodgett@usgs.gov)
+Notes: mod for AK, nov 2020, by MS -- debug by dblodgett 2/2022
+
+```{r setup_0, echo=FALSE, cache=FALSE}
+knitr::opts_chunk$set(
+  collapse = TRUE,
+  comment = "#>",
+  fig.width=6, 
+  fig.height=4)
+```
+
+Load functions and configuration. 
+See `hyRefactory_config.R` for all constants.
+```{r setup}
+library(hyfabric)
+rpu_code <- "19a"
+VPU <- substr(rpu_code, 1, 2)
+out_agg_gpkg <- file.path("cache", paste0(rpu_code, ".gpkg"))
+
+source("R/utils.R")
+source("R/hyRefactor_funs.R")
+source("R/config.R")
+source("R/ExCONUS_config.R")
+
+data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
+```
+
+Load and filter source data. This section cleanly sets up the environment for later.
+```{r network}
+if(needs_layer(out_refac_gpkg, ak_cats_layer)) {
+  
+  cats <- read_sf(nav_gpkg, ak_cats_layer)
+
+  write_sf(cats, out_refac_gpkg, ak_cats_layer)
+} else {
+  cats <- read_sf(out_refac_gpkg, ak_cats_layer)
+}
+
+# from reference fabric
+nhd <- read_sf(nav_gpkg, nhd_flowline)
+
+# Thematic POIs built in Navigate
+POIs <- read_sf(nav_gpkg, nav_poi_layer) %>%
+  inner_join(select(st_drop_geometry(nhd), COMID, DnHydroseq), 
+             by = "COMID")
+
+# build final outlets set
+outlets <- POIs %>%
+  mutate(type = ifelse(!is.na(Type_Term), "terminal", "outlet")) %>%
+  filter(COMID %in% cats$COMID)
+
+if(needs_layer(out_refac_gpkg, outlets_layer)) {
+  write_sf(outlets, out_refac_gpkg, outlets_layer)
+}
+
+if(needs_layer(out_refac_gpkg, nhd_flowline)) {
+  
+  # attribute flowline network used for refactor
+  nhd <- mutate(nhd, refactor = ifelse(
+    TerminalPa %in% unique(filter(nhd, COMID %in% outlets$COMID)$TerminalPa), 
+    1, 0)) %>% # Get Events for catchment splitting
+    mutate(REACHCODE = COMID, FromMeas = 0, ToMeas = 100)
+  
+  write_sf(nhd, out_refac_gpkg, nhd_flowline)
+  
+} else {
+  
+  nhd <- read_sf(out_refac_gpkg, nhd_flowline)
+  
+}
+
+POI_downstream <- filter(nhd, Hydroseq %in% POIs$DnHydroseq, AreaSqKM > 0)$COMID
+avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)$COMID
+
+events <- read_sf(nav_gpkg, events_DF_layer) %>%
+  left_join(select(st_drop_geometry(nhd), 
+                    COMID, LENGTHKM, FromMeas, ToMeas), by = "COMID") %>%
+  filter(REACH_meas < 90 & REACH_meas > 10 & LENGTHKM > (combine_meters / 1000)) %>%
+  # Events cannot be in terminal POIs, code seems to ignore the command not to 
+  # split/combine those flowlines
+  filter(!COMID %in% filter(outlets, type == "terminal")$COMID) %>%
+  filter(!COMID %in% avoid)
+
+if(needs_layer(out_refac_gpkg, split_layer)) {
+  write_sf(events, out_refac_gpkg, split_layer)
+}
+
+```
+
+Run refactor / reconcile
+```{r refactor}
+if(needs_layer(out_refac_gpkg, reconciled_layer)){
+
+  tf <- paste0("temp/refactored_", rpu_code,".gpkg")
+  tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
+  
+  # see hyRefactor_funs.R for these parameters.
+  refactor_nhdplus(nhdplus_flines = nhd, 
+                   split_flines_meters = 10000000, 
+                   split_flines_cores = 6, 
+                   collapse_flines_meters = combine_meters,
+                   collapse_flines_main_meters = combine_meters,
+                   out_refactored = tf, 
+                   out_reconciled = tr, 
+                   three_pass = TRUE, 
+                   purge_non_dendritic = FALSE, 
+                   exclude_cats = c(avoid, POI_downstream),
+                   events = events,
+                   warn = TRUE)
+  
+  write_sf(read_sf(tf), out_refac_gpkg, refactored_layer)
+  write_sf(read_sf(tr), out_refac_gpkg, reconciled_layer)
+  
+  unlink(list(tf, tr))
+  rm(tf, tr)
+} else {
+  # reconciled have no ID applied in refactoring flowlines
+  reconciled <- read_sf(out_refac_gpkg, reconciled_layer)
+  # refactored has the original COMIDs and other attributes
+  refactored <- read_sf(out_refac_gpkg, refactored_layer)
+}
+```
+
+```{r lookup}
+# create lookup for ref flowlines to use in the non-dendritic steps
+lookup_table <- dplyr::select(st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)), ID, member_COMID) %>%
+  dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
+  tidyr::unnest(cols = member_COMID) %>%
+  dplyr::mutate(merit_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
+  dplyr::rename(reconciled_ID = ID)
+
+if(is.character(lookup_table$reconciled_ID)) lookup_table$reconciled_ID <- as.integer(lookup_table$reconciled_ID)
+
+lookup_table <- tibble::tibble(merit_COMID = unique(as.integer(lookup_table$member_COMID))) %>%
+    dplyr::left_join(lookup_table, by = "merit_COMID") 
+
+
+readr::write_csv(lookup_table, lookup_table_file)
+write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
+```
+
+
+For Alaska, we won't split catchments for now. So no need for elevation data.
+```{r}
+if(needs_layer(out_refac_gpkg, divide_layer)) {
+  # reconciled have no ID applied in refactoring flowlines
+  reconciled <- read_sf(out_refac_gpkg, reconciled_layer)
+  # refactored has the original COMIDs and other attributes
+  refactored <- read_sf(out_refac_gpkg, refactored_layer)
+
+  fdr_temp <- terra::rast(data_paths$merit_fdr)
+  r_crs <- terra::crs(fdr_temp)
+  
+  refactored <- sf::st_transform(refactored, r_crs)
+  reconciled <- sf::st_transform(reconciled, r_crs)
+  
+  divides <- st_transform(cats, r_crs) %>%
+    rename(FEATUREID = COMID) %>%
+    reconcile_catchment_divides(fline_ref = refactored,
+                                fline_rec = reconciled,
+                                fdr = data_paths$merit_fdr,
+                                fac = data_paths$merit_fac,
+                                para = 6,
+                                cache = cache_split,
+                                vector_crs = crs, 
+                                keep = NULL, fix_catchments = TRUE)
+
+  if(exists("divides")) unlink(cache_split)
+  
+  divides <- sf::st_transform(divides, crs) %>%
+    rename_geometry("geometry")
+    
+  write_sf(divides, out_refac_gpkg, divide_layer)
+} else {
+  reconciled <- read_sf(out_refac_gpkg, reconciled_layer)
+  refactored <- read_sf(out_refac_gpkg, refactored_layer)
+  divides <- read_sf(out_refac_gpkg, divide_layer)
+}
+```
+
+This next section maps outlets to the approprirate reconciled flowline.  
+There are three classes of outlets:
+1) Initial POI set (derived using hyRefactor::map_outlet_ids); lines 223 - 232
+2) Events which were on the upstream portion of a split flowline; 
+  lines 234 - 244
+3) POIs downstream of an event, which need to remapped to the downstream 
+   portion of a split flowline; lines 247-261
+   
+If no events, the section skips to line 282
+
+```{r POIs/outlets}
+# Read lookup for original COMID
+lookup_table <- read_sf(out_refac_gpkg, lookup_table_refactor)
+
+# Create unified POI/outlet layer
+if(needs_layer(out_agg_gpkg, mapped_outlets_layer)) {
+
+  # Get the end of the reconciled flowlines and bind with data frame
+  reconciled_end <- get_node(reconciled, "end") %>%
+    cbind(st_drop_geometry(reconciled)) %>%
+    # Join with divides to get the total drainage area
+    left_join(select(st_drop_geometry(divides), ID, area = areasqkm), by = "ID") %>%
+    mutate(area = ifelse(is.na(area), 0, area)) %>%
+    mutate(TotDASqKM = calculate_total_drainage_area(select(., ID, toID, area)))
+  
+  # If there are events on split flowlines
+  if(!needs_layer(out_refac_gpkg, split_layer)){
+    # Read in events
+    events <- read_sf(out_refac_gpkg, split_layer) 
+    
+    # Get rows of lookup and reconciled flowlines for split events
+    lookup_events <- filter(lookup_table, merit_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, !COMID %in% events$COMID) %>%
+      # Reset Type_Gages attrbute to NA if event exists
+      distinct(.keep_all = TRUE) 
+    
+    # Map non-event outlets ID to reconciled ID 
+    poi_outlets_mapped <- hyRefactor::map_outlet_ids(select(st_drop_geometry(poi_outlets), COMID, type), reconciled) %>%
+      mutate(COMID = as.integer(COMID)) %>%
+      inner_join(select(outlets, -type), by = "COMID") %>%
+      st_as_sf() 
+
+    # Generate outlet ID data frame for split events
+    events_mapped <- select(st_drop_geometry(events), COMID, POI_identifier, type) %>%
+      #IDs that events are associated with
+      mutate(ID = reconciled_end[st_nearest_feature(st_transform(events, st_crs(reconciled_end)), 
+                                                 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") %>%
+      st_as_sf() %>%
+      unique() %>%
+      select(COMID, ID, type, POI_identifier) 
+    
+    # Pois on downstream flowlines split by events
+    poi_DS_events_mapped <- filter(outlets, COMID %in% events$COMID) %>%
+      distinct() %>%
+      # Subset to POIs that have other POI flags than Type_Gage
+      filter_at(vars(-c(COMID, type, geom)), any_vars(!is.na(.))) %>%
+      left_join(lookup_table, by = c("COMID" = "merit_COMID")) %>%
+      inner_join(select(st_drop_geometry(reconciled_end), ID, TotDASqKM, member_COMID), 
+                 by = "member_COMID") %>% 
+      group_by(COMID) %>% 
+      # Max DA associated with each orig_COMID will be the proper reconciled ID value
+      filter(TotDASqKM == max(TotDASqKM)) %>%
+      select(-TotDASqKM) %>%
+      st_compatibalize(., events_mapped) 
+    
+    # Bind together three outlet categories above
+    outlets_POI <- data.table::rbindlist(list(poi_outlets_mapped, events_mapped, poi_DS_events_mapped),
+                                         fill = TRUE) %>%
+      st_as_sf()
+    
+    # combine duplicates if they exist and put back in 
+    if(nrow(outlets_POI) != length(unique(outlets_POI$ID))){
+      duplicates <- outlets_POI %>%
+        group_by(ID) %>%
+        filter(n() > 1) %>%
+        mutate_all(funs(zoo::na.locf(., na.rm = FALSE, fromLast = FALSE))) %>%
+        filter(row_number() == n()) %>%
+        ungroup()
+      
+      outlets_POI <- filter(outlets_POI, !ID %in% duplicates$ID) %>%
+        rbind(duplicates)
+    }
+
+    
+  } else {
+    
+    # Map non-event outlets ID to reconciled ID 
+    outlets_POI<- map_outlet_ids(select(outlets, COMID, type), reconciled) %>%
+      inner_join(mutate(select(outlets, -type), COMID = as.character(COMID)),
+                 by = "COMID") %>%
+      st_as_sf()
+    
+  }
+
+  # Write out combined POIs table
+  write_sf(outlets_POI, out_agg_gpkg, mapped_outlets_layer)
+  
+} else {
+  outlets_POI <- read_sf(out_agg_gpkg, mapped_outlets_layer)
+}
+
+sf::st_layers(out_refac_gpkg)
+```
+
+```{r}
+if(needs_layer(out_agg_gpkg, agg_cats_layer)) {
+  
+  outlets <- dplyr::select(st_drop_geometry(outlets_POI), ID, type)
+  
+  #missing_outlets_rec <- reconciled$ID[is.na(reconciled$toID) &
+  #                                   !(reconciled$ID %in% outlets$ID[outlets$type == "terminal"])]
+
+  #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")
+  #}
+  
+  agg_cats <- hyRefactor::aggregate_catchments(flowpath = reconciled, 
+                                               divide = divides, 
+                                               outlets = outlets,
+                                               da_thresh = 1, 
+                                               only_larger = FALSE, 
+                                               # could be TRUE but geometry issues causes problems here.
+                                               # should cross check what outlets get added with what POIs 
+                                               # get passed in.
+                                               keep = NULL)
+  
+  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 = ",")
+    
+  missing_outlets_agg <- filter(st_drop_geometry(agg_cats$fline_sets), !ID %in% outlets$ID) %>%
+    cbind(get_node(filter(agg_cats$fline_sets, !ID %in% outlets$ID))) %>%
+    st_as_sf() %>%
+    st_compatibalize(., outlets_POI)
+  
+  # Get physical geometry  of reconciled FL end node whose ID is an aggregated outlet
+  rec_outlets <- filter(st_drop_geometry(reconciled), ID %in% agg_cats$fline_sets$ID) %>%
+    cbind(get_node(filter(reconciled, ID %in% agg_cats$fline_sets$ID))) %>%
+    st_as_sf()
+
+  outlets_POI <- data.table::rbindlist(list(outlets_POI, missing_outlets_agg), fill = T) %>%
+    st_as_sf()
+  
+  # Build final POI set
+  final_POIs <- st_drop_geometry(agg_cats$fline_set) %>%
+    # Bring over POI information
+    left_join(st_drop_geometry(outlets_POI), by = "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)) %>%
+    st_as_sf()
+    
+  # Write out
+  write_sf(final_POIs, out_agg_gpkg, mapped_outlets_layer)
+  
+  write_sf(agg_cats$cat_sets, out_agg_gpkg, agg_cats_layer)
+  write_sf(agg_cats$fline_sets, out_agg_gpkg, agg_fline_layer)
+} else {
+  agg_cats <- list(cat_sets = read_sf(out_agg_gpkg, agg_cats_layer),
+                   fline_sets = read_sf(out_agg_gpkg, agg_fline_layer))
+}
+
+agg_cats$cats_sets <- read_sf(out_agg_gpkg, agg_cats_layer)
+agg_cats$cat_sets$set <- lapply(strsplit(agg_cats$cat_sets$set, split = ","), as.integer)
+agg_cats$fline_sets <- read_sf(out_agg_gpkg, agg_fline_layer)
+agg_cats$fline_sets$set <- lapply(strsplit(agg_cats$fline_sets$set, split = ","), as.integer)
+
+```
+
+```{r}
+refactor_lookup <- dplyr::select(st_drop_geometry(reconciled), ID, member_COMID) %>%
+  dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
+  tidyr::unnest(cols = member_COMID) %>%
+  dplyr::mutate(merit_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
+  dplyr::rename(reconciled_ID = ID)
+
+aggregate_lookup_fline <- dplyr::select(st_drop_geometry(agg_cats$fline_sets), ID, set) %>%
+  tidyr::unnest(cols = set) %>%
+  dplyr::rename(aggregated_flowpath_ID = ID, reconciled_ID = set)
+
+aggregate_lookup_catchment <- dplyr::select(st_drop_geometry(agg_cats$cat_sets), ID, set) %>%
+  tidyr::unnest(cols = set) %>%
+  dplyr::rename(aggregated_divide_ID = ID, reconciled_ID = set)
+
+if(is.character(aggregate_lookup_catchment$reconciled_ID)) 
+  aggregate_lookup_catchment$reconciled_ID <- as.integer(aggregate_lookup_catchment$reconciled_ID)
+if(is.character(aggregate_lookup_fline$reconciled_ID)) 
+  aggregate_lookup_fline$reconciled_ID <- as.integer(aggregate_lookup_fline$reconciled_ID)
+
+lookup_table <- tibble::tibble(merit_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
+    dplyr::left_join(refactor_lookup, by = "merit_COMID") %>%
+    dplyr::left_join(aggregate_lookup_fline, by = "reconciled_ID") %>%
+    dplyr::left_join(aggregate_lookup_catchment, by = "reconciled_ID")
+
+if(!file.exists(lookup_table_file)) {
+  readr::write_csv(lookup_table, lookup_table_file)
+  write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)
+}
+```