05_hyRefactor_flines.Rmd

---
title: "NHD Navigate"
output: html_document
editor_options:
  chunk_output_type: console
---

Project: GFv2.0  
Script purpose: Run hyRefactor workflow on flowline network
Date: 2/29/2019  
Author: [Dave Blodgett](dblodgett@usgs.gov)
Notes: 

```{r setup_0, echo=FALSE, cache=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=6, 
  fig.height=4,
  cache=FALSE)
```

Load functions and configuration. 
See `hyRefactory_config.R` for all constants.
```{r setup}
source("R/utils.R")
source("R/hyRefactor_funs.R")
source("R/config.R")
```

Load and filter source NHD Flowline network. Bring in POIs
```{r flowlines}
nhd <- read_sf(out_refac_gpkg, nhd_flowline)

cats <- read_sf(out_refac_gpkg, nhd_catchment)
  
POIs_data <- read_sf(nav_gpkg, poi_data_table)

POIs <- read_sf(nav_gpkg, poi_geometry_table)

POI_id <- st_drop_geometry(POIs)
```

Read POIs and add nhd outlets. Save to a non-spatial table.
```{r refactor}
POIs_ref <- POIs %>%
  inner_join(select(st_drop_geometry(nhd), TotDASqKM, COMID, DnHydroseq), by = c("hy_id" = "COMID"))

# Also need to avoid modification to flowlines immediately downstream of POIs
#      This can cause some hydrologically-incorrect catchment aggregation
POI_downstream <- filter(nhd, Hydroseq %in% POIs_ref$DnHydroseq, AreaSqKM > 0)

# build final outlets set
term_poi <- filter(POIs_data, hl_reference == "Type_Term")

outlets <- POIs_ref %>%
  mutate(type = ifelse(poi_id %in% term_poi$poi_id, "terminal", "outlet")) 

# derive list of unique terminal paths
TerminalPaths <- unique(nhd$TerminalPa)

# Need to avoid refactoring catchments that are long and thing and 
# reasonably large. They cause problems in a downstream step.
avoid <- dplyr::filter(nhd, (sqrt(AreaSqKM) / LENGTHKM) > 3 & AreaSqKM > 1)

# attribute flowline network used for refactor
nhd <- mutate(nhd, refactor = ifelse(TerminalPa %in% TerminalPaths, 1, 0))
write_sf(nhd, out_refac_gpkg, nhd_flowline)

# Prep flowlines for refactor
nhdplus_flines <- sf::st_zm(filter(nhd, refactor == 1)) %>%
  st_as_sf() 

if(!needs_layer(nav_gpkg, event_geometry_table)){
  events <- read_sf(nav_gpkg, event_geometry_table) %>%
    mutate(event_identifier = as.character(row_number()))
  
  events_ref <- filter(events, hy_id %in% nhdplus_flines$COMID) %>%
    rename(COMID = hy_id) %>%
    distinct(COMID, REACHCODE, REACH_meas, event_identifier)
  
  outlets <- filter(outlets, !poi_id %in% events$poi_id) %>%
    rename(COMID = hy_id)
} else {
  events_ref <- NULL
}

if(needs_layer(out_refac_gpkg, refactored_layer)) {
  
  tf <- paste0("temp/refactored_", rpu_code,".gpkg")
  tr <- paste0("temp/reconciled_", rpu_code, ".gpkg")
    
  unlink(list(tf, tr))
  # see hyRefactor_funs.R for these parameters.
  # Select controls whether prepare_nhdplus is used. 
  refactor_nhdplus(nhdplus_flines = dplyr::select(nhdplus_flines, -FTYPE), 
                   split_flines_meters = split_meters, 
                   split_flines_cores = para_split_flines, 
                   collapse_flines_meters = combine_meters,
                   collapse_flines_main_meters = combine_meters,
                   out_refactored = tf, 
                   out_reconciled = tr, 
                   three_pass = TRUE, 
                   purge_non_dendritic = FALSE, 
                   events = events_ref,
                   exclude_cats = c(outlets$COMID, avoid$COMID, POI_downstream$COMID),
                   warn = TRUE)
  
  write_sf(st_transform(read_sf(tf), crs), out_refac_gpkg, refactored_layer)
  write_sf(st_transform(read_sf(tr), crs), out_refac_gpkg, reconciled_layer)
  
  unlink(list(tf, tr))
  rm(tf, tr)
}
```

```{r lookup}
# create lookup for ref flowlines to use in the non-dendritic steps
refactor_lookup <- st_drop_geometry(read_sf(out_refac_gpkg, reconciled_layer)) %>%
  dplyr::select(ID, member_COMID) %>%
  dplyr::mutate(member_COMID = strsplit(member_COMID, ",")) %>%
  tidyr::unnest(cols = member_COMID) %>%
  dplyr::mutate(NHDPlusV2_COMID = as.integer(member_COMID)) %>% # note as.integer truncates
  dplyr::rename(reconciled_ID = ID)

if(is.character(refactor_lookup$reconciled_ID)) refactor_lookup$reconciled_ID <- as.integer(refactor_lookup$reconciled_ID)

lookup_table <- tibble::tibble(NHDPlusV2_COMID = unique(as.integer(refactor_lookup$member_COMID))) %>%
    dplyr::left_join(refactor_lookup, by = "NHDPlusV2_COMID") 

readr::write_csv(lookup_table, lookup_table_file)
write_sf(lookup_table, out_refac_gpkg, lookup_table_refactor)

if(needs_layer(out_refac_gpkg, outlets_layer)) {
  
  # Join refactored to original NHD
  refactored <- read_sf(out_refac_gpkg, refactored_layer) %>%
    select(member_COMID = COMID, Hydroseq, event_identifier, event_REACHCODE) %>%
    inner_join(select(st_drop_geometry(nhd), orig_COMID = COMID, Hydroseq), by = "Hydroseq") 
  
  if(!is.null(events_ref)){
    # Subset for events
    refactored_events <- refactored %>%
      filter(!is.na(event_REACHCODE), !is.na(event_identifier))
  
    event_outlets <- events %>%
      inner_join(st_drop_geometry(refactored_events), by = "event_identifier") %>%
      select(hy_id, event_id, poi_id, member_COMID)
    
    # subset for refactored outlets (non-events)
    refactored_outlets <- filter(refactored, !member_COMID %in% event_outlets$member_COMID)
    
    # get ref_COMId for other outlets
    outlets_ref <- outlets %>%
      left_join(select(st_drop_geometry(refactored_outlets), member_COMID, orig_COMID), 
                by = c("COMID" = "orig_COMID")) %>%
      group_by(COMID) %>%
      filter(member_COMID == max(member_COMID)) %>%
      select(hy_id = COMID, poi_id, member_COMID, type) 
    
    outlets_ref_COMID <- data.table::rbindlist(list(outlets_ref, event_outlets), fill = TRUE) %>%
      st_as_sf()
  } else {
    # get ref_COMId for other outlets
    outlets_ref_COMID <- outlets %>%
      left_join(select(st_drop_geometry(refactored), member_COMID, orig_COMID), 
                by = c("COMID" = "orig_COMID")) %>%
      group_by(COMID) %>%
      filter(member_COMID == max(member_COMID)) %>%
      select(hy_id = COMID, poi_id, member_COMID, type) 
  }

  final_outlets <-  outlets_ref_COMID %>%
    st_as_sf() %>%
    inner_join(select(lookup_table, member_COMID, reconciled_ID), 
                      by = "member_COMID") 
  
  write_sf(final_outlets, out_refac_gpkg, outlets_layer)
} else {
  final_outlets <- read_sf(out_refac_gpkg, outlets_layer)
}

check_dups_poi <- final_outlets %>%
  group_by(reconciled_ID) %>%
  filter(n() > 1) %>%
  ungroup()

if(nrow(check_dups_poi) > 1){
  print("Double-check for double POIs")
  write_sf(check_dups_poi, out_refac_gpkg, paste0(dup_pois, "_", rpu_code))
} else {
  print("no double POIs detected")
}
```