In NOAA-OWP/hyAggregate: Hydrologic and Hydrographic Network Aggregating Tools Based on HY_Features

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=7,
  fig.height=4,
  warning = FALSE,
  message = FALSE
)

library(hyAggregate)
library(dplyr)


network_file         = "/Volumes/Transcend/ngen/enhd_nhdplusatts.parquet"
reference_flowlines  = "/Volumes/Transcend/ngen/reference_CONUS.gpkg"
reference_catchments = "/Volumes/Transcend/ngen/reference_catchments.gpkg"
facfdr               = "/Volumes/Transcend/ngen/fdrfac_cog"
outfile              = "/Volumes/Transcend/ngen/tmp.gpkg"



source('/Users/mjohnson/github/hyRelease/R/update_functions.R')


outfile     = get_UT_reference(
                            network= arrow::read_parquet(network_file),
                            reference_flowlines =  reference_flowlines,
                            comid = 3766334,
                            outfile = outfile)

og  = sf::read_sf(outfile, "reference_flowpaths")


refactor_wrapper(flowpaths  = sf::read_sf(outfile, "reference_flowpaths"),
                 catchments = NULL,
                 facfdr     = facfdr,
                 outfile    = outfile)


ref  = sf::read_sf(outfile, "refactored_flowpaths")


agg = aggregate_by_thresholds(gpkg = outfile,
                              fl_name = "refactored_flowpaths",
                              min_length_km = 1,
                              outfile = outfile)


og$LENGTHKM |> range()
ref$LENGTHKM |> range()
agg$flowpaths$lengthkm |> range()

agg$flowpaths  |> 
  mapview::mapview()

nex = hyAggregate::get_nexus_locations(agg$flowpaths)

# mapview::mapview(agg) + nex
# 
# mean(base$cats$areqsqkm)
# mean(base$fps$LENGTHKM)
# 
# mean(rec$catchments$areasqkm)
# mean(rec$flowpaths$LENGTHKM)
# 
# mean(agg$catchments$areasqkm)
# mean(agg$flowpaths$lengthkm)
# 
# sum(st_length(base$fps)) /1e6
# sum(st_length(rec$flowpaths)) /1e6
# sum(st_length(agg$flowpaths)) /1e6
# 
# sum(base$cats$areqsqkm < 3)
# sum(rec$catchments$areasqkm < 3)
# sum(agg$catchments$areasqkm < 3)
# 
# sum(base$fps$LENGTHKM < 1)
# sum(rec$flowpaths$lengthkm < 1)
# sum(agg$flowpaths$lengthkm < 1)

cat = filter(agg$catchments, ID %in% c(182, 183,184,187,188)) %>% 
  mutate(ID = case_when(ID == 182 ~ 5,
                        ID == 183 ~ 3,
                        ID == 184 ~ 1,
                        ID == 187 ~ 2,
                        ID == 188 ~ 4)) %>% 
  mutate(ID = paste("cat-", ID))

fp  = filter(agg$flowpaths, ID %in% c(182, 183,184,187,188)) %>% 
    mutate(toID = case_when(toID == 182 ~ 2,
                            toID == 183 ~ 1,
                            toID == 29  ~ 3)) %>% 
  mutate(toID = paste0("toID: nex-", toID))

fp$toID

nex2 = nex %>% 
  filter(ID %in% c(182, 183)) %>% 
  mutate(ID = case_when(ID == 182 ~ 2,
                        ID == 183 ~ 1)) %>% 
  mutate(ID = paste0("nex-", ID))

write_sf(fp, "/Users/mjohnson/Desktop/example.gpkg", "fps")

write_sf(cat , "/Users/mjohnson/Desktop/example.gpkg", "cats")

write_sf(nex2 , "/Users/mjohnson/Desktop/example.gpkg", "nex")



agg$flowpaths$ID = paste0("ID: ", agg$flowpaths$ID)
agg$flowpaths$toID = paste0("toID: ", agg$flowpaths$toID)

library(ggplot2)
library(patchwork)

base_ln = ggplot(data = base$fps, aes(x = LENGTHKM)) +
  geom_histogram(bins= 100, fill="#69b3a2", color="#69b3a2", alpha=0.9)+
  geom_vline(xintercept = 1, lwd = 1) +
  labs(x = "Length (KM)", title = paste0("Reference Flowlines\n (", nrow(base$fps), ")"))

rfc_ln = ggplot(data = ref$fps, aes(x = LENGTHKM)) +
  geom_histogram(bins= 100, fill="#69b3a2", color="#69b3a2", alpha=0.9) +
  geom_vline(xintercept = 1, lwd = 1) +
  labs(x = "Length (KM)", title = paste0("Refactored Flowlines\n (",nrow(ref$fps), ")"))

agg_ln = ggplot(data = agg$flowpaths, aes(x = lengthkm)) +
  geom_histogram(bins= 100, fill="#69b3a2", color="#69b3a2", alpha=0.9)+
  geom_vline(xintercept = 1, lwd = 1) +
  labs(x = "Length (KM)", title = paste0("Aggregated Flowlines\n (", nrow(agg$flowpaths), ")"))

#####

base_ar = ggplot(data = base$cats, aes(x = areqsqkm)) +
  geom_histogram(bins= 100, fill="#404080", color="#404080", alpha=0.9)+
  geom_vline(xintercept = 3, lwd = 1) +
  geom_vline(xintercept = 10, lwd = 1) +
  geom_vline(xintercept = 20, lwd = 1) +
  xlim(0,30) +
  labs(x = "Area (KM2)", title = paste0("Refactored Catchments\n (", nrow(base$cats), ")"))

rfc_ar = ggplot(data = ref$cats, aes(x = areasqkm)) +
  geom_histogram(bins= 100, fill="#404080", color="#404080", alpha=0.9)+
  geom_vline(xintercept = 3, lwd = 1) +
  geom_vline(xintercept = 10, lwd = 1) +
  geom_vline(xintercept = 20, lwd = 1) +
  xlim(0,30) +
  labs(x = "Area (KM2)", title = paste0("Aggregated Catchments\n (", nrow(ref$cats), ")"))

agg_ar = ggplot(data = agg$catchments, aes(x = areasqkm)) +
  geom_histogram(bins= 100, fill="#404080", color="#404080", alpha=0.9)+
  geom_vline(xintercept = 3, lwd = 1) +
  geom_vline(xintercept = 10, lwd = 1) +
  geom_vline(xintercept = 20, lwd = 1) +
  xlim(0,30) +
  labs(x = "Area (KM2)", title = paste0("Aggregated Catchments\n (", nrow(agg$catchments), ")"))

base_ln + rfc_ln + agg_ln


base_ar + rfc_ar + agg_ar

NOAA-OWP/hyAggregate documentation built on Aug. 13, 2022, 1:36 p.m.