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")
}

## check or re-initiate authorizatoin
#' Retrieves files from SB in facet file_structure
#' @param self
check_auth <- function(){
  if(!sbtools::is_logged_in()) {
    try(sbtools::initialize_sciencebase_session())
    if(!sbtools::is_logged_in()) stop("stale sciencebase login")
  }
  TRUE
}


#' Documents why streamgage is excluded from POIs
#' @param source (character) data source
#' @param reason (character) reason gage is excluded
#' @return writes output to file
#' 
document_dropped_gage <- 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)
  }
}


#' 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("POIs", "nhd_flowline", "unassigned_gages", "unassigned_TE",
                                    "reference_flowline", "reference_catchment"), 
                                    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()
  
  # Reference catchments vs. island catchments
  if(basename(fcats) == "reference_catchments.gpkg"){
    cats <- read_sf(fcats) %>%
      st_as_sf() %>%
      st_drop_geometry() %>%
      dplyr::select(FEATUREID = featureid)
  } else {
    cats <- read_sf(fcats, "CatchmentSP") %>%
      st_as_sf() %>%
      st_drop_geometry() 
  }
  
  # 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)
}

merge_conus <- function(rpu_vpu, rpu_vpu_out,
                        lookup_table_refactor, 
                        reconciled_layer, 
                        divide_layer, 
                        agg_fline_layer,
                        agg_cats_layer,
                        mapped_outlets_layer) {
  
  vpus <- unique(rpu_vpu$vpuid)
  
  #gf_gpkgs <- paste0("cache/refactor_", vpus, ".gpkg")
  gf_gpkgs <- paste0("cache/GF_", vpus, ".gpkg")
  
  hydrofab::assign_global_identifiers(
    gpkgs = gf_gpkgs, 
    outfiles = gsub("cache", "temp", gf_gpkgs), 
    flowpath_layer = agg_fline_layer, #reconciled_layer, 
    mapped_POI_layer = mapped_outlets_layer, 
    divide_layer = agg_cats_layer, #divide_layer, 
    lookup_table_layer = lookup_table_refactor, 
    flowpath_edge_list = catchment_network_table, 
    update_terminals = TRUE, 
    return_lookup = TRUE, 
    verbose = TRUE)
  
}

# FIXME: THESE FUNCTIONS SHOULD COME FROM hyFab
#' Extract nexus locations for Reference POIs
#' @param gpkg a reference hydrofabric gpkg
#' @param type the type of desired POIs
#' @param verbose should messages be emitted?
#' @return data.frame with ID, type columns
#' @export
nexus_from_poi = function(mapped_POIs,
                          type = c('HUC12', 'Gages', 'TE', 'NID', 'WBIn', 'WBOut'),
                          verbose = TRUE){
  
  # Will either be inferred - or - reference... use "unique" ID (comid facepalm).
  
  valid_types = c('HUC12', 'Gages', 'TE', 'NID', 'WBIn', 'WBOut', "Conf", "Term", "Elev", "Travel", "Con")
  
  if(!all(type %in% valid_types)){
    bad_ids = type[!which(type %in% valid_types)]
    stop(bad_ids, " are not valid POI types. Only ", paste(valid_types, collapse = ", "), " are valid")
  }
  
  if(is.null(type)){
    type = valid_types
  }
  
  nexus_locations  = mapped_POIs %>%
    st_drop_geometry()  %>%
    select(ID, identifier, paste0("Type_", type)) %>%
    mutate_at(vars(matches("Type_")), as.character) %>%
    mutate(poi_id = as.character(identifier),
           identifier = NULL) %>%
    group_by(ID, poi_id) %>%
    ungroup() %>%
    pivot_longer(-c(poi_id, ID)) %>%
    filter(!is.na(value)) %>%
    mutate(type = gsub("Type_", "", name),
           name = NULL)
  
  
  nexus_locations
}

# deprecated function from nhdplusTools -- should remove from package but here for compatibility
stage_national_data <- function(include = c("attribute",
                                            "flowline",
                                            "catchment"),
                                output_path = NULL,
                                nhdplus_data = NULL,
                                simplified = TRUE) {
  
  if (is.null(output_path)) {
    output_path <- dirname(nhdplus_data)
    warning(paste("No output path provided, using:", output_path))
  }
  
  if (is.null(nhdplus_data)) {
    stop("must provide nhdplus_data")
  }
  
  allow_include <- c("attribute", "flowline", "catchment")
  
  if (!all(include %in% allow_include)) {
    stop(paste0("Got invalid include entries. Expect one or more of: ",
                paste(allow_include, collapse = ", "), "."))
  }
  
  outlist <- list()
  
  if (any(c("attribute", "flowline") %in% include)) {
    
    out_path_attributes <- file.path(output_path,
                                     "nhdplus_flowline_attributes.rds")
    out_path_flines <- file.path(output_path, "nhdplus_flowline.rds")
    
    if (!(file.exists(out_path_flines) | file.exists(out_path_attributes))) {
      fline <- sf::st_zm(sf::read_sf(nhdplus_data,
                                     nhdplusTools:::get_flowline_layer_name(nhdplus_data)))
    }
    
    if ("attribute" %in% include) {
      if (file.exists(out_path_attributes)) {
        warning("attributes file exists")
      } else {
        saveRDS(sf::st_set_geometry(fline, NULL), out_path_attributes)
      }
      outlist["attributes"] <- out_path_attributes
    }
    
    if ("flowline" %in% include) {
      if (file.exists(out_path_flines)) {
        warning("flowline file exists")
      } else {
        saveRDS(fline, out_path_flines)
      }
      outlist["flowline"] <- out_path_flines
    }
  }
  
  if (exists("fline")) rm(fline)
  
  if ("catchment" %in% include) {
    out_path_catchments <- file.path(output_path, "nhdplus_catchment.rds")
    if (file.exists(out_path_catchments)) {
      warning("catchment already exists.")
    } else {
      
      layer_name <- nhdplusTools:::get_catchment_layer_name(simplified, nhdplus_data)
      
      saveRDS(sf::st_zm(sf::read_sf(nhdplus_data, layer_name)),
              out_path_catchments)
    }
    outlist["catchment"] <- out_path_catchments
  }
  
  return(outlist)
}