01_NHD_prep.Rmd

---
title: "NHD Prep"
output: html_document
editor_options:
  chunk_output_type: console
---
Prepares the NHD network for navigate and refactor. Adds base layers to the reference fabric.

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

```{r setup}
# Load custom functions
source("R/utils.R")
source("R/NHD_navigate.R")

# Load Configuration of environment
source("R/config.R")

if(!exists("vpu_codes")) {
  vpu_codes <- unique(rpu_vpu$vpuid)
} else {
  rpu_vpu <- filter(rpu_vpu, vpuid %in% vpu_codes)
}

all_rpu_codes <- unique(rpu_vpu$rpuid)
full_cat_table <- readRDS(data_paths$fullcats_table)

out_vpu <- rep(list(list()), length(vpu_codes))
names(out_vpu) <- vpu_codes

out_rpu <- rep(list(list()), length(all_rpu_codes))
names(out_rpu) <- all_rpu_codes

fline <- sf::read_sf(data_paths$ref_flowline)

fline <- sf::st_cast(fline, "LINESTRING")

catchment <- sf::read_sf(file.path(data_paths$nhdplus_dir, "reference_catchments.gpkg"))

# we can remove truely degenerate COMIDs 
# for 0 upstream area and no catchment area
degen_comid <- fline[fline$TotDASqKM == 0 & 
                       !fline$COMID %in% catchment$featureid, ]$COMID

# need to make sure we don't disconnect anything.
keep_tocomid <- fline$toCOMID[!fline$COMID %in% degen_comid]

if(length(degen_comid[degen_comid %in% keep_tocomid]) > 0) stop("this will break the network")

fline <- fline[!fline$COMID %in% degen_comid, ]

for(VPU in vpu_codes) {
  
  rm(rpu_code)
  
  source("R/config.R")
  
  if (needs_layer(nav_gpkg, nhd_flowline)) {
    
    # Subset NHD by VPU
    nhd <- subset_vpu(fline, vpu, run_make_standalone = FALSE) 
    
    # there are some duplicates with enhd
    nhd <- nhd %>%
      group_by(COMID) %>%
      filter(row_number() == 1) %>%
      ungroup()
    
    nhd <- make_standalone(nhd)
    
    # Filter and write dendritic/non-coastal subset to gpkg
    # This will be iterated over to produce the final network after POIs identified
    zero_order <- filter(nhd, TerminalFl == 1 & TotDASqKM < min_da_km)
    
    non_dend <- unique(unlist(lapply(zero_order$COMID, NetworkNav, nhdDF = st_drop_geometry(nhd))))
    
    # Add fields to note dendritic and POI flowlines
    nhd <- nhd %>% 
      mutate(dend = ifelse(!COMID %in% non_dend, 1, 0),
             poi = ifelse(!COMID %in% non_dend & TotDASqKM >= min_da_km, 1, 0)) 
    
    write_sf(nhd, nav_gpkg, nhd_flowline)
    
    cats_rpu <- full_cat_table %>%
      filter(RPUID %in% rpu_codes$rpuid)
    
    cats <- catchment %>% 
      rename(FEATUREID = featureid, 
             AREASQKM = areasqkm, 
             RPUID = rpuid, 
             VPUID = vpuid) %>%
      filter(FEATUREID %in% 
               unique(c(nhd$COMID, 
                        full_cat_table$FEATUREID[full_cat_table$RPUID %in% rpu_codes$rpuid]))) %>%
      mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
             full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
      filter(full_cats == 1 | hy_cats == 1) %>%
      st_sf()
    
    write_sf(cats, nav_gpkg, nhd_catchment)
    
  } else {
    
    nhd <- read_sf(nav_gpkg, nhd_flowline)
    cats <- read_sf(nav_gpkg, nhd_catchment)
    
  }
  
  out_vpu[[VPU]] <- nhd %>%
    sf::st_drop_geometry() %>%
    select(COMID, toCOMID) %>%
    filter(!toCOMID %in% COMID & !toCOMID == 0)
  
  for(rpu_code in rpu_codes$rpuid) {
    source("R/config.R")
    
    if(needs_layer(out_refac_gpkg, nhd_flowline)) {
      nhd_sub <- subset_rpu(nhd, rpu_code, run_make_standalone = TRUE) %>%
        st_sf()
      
      write_sf(nhd_sub, out_refac_gpkg, nhd_flowline)
      
      cats_rpu <- full_cat_table %>%
        filter(RPUID == rpu_code)
      
      cats %>%
        mutate(hy_cats = ifelse(FEATUREID %in% nhd_sub$COMID, 1, 0),
               full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
        filter(full_cats == 1 | hy_cats == 1) %>%
        st_sf() %>%
        write_sf(out_refac_gpkg, nhd_catchment)
      
    } else {
      nhd_sub <- read_sf(out_refac_gpkg, nhd_flowline)
    }
    
    out_rpu[[rpu_code]] <- nhd_sub %>%
      sf::st_drop_geometry() %>%
      select(COMID, toCOMID) %>%
      filter(!toCOMID %in% COMID & !toCOMID == 0)
  }
}
```

```{r}

if(!file.exists("cache/rpu_vpu_out.csv")) {
  make_df <- function(x, d, n) {
    y <- d[[x]]
    nr <- nrow(y)
    na <- names(d)[x]
    o <- data.frame(d = rep(na, nr), 
                    COMID = d[[x]]$COMID, 
                    toCOMID = d[[x]]$toCOMID)
    names(o) <- c(n, "COMID", "toCOMID") 
    o
  }
  
  rpu <- do.call(rbind, lapply(1:length(out_rpu), make_df, d = out_rpu, n = "rpu"))
  
  vpu <- do.call(rbind, lapply(1:length(out_vpu), make_df, d = out_vpu, n = "vpu"))
  
  out_rpu_vpu <- left_join(rpu, vpu, by = "COMID") 
  
  out_rpu_vpu <- select(out_rpu_vpu, RPUID = rpu, VPUID = vpu, 
                        COMID = COMID, toCOMID = toCOMID.x)
  
  out_rpu_vpu <- left_join(out_rpu_vpu, select(sf::st_drop_geometry(fline), COMID, toRPUID = RPUID, toVPUID = VPUID), by = c("toCOMID" = "COMID"))
  
  readr::write_csv(out_rpu_vpu, "cache/rpu_vpu_out.csv")
}

```