diff --git a/workspace/03_hyRefactor_AK.Rmd b/workspace/03_hyRefactor_AK.Rmd
new file mode 100644
index 0000000000000000000000000000000000000000..9e2c113e433b4d0015c67b54750f66a994b92352
--- /dev/null
+++ b/workspace/03_hyRefactor_AK.Rmd
@@ -0,0 +1,398 @@
+---
+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)
+}
+```