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_rds, new_atts = NULL, 
                           nhdpath = nhdplusTools::nhdplus_path(), 
                           par = 6) {
  
  if(!file.exists(out_rds)) {
  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))
  fline <- 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, StreamOrde, 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
  
  saveRDS(nhd, out_rds)
  
  ### ONLY RUN IF CHANGED ###
  # sbtools::authenticate_sb()
  # sbtools::item_replace_files("5dcd5f96e4b069579760aedb", "data/NHDPlusNationalData/nhdplus_flowline_update.rds")
  }
  return(out_rds)
}

#' 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)
  }
}

#' Merges geopackages together to create CONUs geopackage of features
#' @param VPU VPU from NHDPlusV2
#' @param full_cats  (sf data.frame) all catchments (sinks/flowline) within given VPU
#' @param divides (sf data.frame) divides layer for a given VPU
#' @param wbd  (data.frame) HUC12 X-walk table, see get_data for full citation
#' @param wbd_SF (sf data.frame) HUC12 layer from NHDPlusv2.1 national geodatabase
#' @return (sf data.frame) intersection of nhd catchments and divides layer with HUC12
#'
get_HUC12_Xwalk <- function(VPU, full_cats = FALSE, divides = FALSE, wbd = FALSE, wbd_SF = FALSE){

  # HUC02 includes some 
  if(VPU == "02"){
    grep_exp <-"^02|^04"
  } else if (VPU == "08") {
    grep_exp <- "^03|^08"
  } else {
    grep_exp <- paste0("^", substr(VPU, start = 1, stop = 2))
  }
  
  # crosswalk to HU12, filter by VPU
  nhd_to_HU12 <- read.csv(wbd, colClasses = c("character", "integer", "character")) %>% 
    filter(grepl(grep_exp, .data$HUC_12)) %>%
    mutate(HUC_8 = substr(HUC_12, 1, 8))
  
  # Read in the NHDPlus National GDB HUC12 layer
  HUC12_lyr <- readRDS(file.path(wbd_SF)) %>% 
    #filter(grepl(paste0("^", VPU,".*"), .data$HUC_12)) %>%
    filter(grepl(grep_exp, .data$HUC_12)) %>% #, 
    st_set_precision(10) %>%
    st_transform(st_crs(full_cats)) %>%
    st_make_valid() %>%
    select(HUC_12, ACRES, HU_10_TYPE, HU_12_DS, HU_12_MOD, HU_12_TYPE)

  reg_cats <- full_cats %>% 
    left_join(select(nhd_to_HU12, c(FEATUREID, HUC_12)), by = "FEATUREID")
  
  # Intersect the HUC12s and catchments
  cats_HUC12 <- sf::st_intersection(reg_cats, HUC12_lyr) %>%
    # Convert area to squaker kilometers 
    mutate(intArea = as.numeric(st_area(.)) * 0.000001,
           AreaHUC12 = ACRES * .00404686)
  
  divides_HUC12 <- sf::st_intersection(divides, HUC12_lyr) %>%
    # Convert area to squaker kilometers 
    mutate(intArea = as.numeric(st_area(.)) * 0.000001,
           AreaHUC12 = ACRES * .00404686)
  
  
  return(list(cats_HUC12 = cats_HUC12, divides_HUC12 = divides_HUC12))
}

#' Fix flow direction
#' @description If flowlines aren't digitized in the expected direction, 
#' this will reorder the nodes so they are.
#' @param x The COMID of the flowline to check
#' @param r The entire network to check from. Requires a "toCOMID" field.
#' @return a geometry for the feature that has been reversed if needed.
fix_flowdir <- function(comid, network) {
try({
  library(sf)
  
  f <- network[network$COMID == comid, ]
  
  if(is.na(f$toCOMID)) {
    
    check_line <- network[network$toCOMID == f$COMID, ][1, ]
    
    check_position <- "start"
    
  } else {
    
    check_line <- network[network$COMID == f$toCOMID, ][1, ]
    
    check_position <- "end"
    
  }
  
  suppressMessages(check_end <- sf::st_join(nhdplusTools::get_node(f, position = check_position), 
                                            dplyr::select(check_line, check_COMID = COMID)))
  
  reverse <- is.na(check_end$check_COMID)
  
  if(reverse) {
    sf::st_geometry(f)[reverse] <- sf::st_reverse(sf::st_geometry(f)[reverse])
  }
  
  return(sf::st_geometry(f))
})
}


#' 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, vpu){
  
  fl <- read_sf(nhd_gdb, "NHDFlowline_Network") %>%
    align_nhdplus_names()
  
  # read the CatchmentSP
  cats <- 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)
}

### DEPRECATED moved to nhdplusTools as subset_vpu with different methods###
#' Subsets VPU for terminal outlets
#' @param nhd_rds path to staged nhd flowline rds file
#' @param VPU VPU from NHDPlusV2
#' @return sf data.frame containing subset VPU
#' 
VPU_Subset <- function(nhd_rds, VPU){
  
  # Read in COMIDs of outlets
  regoutlets <- list(`05` = 1844789L, 
                     `06` = 1861888L, 
                     `07` = 5093446L, 
                     `10U` = 17244390L, 
                     `10L` = 6018266L, 
                     `14` = c(20734041L, 18267749L), 
                     `11` = c(14320629L, 941140164L, 15334434L))
  
  nhd <- readRDS(nhd_rds)
  
  if (!VPU %in% names(regoutlets[!names(regoutlets) %in% c("10L", "10U")])){
    
    # For regions that terminate to the NHDPlus domain
    # 1,2,3,4,8,9,12,15,17,18
    # Read in NHDPlusV2 flowline simple features and filter by vector processing unit (VPU)
    nhd <- nhd %>% 
      filter(grepl(paste0("^", VPU, ".*"), .data$VPUID)) %>% 
      filter(FTYPE != "Coastline")
    
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE, skip_toCOMID = TRUE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID)
    
  } else if (VPU %in% c("05", "06", "07", "10U", "14")){
    
    # Subset VPUs that are single outlet
    outlet <- regoutlets[VPU]
    
    nhd <- filter(nhd, COMID %in% unlist(get_UT(nhd, outlet)))
    
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE, skip_toCOMID = TRUE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID, VPUID == VPU)
    
  } else {
    
    outlets <- unlist(regoutlets[VPU])
    
    for (out in outlets){
      
      nhd_US <- get_UT(nhd, out)
      
      keep <- prepare_nhdplus(filter(nhd, COMID %in% unlist(nhd_US)), 
                              0, 0, 0, FALSE, 
                              skip_toCOMID = TRUE)
      
      ifelse( exists("final") , 
              final <- rbind(final, keep), 
              final <- keep)
      
    }
    
    nhd <- filter(nhd, COMID %in% final$COMID)
    
  }
  
  nhd <- select(nhd, -VPUID, -FromNode, -ToNode, 
                -StartFlag, -StreamOrde, -StreamCalc)
  
  return(st_as_sf(nhd))
  
}

### DEPRECATED moved to nhdplusTools as is ####
# make sf1 compatible with sf2
st_compatibalize <- function(sf1, sf2) {
  sf1 <- st_transform(sf1, st_crs(sf2))
  
  g <- attr(sf1, "sf_column")
  gp <- attr(sf2, "sf_column")
  names(sf1)[names(sf1) == g] <- gp
  attr(sf1, "sf_column") <- gp
  sf1
}

### DEPRECATED moved to nhdplusTools as is ####
rename_geometry <- function(g, name){
  current = attr(g, "sf_column")
  
  names(g)[names(g)==current] = name
  
  st_geometry(g) <- name
  
  g
}