utils.R

needs_layer <- function(db, layer) {
  if(file.exists(db)) {
    layers <- st_layers(db)
    if(layer %in% layers$name)
      return(FALSE)
  }
  TRUE
}

layer_exists <- function(db, layer) {
  if(file.exists(db)) {
    layers <- st_layers(db)
    if(layer %in% layers$name)
      return(TRUE)
  }
  FALSE
}

# Subset NHD RDS by hydrologic region
VPU_Subset <- function(nhdPath, VPU){
  # Read in COMIDs of outlets
  RegOutlets <- jsonlite::read_json(file.path("cache", "RegOutlets.json"))
  print (paste0("subsetting for VPU ", VPU))
  
  if (!VPU %in% names(RegOutlets)){
    
    # 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 <- readRDS(nhdPath) %>% 
      filter(grepl(paste0("^", VPU, ".*"), .data$VPUID)) %>% filter(FTYPE != "Coastline")
    
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE)
    
    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]
    #print (outlet)
    
    nhd <- readRDS(nhdPath)
    
    #nhd_US <- get_UT(nhd, outlet)
    
    nhd <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, outlet)))
     
    keep <- prepare_nhdplus(nhd, 0, 0, 0, FALSE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID, VPUID == VPU)
    
  } else if (VPU == "10L"){

    # Subset by flowlines in Region 10
    nhd <- readRDS(nhdPath) %>% 
      filter(grepl(paste0("^",substr(VPU,1,2),".*"), .data$VPUID)) # take out this VPUID
    
    # R10U upstream
    nhd_US_10U <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, RegOutlets["10U"])))
    
    # R10L upstream
    nhd_US_10L <- nhd %>% filter(COMID %in% unlist(get_UT(nhd, RegOutlets["10L"])))
    
    nhd_10L <- nhd_US_10L %>% filter(!COMID %in% nhd_US_10U$COMID)
    
    keep <- prepare_nhdplus(nhd_10L, 0, 0, 0, FALSE)
    
    nhd <- filter(nhd, COMID %in% keep$COMID)
    
  } else {
    
    outlets <- unlist(RegOutlets[VPU])
    
    nhd <- readRDS(nhdPath) 
    
    #outlets <- c(14320629, 941140164, 15334434)
    
    for (out in outlets){
      #print (out)
      nhd_US <- get_UT(nhd, out)
      keep <- prepare_nhdplus(nhd %>% filter(COMID %in% unlist(nhd_US)), 0, 0, 0, FALSE)
      #print (dim(keep))
      ifelse( exists("final") , final <- rbind(final, keep), final <- keep)
    }
    
    nhd <- nhd %>% filter(COMID %in% final$COMID)
  }
  return (nhd)
}

return_NonDendirtic <- function(VPU, nhd, nhdPath){
  # Returns non-dendritic flowlines for assignment of coastal/non-dendritic catchments/flowlines
  nhdAll <- readRDS(nhdPath) %>%
    filter(grepl(paste0("^",VPU,".*"), .data$VPUID))
    nhd_ND <- nhdAll %>% filter(!COMID %in% nhd$COMID) %>% filter(TerminalFl == 1)
    
  return(nhd_ND)
  }

Merge_hydReg <- function(feat){
  out_gpkg <-"cache/GF_CONUS.gpkg"
  if(needs_layer(out_gpkg, feat)) {
    feat_01 <- read_sf("cache/GF_01.gpkg", paste0(feat,"_01")) 
    feat_02 <- read_sf("cache/GF_02.gpkg", paste0(feat,"_02"))
    feat_03 <- read_sf("cache/GF_03.gpkg", paste0(feat,"_03"))
    feat_04 <- read_sf("cache/GF_04.gpkg", paste0(feat,"_04"))
    feat_05 <- read_sf("cache/GF_05.gpkg", paste0(feat,"_05"))
    feat_06 <- read_sf("cache/GF_06.gpkg", paste0(feat,"_06"))
    feat_07 <- read_sf("cache/GF_07.gpkg", paste0(feat,"_07"))
    feat_08 <- read_sf("cache/GF_08.gpkg", paste0(feat,"_08"))
    feat_09 <- read_sf("cache/GF_09.gpkg", paste0(feat,"_09"))
    feat_10U <- read_sf("cache/GF_10U.gpkg", paste0(feat,"_10U"))
    feat_10L <- read_sf("cache/GF_10L.gpkg", paste0(feat,"_10L"))
    feat_11 <- read_sf("cache/GF_11.gpkg", paste0(feat,"_11"))
    feat_12 <- read_sf("cache/GF_12.gpkg", paste0(feat,"_12"))
    feat_13 <- read_sf("cache/GF_13.gpkg", paste0(feat,"_13"))
    feat_14 <- read_sf("cache/GF_14.gpkg", paste0(feat,"_14"))
    feat_15 <- read_sf("cache/GF_15.gpkg", paste0(feat,"_15"))
    feat_16 <- read_sf("cache/GF_16.gpkg", paste0(feat,"_16"))
    feat_17 <- read_sf("cache/GF_17.gpkg", paste0(feat,"_17"))
    feat_18 <- read_sf("cache/GF_18.gpkg", paste0(feat,"_18"))

    allfeat <- list(feat_01, feat_02, feat_03, feat_04, feat_05, feat_06, feat_07, feat_08, feat_09,
                    feat_10L, feat_10U, feat_11, feat_12, feat_13, feat_14, feat_15, feat_16, feat_17,
                    feat_18)
                  
  
    all_sf<-do.call(rbind, allfeat)
    write_sf(all_sf, out_gpkg, feat)
  }
}

siteAtt <- function(out_gpkg, hydReg){
  
  data_paths <- jsonlite::read_json(file.path("cache", "data_paths.json"))
  
  # Layers
  out_gpkg <- paste0("cache/GF_", hydReg, ".gpkg")
  gagesIII_pois <- read_sf(out_gpkg, paste0("GagesIII_", hydReg)) # GAGESIII POIs
  TE_pois <- read_sf(out_gpkg,paste0("TE_", hydReg)) # Thermoelectric POIs
  NID_pois <- read_sf(out_gpkg,paste0("NID_", hydReg)) # NID POIs
  nsegment_raw <- read_sf(out_gpkg,paste0("nsegment_raw_", hydReg)) # Minimally-sufficient network attributed with POI_ID
  nsegment <- read_sf(out_gpkg,paste0("Nsegment_", hydReg)) # Minimally-sufficient network dissolved by POI_ID
  WBs <- read_sf(out_gpkg,paste0("WB_", hydReg))
  WBout_POIs <- read_sf(out_gpkg,paste0("WBout_", hydReg))
  WBin_POIs <- read_sf(out_gpkg,paste0("WBin_", hydReg))
  
  #*************************************************************************
  # Attribute gagesIII POIs with adjusted measure and drainage area information
  gagesIII <- read_sf(data_paths$gagesiii_points_path, "GagesIII_gages") %>% 
    filter(grepl(paste0("^", hydReg, ".*"), .data$NHDPlusReg)) 
  
  # Load GAGESIII POIs from gpkg at cache dir
  gagesIIIPOIs <- gagesIII_pois %>% st_drop_geometry(.) %>% select(Type_GagesIII, COMID) %>%
    inner_join(gagesIII, by = "COMID") %>% select(Type_GagesIII, COMID, COMID_meas) 
  
  # Derive segment measure attributes
  nsegment_gages <- st_drop_geometry(nsegment_raw) %>% left_join(gagesIIIPOIs, by = "COMID") %>%
    select(Hydroseq, Type_GagesIII, COMID_meas, LENGTHKM, POI_ID, TotDASqKM) %>%
    mutate(COMID_meas = 1 - replace_na(COMID_meas/100,0), Seg_meas = COMID_meas * LENGTHKM) %>% group_by(POI_ID) %>%
    summarize(HS = min(Hydroseq), Length = sum(LENGTHKM), Percent = sum(Seg_meas)/sum(LENGTHKM), 
              TotDASqKM = max(TotDASqKM), Type_GagesIII =  first(na.omit(Type_GagesIII))) %>% filter(!is.na(Type_GagesIII))
  
  # Calcualte estimated proporational drainage area
  gagesIIIPOI_DA <- gagesIIIPOIs %>% left_join(nsegment_raw, by = "COMID") %>%
    select(Type_GagesIII, COMID, COMID_meas, AreaSqKM, TotDASqKM) %>%
    mutate(Inc_DA = (1-(COMID_meas / 100)) * AreaSqKM) %>%
    mutate(Total_DA = TotDASqKM - (AreaSqKM - Inc_DA)) 

  # PUt everything together
  gagesIII_stats <- gagesIIIPOI_DA %>% inner_join(nsegment_gages, by = "Type_GagesIII") %>%
    select(Type_GagesIII, COMID, Length, Percent, TotDASqKM.y) %>%
    rename(Seg_Length = Length, Seg_Measure = Percent, Adj_NHD_DA = TotDASqKM.y) %>%
    mutate(VPU = hydReg)
  
  # Compile data into GagesIII folder
  if (file.exists("data/GAGESIII_gages/gagesIII_MDA.rds")){
    # open file and bind to it
    curRDS_gages <- readRDS("data/GAGESIII_gages/gagesIII_MDA.rds") #%>% filter(VPU != hydReg) %>%
    bind_rows(gagesIII_stats)
    saveRDS(curRDS_gages, "data/GAGESIII_gages/gagesIII_MDA.rds")
  } else {
    saveRDS(gagesIII_stats, "data/GAGESIII_gages/gagesIII_MDA.rds")
  }
  #************************************************************************

  #************************************************************************
  # Read in TE shapefile
  TE_fac <- st_read(file.path(data_paths$TE_points_path, "/2015_TE_Model_Estimates_lat.long_COMIDs.shp")) %>%
    filter(COMID %in% TE_pois$COMID)

  # Cast TE as points and project to USGS albers
  TE_fac_alb <- TE_fac %>% st_cast('POINT') %>% st_transform(., 6703)

  # Subset segs based on TE and project to USGS albs
  TE_seg_pnts <- nsegment %>% filter(POI_ID %in% TE_fac_alb$COMID) %>% st_transform(., 6703) %>% group_by(POI_ID) %>%
    st_line_merge(.) %>% st_cast("POINT")

  #for (POI_ID in TE_fac_alb$COMID){
  TE_fac_alb$TE_measure <- unlist(lapply(1:nrow(TE_fac_alb), function(x){
    COM <- TE_fac_alb$COMID[x]
    fac <- TE_fac_alb %>% filter(COMID == COM)
    seg <- TE_seg_pnts %>% filter(POI_ID == COM)

    # Get index of nearest vertex on segment and rebuild two line IDs, then get length of each
    nearVert <- st_nearest_feature(fac, seg)
    # Pull out the line upstream and downstream of snapping point
    Upline <- seg[1:nearVert,] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)
    Downline <- seg[nearVert:nrow(seg),] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)

    # Return the measure
    return (Upline / (Upline + Downline))
  }))

  # Derive proportional total (probably need proporational incremental ratio too)
  TE_data <- st_drop_geometry(TE_fac_alb) %>%
    left_join(st_drop_geometry(nsegment_raw) %>% select(COMID, AreaSqKM, TotDASqKM), by = "COMID") %>%
    mutate(Inc_DA = TE_measure * AreaSqKM, Total_DA = TotDASqKM - (AreaSqKM - Inc_DA)) %>%
    select(EIA_PLANT_, LATITUDE, LONGITUDE, COMID, TE_measure, Inc_DA, Total_DA)

  TE_data_upstream <- TE_data %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("COMID" = "To_POI_ID")) %>%
    group_by(COMID) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " "))

  TE_data_final <- TE_data_upstream %>%
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("COMID" = "POI_ID")) %>%
    group_by(COMID) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " ")) %>%
    rename(Downstream_POI = To_POI_ID) %>% select(-POI_ID) %>% distinct() %>% mutate(VPU = hydReg)

  if (file.exists("data/TE_points/TE_adj.rds")){
    # open file and bind to it
    curRDS_TE <- readRDS("data/TE_points/TE_adj.rds") %>% filter(VPU != hydReg) %>%
      bind_rows(TE_data_final)
    saveRDS(curRDS_TE, "data/TE_points/TE_adj.rds")
  } else {
    saveRDS(TE_data_final, "data/TE_points/TE_adj.rds")
  }
  #************************************************************************
  
  #************************************************************************
  # Read in NID shapefile
  NID_fac <- read.csv(file.path(data_paths$NID_points_path,"/NID_attributes_20170612.txt"), stringsAsFactors = F) %>%
    filter(FlowLcomid %in% NID_pois$COMID)
  
  NID_snap_alb <- readRDS(file.path(data_paths$NID_points_path,"/NAWQA_NID_snap.rds")) %>%
    filter(Source_FeatureID %in% NID_fac$Source_FeatureID) %>% st_transform(., 6703) %>%
    inner_join(NID_fac, by = "Source_FeatureID")
  
  # Subset segs based on TE and project to USGS albs
  NID_seg_pnts <- nsegment %>% filter(POI_ID %in% unique(NID_fac$FlowLcomid)) %>% st_transform(., 6703) %>% group_by(POI_ID) %>%
    st_line_merge(.) %>% st_cast("POINT")
  
  #options(warn = -1)
  NID_snap_alb$NID_measure <- unlist(lapply(1:nrow(NID_snap_alb), function(x){
    fac <- NID_snap_alb[x,]
    COM <- fac$FlowLcomid
    seg <- NID_seg_pnts %>% filter(POI_ID == COM)
    
    # Get index of nearest vertex on segment and rebuild two line IDs, then get length of each
    nearVert <- st_nearest_feature(fac, seg)
    # Pull out the line upstream and downstream of snapping point
    Upline <- seg[1:nearVert,] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)
    Downline <- seg[nearVert:nrow(seg),] %>% summarise(geom = st_combine(geom)) %>% st_cast("LINESTRING") %>% st_length(.)
    
    # Return the measure
    return (Upline / (Upline + Downline))
  }))
  #options(warn = oldw)
  
  # Derive proportional total (probably need proporational incremental ratio too)
  NID_data <- st_drop_geometry(NID_snap_alb) %>% 
    left_join(st_drop_geometry(nsegment_raw) %>% select(COMID, AreaSqKM, TotDASqKM), by = c("FlowLcomid" = "COMID")) %>%
    mutate(Inc_DA = NID_measure * AreaSqKM, Total_DA = TotDASqKM - (AreaSqKM - Inc_DA)) %>%
    select(Source_FeatureID, VPU.x, Outlet, FlowLcomid, WBAreaSQKM, NIDSASQKM, NLCDSASQKM, MAX_SA, MultiMain, Catchcomid, FalconeDate, 
           YEAR_COMPLETED, MAX_DISCHARGE, NID_STORAGE, NID_measure, Inc_DA, Total_DA)
  
  NID_data_upstream <- NID_data %>% 
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("FlowLcomid" = "To_POI_ID")) %>%
    group_by(FlowLcomid) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " ")) 
  
  NID_data_final <- NID_data_upstream %>% 
    left_join(st_drop_geometry(nsegment) %>% select(POI_ID, To_POI_ID), by = c("FlowLcomid" = "POI_ID")) %>%
    group_by(FlowLcomid) %>% mutate(Upstream_POIs = paste0(unique(POI_ID), collapse = " ")) %>%
    rename(Downstream_POI = To_POI_ID) %>% select(-POI_ID) %>% distinct() %>% mutate(VPU = hydReg)
  
  if (file.exists("data/NID_points/NID_adj.rds")){
    # open file and bind to it
    curRDS_NID <- readRDS("data/NID_points/NID_adj.rds") #%>% filter(VPU != hydReg) %>%
    bind_rows(NID_data_final)
    saveRDS(curRDS_NID, "data/NID_points/NID_adj.rds")
  } else {
    saveRDS(NID_data_final, "data/NID_points/NID_adj.rds")
  }
}