utils.R

hyfabric_source <- list.files(pattern = "hyfabric_.*.tar.gz")
hyfabric_version <- substr(hyfabric_source, 10, 14)

if(!(require(hyfabric, quietly = TRUE) && packageVersion("hyfabric") == hyfabric_version)) {
  dir.create(path = Sys.getenv("R_LIBS_USER"), showWarnings = FALSE, recursive = TRUE)
  
  .libPaths(Sys.getenv("R_LIBS_USER"))
  
  install.packages(hyfabric_source, lib = Sys.getenv("R_LIBS_USER"),
                    repos = NULL, type = "source")
}

fix_headwaters <- function(nhd_rds, out_gpkg, new_atts = NULL, 
                           nhdpath = nhdplusTools::nhdplus_path(), 
                           par = 6) {
  
  if(!file.exists(out_gpkg)) {
  nhd <- readRDS(nhd_rds)
  
  # Need to replace headwater flowlines with the trimmed burn line event.
  # Some headwater flowlines are not fully within the catchment they go with.
  ble <- sf::read_sf(nhdpath, "BurnLineEvent")
  
  ble <- dplyr::left_join(dplyr::select(sf::st_drop_geometry(nhd), COMID), 
                          ble, by = "COMID")
  
  ble <- sf::st_zm(sf::st_as_sf(ble))
  nhd <- sf::st_zm(nhd)
  
  sf::st_geometry(nhd)[nhd$StartFlag == 1] <- sf::st_geometry(ble)[nhd$StartFlag == 1]
  
  if(!is.null(new_atts)) {
    new_atts <- data.table::fread(new_atts, 
                                  integer64 = "character", 
                                  showProgress = FALSE)
    
    if("weight" %in% names(new_atts)) new_atts <- rename(new_atts, 
                                                        ArbolateSu = weight)
    
    nhd <- select(nhd, COMID, VPUID, RPUID, FTYPE, FromNode, ToNode, 
                  StartFlag, StreamCalc, Divergence, WBAREACOMI, 
                  VPUID, RPUID, DnMinorHyd) %>%
      right_join(new_atts, by = c("COMID" = "comid")) %>%
      nhdplusTools::align_nhdplus_names()
  }
  
  nhd <- sf::st_sf(nhd)
  
  m_to_km <- (1/1000)
  
  nhd$LENGTHKM <- as.numeric(sf::st_length(sf::st_transform(nhd, 5070))) * m_to_km

  # is a headwater and for sure flows to something.  
  check <- !nhd$COMID %in% nhd$toCOMID & !(nhd$toCOMID == 0 | is.na(nhd$toCOMID) | !nhd$toCOMID %in% nhd$COMID)
  check_direction <- dplyr::filter(nhd, check)$COMID
  
  cl <- parallel::makeCluster(16)
  
  # check and fix these.
  new_geom <- pbapply::pblapply(check_direction, fix_flowdir, 
                                network = nhd[nhd$COMID %in% check_direction, ], cl = cl)
  
  parallel::stopCluster(cl)
  
  # One is empty. Deal with it and remove from the replacement.
  error_index <- sapply(new_geom, inherits, what = "try-error")
  errors <- filter(nhd, COMID %in% check_direction[error_index])
  check[which(nhd$COMID %in% errors$COMID)] <- FALSE
  
  if(!all(sapply(sf::st_geometry(errors), sf::st_is_empty))) {
    stop("Errors other than empty geometry from fix flowdir")
  }
  
  new_geom <- new_geom[!error_index]
  new_geom <- do.call(c, new_geom)
  st_geometry(nhd)[check] <- new_geom
  
  sf::write_sf(nhd, out_gpkg, "reference_flowlines")
  
  }
  return(out_gpkg)
}

## Is available in sbtools now, but will leave till on cran.
#' Retrieves files from SB in facet file_structure
#' @param sbitem (character) ScienceBase item ID
#' @param out_dir (directory) directory where files are downloaded
#' 
# Retrieve files from SB if files not listed under sbtools
retrieve_sb <- function(sbitem, out_dir){
  
  item <- item_get(sbitem)
  
  linked_files <- 
    lapply(item$facets, function(x) 
    {
      list(name = x$name, 
           files = lapply(x$files, function(y, out_path)
           {
             
             out_file <- file.path(out_path, y$name)
             download.file(y$downloadUri, out_file, mode = "wb")
             
             list(fname = y$name,
                  url = y$downloadUri)
             
           }, 
           out_path = file.path(data_dir, x$name)))
    })
}


#' Documents why streamgage is excluded from POIs
#' @param source (character) data source
#' @param reason (character) reason gage is excluded
#' @return writes output to file
#' 
gage_document <- function(data, source, drop_reason){
  if(!file.exists("cache/dropped_gages.csv")){
    write.csv(data %>% mutate(data_source = source, reason = drop_reason),
              "cache/dropped_gages.csv", quote = F, row.names = F)
  } else {
    write.table(data %>%
                  mutate(data_source = source, reason = drop_reason),
                file = "cache/dropped_gages.csv", append = T, sep = ",", 
                col.names = F, row.names = F, quote = F)
  }
}

#' clean network
#' @param net data.frame nhdplus network
#' @param net data.frame new nhdplus network as from e2nhd
#' @param nwm data.frame with national water model new network
#' @return data.frame fully reconciled attributes
clean_net <- function(net, net_new, nwm) {
  
  # NOTE: need to disconnected diverted paths as indicated in the tonode / fromnode attributes.
  
  diverted_heads <- filter(net_new, divergence == 2)
  
  net_new <- filter(net_new, !tocomid %in% diverted_heads$comid)
  
  # NOTE: avoid fcode == 56600 (coastal)
  avoid <- net$comid[net$fcode == 56600]
  
  # NOTE: avoid large groups where things have been redirected in an un-natural way.
  groups <- group_by(nwm, tocomid)
  
  groups <- data.frame(tocomid = summarise(groups, .groups = "drop"), size = group_size(groups))
  
  groups <- filter(groups, size > 10 & tocomid != 0)
  
  # NOTE: also avoid all flowlines that go to coastal flowlines.
  # This is due to unwanted great lakes connections in the NWM.
  avoid <- c(avoid, 
             net_new$comid[net_new$tocomid %in% avoid], 
             net_new$comid[net_new$tocomid %in% groups$tocomid],
             nwm$comid[nwm$tocomid %in% groups$tocomid])
  
  avoid <- unique(avoid)
  
  nwm <- filter(nwm, tocomid %in% c(net_new$comid, 0))
  
  net_new <- left_join(net_new, select(nwm, comid, 
                                       nwm_tocomid = tocomid), 
                       by = "comid") %>% 
    mutate(tocomid = ifelse(is.na(tocomid), 0, tocomid))
  ## NWM already uses 0 for outlets.
  
  # NOTE: just being explicit about this for clarity.
  candidate <- select(net_new, comid, tocomid, nwm_tocomid) %>%
    filter(!comid %in% avoid)
  
  mismatch <- candidate[candidate$tocomid != candidate$nwm_tocomid, ]
  
  net_new <- left_join(net_new, select(mismatch, comid, 
                                       new_tocomid = nwm_tocomid), 
                       by = "comid")
  
  to_change <- filter(net_new, !is.na(new_tocomid))
  
  net_new <- net_new %>%
    mutate(divergence = ifelse(comid %in% to_change$new_tocomid & divergence == 2, 1,
                               ifelse(comid %in% to_change$tocomid & divergence == 1, 2, 
                                      divergence)),
           tocomid = ifelse(comid %in% to_change$comid, new_tocomid, tocomid)) %>%
    select(-nwm_tocomid, -new_tocomid)
  
  net_new <- filter(net_new, fcode != 56600)
  
  net_new <- mutate(net_new, tocomid = ifelse(!tocomid %in% comid, 0, tocomid))
  
  return(net_new)
  
}

#' Merges geopackages together to create CONUs geopackage of features
#' @param feat (character) Type of feature to merge (POIs, segments) character
#' @param in_gpkg  (character) geopackage to merge (GF, ND_, etc.)
#' @param out_gpkg (character) Name of output geopackage
#' @return output geopackage of CONUS for given features
#'
Merge_VPU <- function(feat, in_gpkg, out_gpkg){

  if(needs_layer(out_gpkg, feat)) {
    all_sf <- paste0(feat, "_CONUS")
    #VPUs <- c("03N", "03S", "03W")
    #VPUs <- c("10L", "10U")
    VPUs <-c(paste0("0", c(1:2)), "03N", "03S", "03W", paste0("0", c(4:9)),
             as.character(11:18), "10U", "10L")
    featCON <- do.call("rbind", lapply(VPUs, function(x) {
      tryCatch({
        layer <- ifelse(feat %in% c("nhd_flowline", "unassigned_gages", "unassigned_TE"), 
                                    feat, paste0(feat, "_", x))
        read_sf(file.path("cache", paste0(in_gpkg, x, ".gpkg")), 
                layer)}, 
        error = function(e) stop(paste("didn't find", x,
                                       "\n Original error was", e)))
    }))
    write_sf(featCON, out_gpkg, feat)
  }
}

#' Capture all catchments (flowline, sink, etc.) in a given RPU
#' @param fcats (sf data.frame) NHDPlusv2.1 National geodatabase catchment layer
#' @param the_RPU (character) RPU to retrive full catchment set.
#' @param fl_rds character path to flowline rds file
#' @param nhd_gdb character path to NHD geodatabase
#' @return (sf data.frame) with FEATUREID, RPUID, FTYPE, TERMINAFL 
cat_rpu <- function(fcats, nhd_gdb){
  
  fl <- read_sf(nhd_gdb, "NHDFlowline_Network") %>%
    align_nhdplus_names()
  
  # read the CatchmentSP
  cats <- readRDS(fcats) %>%
    st_as_sf() %>%
    st_drop_geometry() %>%
    dplyr::select(FEATUREID, SOURCEFC)
  
  # read the Flowlines
  flowln_df <- fl %>%
    st_as_sf() %>%
    st_drop_geometry() %>%
    select(COMID, FTYPE, TerminalFl, RPUID)%>%
    rename(FEATUREID = COMID) %>%
    rename(TERMINALFL = TerminalFl)%>%
    filter(FEATUREID %in% cats$FEATUREID)
  
  # Read in sinks
  sink_df <- sf::st_read(nhd_gdb, layer = "Sink") %>%
    st_drop_geometry() %>%
    dplyr::select(SINKID, RPUID = InRPU) %>%
    dplyr::rename(FEATUREID = SINKID) %>%
    dplyr::mutate(FTYPE = "Sink", TERMINALFL = 0) %>%
    dplyr::filter(!FEATUREID %in% flowln_df$FEATUREID)
  
  #combine all the RPUIDs
  lkup_rpu <- rbind(flowln_df, sink_df)
  
  # FEATUREID 10957920, 20131674, 24534636 - this are the ids of the missing, checked on map,
  #     they look like they should not be used

  # add the records for the missing
  missrec_df <- data.frame(FEATUREID=c(10957920, 20131674, 245346360), 
                           RPUID = c("03a", "03d", "17b"),
                           TERMINALFL = c(0, 0, 0)) %>%
    mutate(FTYPE = "")
  
  lkup_rpu2 <- rbind(lkup_rpu, missrec_df)
  return(lkup_rpu2)
}