test_nwm_3dhp.qmd

---
title: "NWM on 3DHP testing"
author: "dan.auerbach@dfw.wa.gov"
date: "`r Sys.Date()`"
format:
  html:
    embed-resources: true
    theme: yeti 
    code-fold: true
    toc: true
    toc-location: left
    grid:
      sidebar-width: 180px
      body-width: 1100px
      margin-width: 20px
---

to-do
  - culverts
    - test length*slope (width?) as non-dollar relative 'cost' indicator
    - hydgeo scaling (on upstream length as DA proxy) for est width & depth - but can/should be able to get directly from DEM?
  - stream type
    - tabulations


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE, fig.width = 8, fig.height = 8)

library("tidyverse", quietly = T)
library("sf")
library("patchwork")
library("gt")
theme_set(theme_minimal()) 

dir_data_common <- "~/T/DFW-Team WDFW Watershed Synthesis - data_common"
epsg <- 2927 #WA state standard; NAD83(HARN)/ft

load("test_nwm_3dhp.RData")
```

# objects

```{r lines_3dhp, eval=FALSE}
# sf::st_layers(file.path(dir_data_common, "3dhp/17110008_EDH_230605.gdb"))

# #after initial read before select()
# as_tibble(lines) |> #glimpse()
#   #count(FClass, EClass) #Fclass all 1, Eclass 0,2,3
#   #count(Source) #all "Lidar YYYY", mostly 2017 and 2013
#   #count(Method) #all "Flow direction and accumulation combined with proprietary methods"
#   #count(UserCode) #all NA except 2 " "
#   #count(Comments) #mostly NA, sort of interesting but probably not standard, in USGS?
#   count(GNIS_Name) #100 names and a few versions of blanks and nulls

# #test if multilinestring is actually more than single line
# #not, anywhere in the full set actually
# any(map_int(st_geometry(lines_3dhp), length) > 1)

#in 2927, 74K, note XYZ-ness
lines_3dhp <- sf::st_read(file.path(dir_data_common, "3dhp/17110008_EDH_230605.gdb"), layer = "Lines") |> 
  select(-c(FClass, Method, UserCode, Comments)) |> 
  sf::st_cast("LINESTRING") 

# lines_3dhp |> 
#   #as_tibble() |> glimpse()
#   ggplot() + geom_sf(aes(color = Source))

# #not clear exactly where this is...
# cat_3dhp <- sf::st_read(file.path(dir_data_common, "3dhp/1711000802_catchments/1711000802_catchments.shp"))

```

could add custom function to `slope_xyz` to calc max gradient along reach rather than length-weighted default

```{r add_slope_est_convex_hull, eval=FALSE}
lines_3dhp <- lines_3dhp |> 
  mutate(
    across(SHAPE, list(
      zslope_est = ~slopes::slope_xyz(., directed = T),
      zmax = ~map_dbl(., ~slopes::z_max(.x))
      ), .names = "{.fn}")
  )

# lines_3dhp |> 
#   st_centroid() |> 
#   ggplot() + geom_sf(aes(color = zslope_est), size = 0.5, alpha = 0.5) + wacolors::scale_color_wa_c("forest_fire")

#slightly tighter than bbox
lines_3dhp_uch <- lines_3dhp |> st_union() |> st_convex_hull()

```

cut nhdp to 3dhp & associate comid

```{r sf_nhdp_wa, eval=FALSE}
# #med res NHDplus, 230K polys including non-WA, no attrib but FEATUREID==COMID for StreamCat join
# sf_nhdp <- sf::read_sf(file.path(dir_data_common, "NHDPlus17/NHDPlusCatchment/Catchment.shp")) |>
#   sf::st_transform(crs = sf::st_crs(epsg)) |>
#   select(COMID = FEATUREID, areasqkm = AreaSqKM)
# 
# #use StreamCat for WA to nonspatially subset catchment polys
# sf_nhdp_wa <- inner_join(
#   sf_nhdp |> mutate(comid = as.character(COMID))
#   ,
#   StreamCatTools::sc_nlcd(state = "WA", year = "2019", aoi = 'riparian_watershed') |>
#     rename_with(tolower) |> as_tibble() |>
#     mutate(comid = as.character(comid)) |> 
#     select(comid, state)
#   ,
#   by = "comid")
# 
# saveRDS(sf_nhdp_wa, "~/T/DFW-Team WDFW Watershed Synthesis - flow_trees_heat/sf_nhdp_wa.rds")

sf_nhdp_wa <- readRDS("~/T/DFW-Team WDFW Watershed Synthesis - flow_trees_heat/sf_nhdp_wa.rds")
```

```{r sf_nhdp_wa_3dhp, eval=FALSE}
#sf_nhdp_wa_3dhp <- sf_nhdp_wa[st_as_sfc(st_bbox(lines_3dhp)),] #1464 of 57313
sf_nhdp_wa_3dhp <- sf_nhdp_wa[lines_3dhp_uch,] #959

#not bad, maybe ~1min?
lines_3dhp <- st_join(
  lines_3dhp,
  sf_nhdp_wa_3dhp |> select(-state),
  join = st_intersects, #st_crosses?
  largest = T
)

```

cut FPDSI culvert subset down to 3dhp

```{r sf_fpb, eval=FALSE}
#also: tbl(con_fpdsi, "SiteFishPassageViewSpatial") for current but need to revisit permissions for full schema

sf_fpb <- sf::read_sf(
  file.path(dir_data_common,"WdfwFishPassage.gdb"), 
  layer = "WDFW_FishPassageSite") |> 
  sf::st_transform(epsg) |> 
  rename_with(.cols = starts_with("FishPassage"), .fn = ~str_replace(., "FishPassage","FP")) |> 
  filter(
    FPFeatureTypeCode==1 #culverts
    # #!is.na(BarrierCorrectionTypeCode)
    # !is.na(BarrierCorrectionYearsText)
  ) |> 
  mutate(
    fish_use = case_when(
      FishUseCode == 10 ~ "yes",
      FishUseCode == 20 ~ "no",
      FishUseCode == 99 ~ "unk"
    )
  ) |> 
  left_join(
    tibble(
      OwnerTypeCode = c(1:7,9,10,12), 
      owner = c("city", "county", "federal", "private", "state", "tribal", "other", "drainage_dist", "diking_dist", "unk")
    ),
    by = "OwnerTypeCode"
  )

```

```{r sf_fpb_3dhp, eval=FALSE}
#sf_fpb_3dhp <- sf_fpb[st_as_sfc(st_bbox(lines_3dhp)),] #this is 2221
sf_fpb_3dhp <- sf_fpb[lines_3dhp_uch,] #1374

#as_tibble(sf_fpb_3dhp) |> count(fish_use)
#sf_fpb_3dhp |> ggplot() + geom_sf(aes(shape = fish_use, color = fish_use)) + geom_sf(data = lines_3dhp_uch, fill = NA, color = "hotpink")

```

add a basic barrier attribute count to NHDcats

```{r sf_nhdp_wa_3dhp_fpb, eval=FALSE}
sf_nhdp_wa_3dhp_fpb <- sf_nhdp_wa_3dhp |> 
  select(-state) |> 
  left_join(
    st_join(sf_fpb_3dhp, sf_nhdp_wa_3dhp) |> #pt in poly, not expanding
      as_tibble() |> 
      count(comid, fish_use) |> 
      pivot_wider(names_from = fish_use, values_from = n)
    ,
    by = "comid"
  )

# {ggplot() + geom_sf(data = sf_nhdp_wa_3dhp_fpb, aes(fill = yes), color = NA) + 
#     geom_sf(data = lines_3dhp_uch, fill = NA, color = "hotpink") +
#     wacolors::scale_fill_wa_c("forest_fire", na.value = "grey90")
#   }+{ggplot() + geom_sf(data = sf_nhdp_wa_3dhp_fpb, aes(fill = no), color = NA) + 
#       geom_sf(data = lines_3dhp_uch, fill = NA, color = "hotpink") +
#       wacolors::scale_fill_wa_c("ferries", na.value = "grey90")} +
#   plot_layout(ncol = 1) + 
#   plot_annotation(subtitle = "Fish passage barriers per 1:100NHDplus catchment by 'fish use' \n pink hull shows 3DHP extent")

```

trim NHDplusHighRes down to 3DHP general extent (convex hull of unioned lines) 

```{r sf_nhdphr_flw_3dhp, eval=FALSE}
sf_nhdphr_flw_3dhp <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "NHDFlowline") |> 
  select(Permanent_Identifier, GNIS_Name:ReachCode, FCode, Shape_Length, NHDPlusID) |> 
  st_cast("LINESTRING") |> 
  left_join(
    read_csv(file.path(dir_data_common, "nhdphr/brinkerhoffetal2024_results_archive/results_1711.csv"))
    ,
    by = "NHDPlusID") |> 
  drop_na(percQEph_reach) |> 
  st_zm() |> 
  st_transform(st_crs(lines_3dhp))

sf_nhdphr_flw_3dhp <- sf_nhdphr_flw_3dhp[lines_3dhp_uch,] #23785

```


cut 3dhp to focal comid(s) & build test network

```{r lines_3dhp_to_lsfn, eval=FALSE}
# #where is the greatest density of 3dhp lines in NHDplus MR catchments?
# sf_nhdp_wa_3dhp |> 
#   left_join(
#     lines_3dhp |> as_tibble() |> count(comid, name = "n_3dhp") |> arrange(desc(n_3dhp))    
#     , by = "comid"
#   ) |> 
#   drop_na(n_3dhp) |> 
#   ggplot() + geom_sf(aes(color = n_3dhp, fill = n_3dhp)) + 
#   wacolors::scale_color_wa_c("sea", reverse = T, aesthetics = c("color","fill")) + 
#   labs(title = "Pilot 3DHP lines per med res NHDplus catchments")
# 
# #for now just take the comid with the greatest number of associated 3dhp line units: "24275471"
# lines_3dhp |> as_tibble() |> count(comid, name = "n_3dhp") |> arrange(desc(n_3dhp)) |> slice(1:10)


# #spatial bounds initially generate a 'rooted forest' with GIS artifact strays
# ggplot() +
#   geom_sf(data = st_as_sf(lsfn, "edges"), col = "lightblue", linewidth = 0.3) +
#   geom_sf(
#     data = lsfn |> 
#       tidygraph::activate(nodes) |> 
#       filter(tidygraph::node_is_root()) |> 
#       st_as_sf("nodes")
#     ,
#     color = "orange"
#   )


#create sfnetwork object with nodes and edges that have from/to attributes added...
#identical feature count, follows from length(geometry()) test
#the convert to_largest_component drops GIS artifact strays from initial 'rooted forest'
#then adds calc'd topological and length weighted distance for nodes which are then associated to 'to' of edges
#map inspection that the ground distance versions are 'from downstream node to downstream [root] node'
lsfn <- lines_3dhp |> 
  filter(comid == "24275471") |> 
  select(uid = UniqueID, shape_length = SHAPE_Length, comid, desc = Desc) |> 
  sfnetworks::as_sfnetwork() |> 
  tidygraph::activate(nodes) |> 
  tidygraph::convert(
    tidygraph::to_largest_component,
    .clean = T #this drops the .tidygraph_node_index/edge_index that no longer apply & are potentially confusing
  ) |> 
  mutate(
    idn = seq_along(SHAPE),
    #raindrop perspective paths from nodes to root
    d_root = tidygraph::node_distance_to(tidygraph::node_is_root()),
    d_root_ft = tidygraph::node_distance_to(tidygraph::node_is_root(), weights = shape_length),
    d_root_km = 0.0003 * d_root_ft,
    #looking 'up' perspective to perimeter from per-node 'me' 
    in_n = tidygraph::local_size(order = nrow(tidygraph::.N()), mode = "in")
    ) 


#associate the immediately upstream streamlength from edges to each 'to' node
#summed over both upstream tribs
#also add 'below_culvert' logical:
#   every culvert edge has a single 'to' node,
#   but each of those 'to' nodes has potentially 2 incoming edges
#   one of which may be a non-culvert
#   and there is a potential case of 2 converging culverts
#   so need logic to define per-node state in 
#   {no-incoming, not, culvert, not+not, culvert+not, culvert+culvert}
lsfn <- lsfn |> 
  tidygraph::activate(nodes) |> 
  left_join(
    select(as_tibble(as_tibble(lsfn, active = "edges")), uid, idn = to, shape_length, desc) |> 
      summarise(
        below_culvert = if_else(any(str_detect(desc, "ulvert")), T, F),
        in_ft_first_order = sum(shape_length),
        .by = "idn"),
    by = c("idn")
    )

#ensure whole thing is directed acyclic to root
tidygraph::with_graph(lsfn, tidygraph::graph_is_dag()) #true, good

#now per-node, define upstream neighborhood and sum over first-order incoming lengths of neighborhood members
#sweet: https://tidygraph.data-imaginist.com/reference/map_local.html
#this may not be the most efficient approach but...
#no appreciable speed diff between order = 3,10,max
lsfn_order <- tidygraph::with_graph(lsfn, tidygraph::graph_order())

lsfn <- lsfn |> 
  tidygraph::activate(nodes) |> 
  mutate(
    in_ft = tidygraph::map_local_dbl(
      order = lsfn_order, mode = "in", 
      .f = function(neighborhood, ...){
        sum(as_tibble(neighborhood, active = 'nodes')$in_ft_first_order, na.rm = T)
        }
      ),
    in_km = 0.0003 * in_ft
  )


#and then add node attributes back to edge by 'to'
lsfn <- lsfn |> 
  tidygraph::activate(edges) |> 
  mutate(culvert = str_detect(desc, "ulvert")) |> 
  left_join(
    select(as_tibble(as_tibble(lsfn, active = "nodes")), to = idn, d_root, ends_with("_km")) 
    ,
    by = c("to")
    ) |> 
  tidygraph::activate(nodes)
  


## unused but keeping handy for now
## temporarily converting edges to nodes 
# lsfn |> 
#   tidygraph::morph(
#     tidygraph::to_linegraph
#   ) |> 
#   mutate(
#     d_root = tidygraph::node_distance_to(tidygraph::node_is_root())
#   ) |> 
#   tidygraph::unmorph() 

```

```{r func_prep_sfn}
#little wrapper to pull out 3dhp for a given catchment
#then make into attributed sfnetwork
#follow above flow in chunk lines_3dhp_to_lsfn
prep_sfn <- function(cat){
  l <- lines_3dhp |> 
    filter(comid == cat) |> 
    select(uid = UniqueID, shape_length = SHAPE_Length, comid, desc = Desc, zslope_est, zmax) |> 
    sfnetworks::as_sfnetwork() |>
    tidygraph::activate(nodes) |>
    tidygraph::convert(
      tidygraph::to_largest_component,
      .clean = T #this drops the .tidygraph_node_index/edge_index that no longer apply & are potentially confusing
    ) |>
    mutate(
      idn = seq_along(SHAPE),
      #raindrop perspective paths from nodes to root
      d_root = tidygraph::node_distance_to(tidygraph::node_is_root()),
      d_root_ft = tidygraph::node_distance_to(tidygraph::node_is_root(), weights = shape_length),
      d_root_km = 0.0003 * d_root_ft,
      #looking 'up' perspective to perimeter from per-node 'me'
      in_n = tidygraph::local_size(order = nrow(tidygraph::.N()), mode = "in")
    )

  #ensure whole thing is directed acyclic to root
  if(!tidygraph::with_graph(l, tidygraph::graph_is_dag())) {stop();print("Not DAG")}

  #associate the immediately upstream streamlength from edges to each 'to' node
  #summed over both upstream tribs; also add 'below_culvert' logical:
  # defining per-node state in {no-incoming, not, culvert, not+not, culvert+not, culvert+culvert}
  l <- l |>
    tidygraph::activate(nodes) |>
    left_join(
      select(as_tibble(as_tibble(l, active = "edges")), uid, idn = to, shape_length, desc, zslope_est, zmax) |>
        summarise(
          below_culvert = if_else(any(str_detect(desc, "ulvert")), T, F),
          in_slope_mean = mean(zslope_est),
          in_ft_first_order = sum(shape_length),
          .by = "idn"),
      by = c("idn")
    )

  #now per-node, define upstream neighborhood and sum over first-order incoming lengths of neighborhood members
  l_order <- tidygraph::with_graph(l, tidygraph::graph_order())
  
  l <- l |>
    tidygraph::activate(nodes) |>
    mutate(
      in_ft = tidygraph::map_local_dbl(
        order = l_order, mode = "in",
        .f = function(neighborhood, ...){
          sum(as_tibble(neighborhood, active = 'nodes')$in_ft_first_order, na.rm = T)
        }
      ),
      in_km = 0.0003 * in_ft
    )

  #and then add node attributes back to edge by 'to'
  l <- l |>
    tidygraph::activate(edges) |>
    mutate(culvert = str_detect(desc, "ulvert")) |>
    left_join(
      select(as_tibble(as_tibble(l, active = "nodes")), to = idn, d_root, ends_with("_km"))
      ,
      by = c("to")
    ) |>
    tidygraph::activate(nodes)
  
  return(l)
}


lsfn <- prep_sfn("24275471")
```


# perfuse fixed root value

## alt1: topology + channel distance

playing with reverse "contributing" + "travel time" to get a cfs-unit value 'up'

the "in_km" built from topology and edge geometry length is proxy, not true 'drainage area', especially over very large domains or other (lithology, land cover) factors that cause adjoining subbasins to have very different drainage density/dissection.

Need to progress iteratively from network perimeter down, at a multi-COMID scale?
This first pass going to create big flow 'jumps' where mainstem catchments join?


```{r dumb_instant_travel_perfuse}

z <- as_tibble(lsfn, active = 'nodes')

#z |> sf::st_zm() |> select(-in_node) |> mapview::mapview(zcol = "in_ft")

#hot damn
ggplot() +
  geom_sf(data = as_tibble(lsfn, active = "edges"), color = "lightblue", linewidth = 0.5) +
  geom_sf(data = z, aes(color = log10(in_ft)), alpha = 0.8) +
  wacolors::scale_color_wa_b()

#try for super dumb demo of 'spreading a fixed volume up network'?
#good: decay is appropriately at power-law rate, 
# so is this capturing nested/recursive/cumulative contribution more or less correctly?
#bad/weird: this is showing a weird thing, unclear what it means in flow units
# "if the flow at the outlet is 1000, then this is the independent portion contributed by each network position" 
z |>
  #as_tibble() |> 
  arrange(desc(in_km)) |> 
  mutate(
    in_km_pct = in_km / sum(in_km),
    q1k =  1000*in_km_pct,
    q1k_log10 = log10(q1k)
  ) |> 
  #arrange(q1k) |> filter(!is.infinite(q1k_log10))
  #pull(q1k) |> sum(na.rm = T)
  ggplot() + 
  geom_sf(data = st_as_sf(lsfn, "edges"), col = "lightblue", linewidth = 0.3) +
  geom_sf(aes(color = q1k_log10)) + wacolors::scale_color_wa_b("ferries", reverse = T)


#check summing upstream in subgraph
#this is somewhat reassuring: the 'outlet' flow is scaling appropriately within-network
lsfn <- lsfn |> 
  tidygraph::activate("nodes") |> 
  mutate(
    in_km_pct = in_km / sum(in_km),
    q1k =  1000*in_km_pct,
    q1k_log10 = log10(q1k)
  )

z <- lsfn |> 
  tidygraph::activate("nodes") |> 
  mutate(
    q_recalc = tidygraph::map_local_dbl(
      order = lsfn_order, mode = "in", 
      .f = function(neighborhood, ...){
        sum(as_tibble(neighborhood, active = 'nodes')$q1k, na.rm = T)
        }
      )
  )


{
as_tibble(z, active = 'nodes') |> #pull(q_recalc) |> summary()
  ggplot() + 
  geom_sf(data = st_as_sf(lsfn, "edges"), col = "lightblue", linewidth = 0.3) +
  geom_sf(aes(color = q_recalc), size = 0.5) + 
  wacolors::scale_color_wa_c("ferries", reverse = T) +
  labs(
    title = "Assuming instantaneous travel time, \n
       if outlet flow is 1000 cfs/cms/etc, and independent network contributions are uniform, \n
       then this is the (instant) value everywhere in the network, \n
       with Q per node summed from upstream subgraph (in-neighborhood)
    ")
}+{
as_tibble(z, active = 'nodes') |> 
  ggplot() + 
  geom_sf(data = st_as_sf(lsfn, "edges"), col = "lightblue", linewidth = 0.3) +
  geom_sf(aes(color = log10(q_recalc)), size = 0.5) + 
  wacolors::scale_color_wa_b("ferries", n.breaks = 6, reverse = T) +
  labs(title = "same, log10 to show power law decay")
} + plot_layout(ncol = 1)

ggsave(filename = "figures/f_ins_tt_unif_contrib_1k.png", width = 7, height = 11, dpi = 150, bg = "white")

```

## alt2a: alt1 + gradient

loop back to, but paused pending further review of [Brinkerhoff et al.](https://github.com/craigbrinkerhoff/CONUS_ephemeral)

# what about desc==culvert (& vs fpb)

> "The elevation derived hydrography specification calls for the segmentation and distinction of culvert features from the rest of the flowline network. Culvert features were automatically identified by comparing the monotonically enforced elevation values to the elevation ground model. Vertices that are misaligned with the DEM surface after monotonic Z smoothing and in close proximity to ancillary road and culvert data were extracted and used to classify culverts within the network. This segmentation process was manually reviewed. Additional culverts were manually extracted as necessary." -- NV5 Geospatial, pilot project report "Stillaguamish Watershed Lidar-derived 3DHP", June 9, 2023

"Connector: culvert" reaches (edges) in terms of the upstream network km

```{r func_gg_culv}
#conn_culvert_upstream

# as_tibble(lsfn, "edges") |> 
#   #as_tibble() |> count(desc)
#   #ggplot() + geom_sf(aes(color = desc)) + scale_color_manual(values = c("grey","orange","lightblue"))
#   ggplot() + geom_sf(aes(color = culvert)) + scale_color_manual(values = c("lightblue","orange"))
# 
# as_tibble(lsfn, active = "nodes") |> ggplot() + geom_sf(aes(color = below_culvert))
# 
# as_tibble(lsfn, "edges") |> st_zm() |> mapview::mapview(zcol = "culvert", burst = T)

# ggplot() +
#   geom_sf(data = as_tibble(lsfn, "edges") |> filter(!culvert), linewidth = 0.2, color = "lightblue") +
#   geom_sf(data = as_tibble(lsfn, "nodes") |> filter(below_culvert), aes(color = in_km, size = in_km)) + 
#   scale_color_gradient(low = "white", high = "purple") + scale_size_area(max_size = 4) +
#   theme(line = element_blank(), axis.text = element_blank()) +
#   ggspatial::annotation_scale(location = "bl")

gg_culv <- function(){
  ggplot() +
    geom_sf(data = as_tibble(lsfn, "edges") |> filter(!culvert), linewidth = 0.2, color = "lightblue") +
    geom_sf(data = as_tibble(lsfn, "nodes") |> filter(below_culvert),
            #aes(color = in_km, size = 1 / d_root_km)) + 
            aes(size = in_km, color = 1 / d_root_km)) + 
    #scale_color_gradient(low = "tan", high = "purple") + 
    wacolors::scale_color_wa_b("puget", n.breaks = 8, reverse = T) + 
    scale_size_area(max_size = 5) +
    theme(line = element_blank(), axis.text = element_blank()) +
    ggspatial::annotation_scale(location = "bl") +
    labs(subtitle = paste("3DHP for NHD comid", as_tibble(lsfn, "edges")$comid[1]))
}

#gg_culv()

# lsfn |>
#   mutate(
#     ac = tidygraph::centrality_alpha()
#   ) |> filter(in_n != ac)

```

```{r culv_eda}
#first couple cats not in 3dhp, train on 24275405 32+5
sf_nhdp_wa_3dhp_fpb[lines_3dhp_uch,] |> arrange(desc(yes+no))

fcat <- "24275405"
lsfn <- prep_sfn(fcat)
#can should instead reverse join poly comid on point...
lsfn_fpb <- sf_fpb_3dhp[filter(sf_nhdp_wa_3dhp_fpb, comid == fcat),]

gg_culv() +
  geom_sf(
    data = lsfn_fpb, aes(shape = fish_use), 
    color = "grey30", alpha = 0.5) + scale_shape_manual(values = c(3,4))

mapview::mapview(st_zm(as_tibble(lsfn, "edges")), zcol = "culvert", burst = T, color = c("lightblue","hotpink")) +
  mapview::mapview(lsfn_fpb, zcol = "fish_use", burst = T)

#compare field assessed slope for a match: uid 65142
#https://apps.wdfw.wa.gov/fishpassagephotos/Reports/LP27_Report.pdf
#slope 1.03%


# #remotes::install_github("ropensci/slopes")
# #from slope_matrix workhorse help "The output value is a proportion representing the change in elevation for a given change in horizontal movement along the linestring. 0.02, for example, represents a low gradient of 2% while 0.08 represents a steep gradient of 8%."
# z <- as_tibble(lsfn, "edges") |> filter(uid == "65142")
# 
# slopes::z_value(z)
# slopes::slope_xyz(z, directed = T)
# #need to confirm that in reaches/edges 
# #first point is upstream/max
# #last point is downstream/min
# z |> mutate(
#   across(SHAPE, list(
#     zslope_est = ~slopes::slope_xyz(., directed = T),
#     z1 = ~slopes::z_start(.),
#     z2 = ~slopes::z_end(.),
#     #zmin = ~slopes::z_min(.), #bug, calls max(), issue created
#     zmin = ~min(slopes::z_value(.)), 
#     zmax = ~slopes::z_max(.)
#   ))
# ) |> select(uid, contains("z"))
# 
# #cannot get to scatter against FPDSI slope without different version of those data (not going to review individual level A reports)
# #but can still start to look at (mis)match...
# #inflections at ~20ft buffer distance 
# tibble(d = seq(0, 50, by = 2)) |> 
#   mutate(
#     n_intsct = d |> 
#       map_int(
#         #st_intersection(
#         ~st_intersects(
#           as_tibble(lsfn, "edges") |> filter(culvert),
#           lsfn_fpb |> st_buffer(dist = .x)
#         ) |> unlist() |> length()
#       )
#   ) |> 
#   ggplot() + geom_col(aes(d, n_intsct))
# 
# d_buff <- 20
# 
# #if just inner_join wanted can use st_intersection
# #but this preserves non-matches as NA
# z <- st_join(
#   as_tibble(lsfn, "edges") |> filter(culvert),
#   lsfn_fpb |> st_buffer(dist = d_buff)
#   ) |> 
#   mutate(
#     across(SHAPE, list(
#       zslope_est = ~slopes::slope_xyz(., directed = T),
#       zmax = ~slopes::z_max(.)
#       ), .names = "{.fn}")
#   )
# 
# z |> 
#   ggplot(aes(fish_use, zslope_est)) +
#   geom_jitter(width = 0.2, alpha = 0.6) +
#   geom_violin(fill = NA)

lines_3dhp_culv <- lines_3dhp |> 
    filter(str_detect(Desc, "ulvert")) 

#full set? (FAST!?!)
sf_3dhp_culv_fpb <- st_join(
  lines_3dhp_culv
  ,
  sf_fpb_3dhp |> st_buffer(dist = d_buff)
  ) |> 
  mutate(
    across(SHAPE, list(
      zslope_est = ~slopes::slope_xyz(., directed = T),
      zmax = ~slopes::z_max(.)
      ), .names = "{.fn}")
  ) |> 
  select(
    uid = UniqueID, shape_length = SHAPE_Length, comid, desc = Desc,
    SiteId, FeatureType, fish_use, FishUseCriteriaCode, FPBarrierStatusCode, FPBarrierReasonCode, PercentFishPassableCode, SurveyDate,
    zslope_est, zmax
    )

sf_3dhp_culv_fpb |> #summary()
  filter(zslope_est <= 0) |>  #drop a single edge with a positive slope (assume this is artifact)
  ggplot(aes(fish_use, zslope_est)) +
  geom_jitter(width = 0.2, alpha = 0.6) +
  geom_violin(fill = NA)

#plot coerced to centroid point for better vis
{
ggplot() +
  geom_sf(
    data = sf_3dhp_culv_fpb |>
      filter(is.na(SiteId), zslope_est <= 0) |> #nrow()
      st_centroid(),
    aes(color = zslope_est), shape = 18, alpha = 0.6
    ) + 
    wacolors::scale_color_wa_c("foothills") +
    labs(subtitle = "8,332 3DHP without (20ft buffer) matches to fish passage inventory culverts")
}+{
ggplot() +
  geom_sf(
    data = sf_3dhp_culv_fpb |>
      drop_na(SiteId) |> 
      filter(zslope_est <= 0) |> #nrow() #drop a single edge with a positive slope (assume this is artifact)
      st_centroid(),
    aes(color = zslope_est, shape = fish_use), alpha = 0.6
    ) +
  wacolors::scale_color_wa_c("foothills") +
  scale_shape_manual(values = c(19,10,15)) +
    labs(subtitle = "462 3DHP with (20ft buffer) matches")
} + {
  ggplot() +
    geom_sf(
      data = sf_fpb_3dhp |> 
        anti_join(as_tibble(sf_3dhp_culv_fpb) |> drop_na(SiteId), by = "SiteId")
      ,
      aes(shape = fish_use), alpha = 0.6
    ) +
    geom_sf(data = lines_3dhp_uch, fill = NA, color = "hotpink") +
    scale_shape_manual(values = c(19,10,15)) +
    labs(subtitle = "915 fish passage inventory culverts without 3DHP barrier matches \n (including those within convex hull but outside current 3DHP linework extent)")
} +
  plot_layout(ncol = 1) +
  plot_annotation(
    title = "Estimated slopes for 3DHP 'connector: culvert' reaches vs fish passage inventory culverts",
    subtitle = "Approx. ~5% of 3DHP 'connector: culvert' (462 of 8794) with barrier matches")

ggsave(filename = "figures/f_culv_slope_compare.png", width = 7, height = 11, dpi = 150, bg = "white")



#examine survey date/source for nonmatches/far distances
sf_fpb_3dhp |> 
  as_tibble() |> 
  #semi_join(as_tibble(sf_3dhp_culv_fpb), by = "SiteId") |> 
  anti_join(as_tibble(sf_3dhp_culv_fpb), by = "SiteId") |> 
  count(DataSource) |> arrange(desc(n))


#Euclidean distances from all fpb to any/culvert 3dhp lines
#takes a minute, returns n-barriers by n-reaches matrix
#need/want to integrate st_nearest_feature?
d_fpb_culv <- st_distance(
  sf_fpb_3dhp,  #barrier points
  lines_3dhp_culv, #culvert line sections only
)

#per barrier, along rows 
sf_fpb_3dhp$d_culv <- round(apply(d_fpb_culv, 1, min), 2)
sum(sf_fpb_3dhp$d_culv <= 20) #matches st_intersect on st_buffer(20) in st_join
plot.ecdf(sf_fpb_3dhp$d_culv)
#worst issues appear due to using convex_hull to constrain barrier set
sf_fpb_3dhp |> 
  mutate(d_in20 = d_culv <= 20) |> 
  #filter(d_culv < 1000) |> 
  ggplot() + geom_sf(aes(color = log10(d_culv))) + facet_wrap(~d_in20, ncol = 1)
#can also test against full linework set for barriers that are on/close to 'stream' reaches (missing in 3DHP?)

#distributions of distances to nearest culvert line by source
sf_fpb_3dhp |> 
  filter(d_culv < 1000) |> 
  ggplot() + 
  geom_boxplot(aes(d_culv, DataSource), varwidth = T)

#unique(nchar(sf_fpb_3dhp$SurveyDate))
sf_fpb_3dhp$date = as.Date(sf_fpb_3dhp$SurveyDate, format = "%m/%d/%Y")
#sf_fpb_3dhp |> filter(is.na(date)) |> select(SiteId, DataSource, SurveyDate)

sf_fpb_3dhp |> 
  ggplot() +
  geom_point(aes(year(date), d_culv, color = DataSource))


```


```{r needs_work}
lines_3dhp_culv <- lines_3dhp |> 
    filter(str_detect(Desc, "ulvert")) 

as_tibble(lines_3dhp_culv) |> #summary()
  rename(length = SHAPE_Length) |> 
  ggplot() + 
  geom_point(aes(length, zmax, color = abs(zslope_est)), alpha = 0.3) +
  scale_x_log10() + scale_y_log10()

#try the whole thing
#goal would be to supplement above ordination of culvert reaches by elevation and slope with up/downstream network distances

#pretty fast, "directed acyclic multigraph with 22 components"
lines_3dhp_sfn <- lines_3dhp |> 
  sfnetworks::as_sfnetwork() 

lines_3dhp_sfn <- lines_3dhp_sfn |>
  tidygraph::activate(nodes) |>
  mutate(ng = tidygraph::group_components())
#non-primary components appear appropriately peripheral?
as_tibble(lines_3dhp_sfn, "nodes") |> 
  #as_tibble() |> count(ng) |> arrange(desc(n)) |> print(n=Inf)
  filter(ng > 1) |> mutate(ng = letters[ng]) |>
  #ggplot() + geom_sf(aes(color = ng))
  st_zm() |> mapview::mapview(zcol = "ng", burst = T)

#this results in 'directed acyclic multigraph with 1 component'
lines_3dhp_sfn <- lines_3dhp_sfn |> 
  tidygraph::convert(
    tidygraph::to_largest_component,
    .clean = T
  ) 
#need to reduce (multiedges) to tree to calc node_distance_to
tidygraph::with_graph(lines_3dhp_sfn, tidygraph::graph_is_tree())
# #have tidygraph option and sfnetwork extension, latter seems preferable for control on 'attribute_summary' arg
# #but other tidygraph 'to_*_tree' options 
# lines_3dhp_sfn |> tidygraph::convert(tidygraph::to_simple)
lines_3dhp_sfn <- lines_3dhp_sfn |> tidygraph::convert(sfnetworks::to_spatial_simple, .clean = T)

#also have https://luukvdmeer.github.io/sfnetworks/reference/simplify_network.html
#but looks like needs package update? sfnetworks::simplify_network(lines_3dhp_sfn)

#but still not a tree
tidygraph::with_graph(lines_3dhp_sfn, tidygraph::graph_is_tree()) 
tidygraph::with_graph(lines_3dhp_sfn, tidygraph::node_is_root()) |> any() 
lines_3dhp_sfn |> filter(tidygraph::node_is_root())

#may need to work with dominator tree? https://tidygraph.data-imaginist.com/reference/node_topology.html
#however, there is an edge with a lowest elevation Z-val...
as_tibble(lines_3dhp_sfn, "edges") |> as_tibble() |> arrange(zmax) 
#the 'to' node is 73095 and/but it has no 'from'
as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 73095)
#looks like maybe still some 'extra' edges? lsfn has 1 more node than edges, this has 41
lines_3dhp_sfn |> 
  tidygraph::activate("nodes") |> 
  mutate(idn = seq_along(SHAPE)) |> 
  arrange(tidygraph::node_topo_order())
#okay so a small number of 'from' nodes appear twice, likely the issue?
as_tibble(lines_3dhp_sfn, "edges") |> as_tibble() |> 
  #count(to) |> arrange(desc(n)) |> count(n) #all either 1 or 2
  count(from) |> arrange(desc(n)) |> 
  #count(n)
  filter(n > 1)
#yes, have same node going to 2 different other nodes
#and have a split path
bind_rows(
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 27590)
  ,
  as_tibble(lines_3dhp_sfn, "edges") |> filter(to == 27590)
  ,
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 25192), 
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 13639),
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 8703),
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 8704)
  ,
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 28273)
) |> 
  select(-COMID) |> st_zm() |> mapview::mapview()

#there are non-junction 'in-line' nodes - indeed, many are 'connector: culvert'
#which explains the count(to) n==1
bind_rows(
  as_tibble(lines_3dhp_sfn, "edges") |> filter(to == 1)
  ,
  as_tibble(lines_3dhp_sfn, "edges") |> filter(from == 1)
  ,
  as_tibble(lines_3dhp_sfn, "edges") |> filter(to == 341),
  as_tibble(lines_3dhp_sfn, "edges") |> filter(to == 340) #from is 25530, which appears to be leaf
  #no rows/edges: as_tibble(lines_3dhp_sfn, "edges") |> filter(to == 25530)
) |> 
  select(-COMID) |> st_zm() |> mapview::mapview()

#for each of the nodes with multiple 'from' drop one of the split paths?
as_tibble(lines_3dhp_sfn, "edges") |> as_tibble() |> 
  semi_join(
    as_tibble(lines_3dhp_sfn, "edges") |> as_tibble() |> count(from) |> arrange(desc(n)) |> filter(n > 1), 
    by = "from"
  )

#could try defining 'root' by lowest elevation node (73095)
#tidygraph::bfs_dist and similar need a root defined...
lines_3dhp_sfn |> 
  mutate(depth = tidygraph::bfs_dist(root = 73095, mode = "in")) |> as_tibble("nodes") |> as_tibble() |> count(depth) |> print(n = 40)

lines_3dhp_sfn <- lines_3dhp_sfn |>  
    mutate(
      idn = seq_along(SHAPE),
      #raindrop perspective paths from nodes to root
      #d_root = tidygraph::node_distance_to(tidygraph::node_is_root()),
      d_root = tidygraph::node_distance_to(73095),
      #d_root_ft = tidygraph::node_distance_to(tidygraph::node_is_root(), weights = shape_length),
      #d_root_km = 0.0003 * d_root_ft,
      #looking 'up' perspective to perimeter from per-node 'me'
      in_n = tidygraph::local_size(order = nrow(tidygraph::.N()), mode = "in")
    )

as_tibble(lines_3dhp_sfn, "nodes") |> arrange(desc(d_root))
#illustrates complexity of lower system distributary/slough structure
as_tibble(lines_3dhp_sfn, "nodes") |> arrange(d_root) |> slice(1:10) |> st_zm() |> mapview::mapview()
as_tibble(lines_3dhp_sfn, "edges") |> select(-COMID) |>
  filter(to %in% {as_tibble(lines_3dhp_sfn, "nodes") |> arrange(d_root) |> slice(1:50) |> pull(idn)}
  ) |> 
  st_zm() |> mapview::mapview()
```


```{r predef_rooted_tree}
#dev 'tree-based' workflow from a priori determined root
#acknowledging that split paths and divergences are biogeo correct
#but need categorically distinct handling?

#def root 
#manual via QGIS inspection; start with Pilchuck Creek
root_edge_uid <- 73853
# #node id in of root_edge in original, pre-convert indexing
# filter(as_tibble(sfn, "edges"), uid == root_edge_uid)$to

#need separate declarations to assure node ID in to_local_neighborhood
#and less importantly the order to search over
#also have sfnetworks::to_spatial_neighborhood,
  #node = filter(as_tibble(sfn, "edges"), uid == root_edge_uid)$to,
  #threshold = , #in sf units; untested whether it would work to use edge attribute as: weights = length_ft

sfn <- lines_3dhp |>
  select(comid, desc = Desc, uid = UniqueID, GNIS_Name, 
         length_ft = SHAPE_Length, zmax, zslope_est) |> 
  sfnetworks::as_sfnetwork() 

sfn <- sfn |> 
  tidygraph::convert(
    tidygraph::to_local_neighborhood,
    node = filter(as_tibble(sfn, "edges"), uid == root_edge_uid)$to,
    mode = "in", 
    order = tidygraph::with_graph(sfn, tidygraph::graph_order())
    ,
    .clean = T
  )
# #already 'directed acyclic simple graph with 1 component'
# #seems unlikely that 3DHP would have loops/multiples?
# #can add back to pipe if seems needed for certainty
# tidygraph::convert(sfnetworks::to_spatial_simple, .clean = T) 

# #however, Pilchuck has a valid split path at node 1033
# as_tibble(sfn, "edges") |> as_tibble() |> count(from) |> arrange(desc(n)) |> filter(n > 1)
# list(
#   edge = as_tibble(sfn, "edges") |> st_zm(),
#   node = as_tibble(sfn, "nodes") |> st_zm() |> rowid_to_column("idn") |> mutate(issplit = if_else(idn==1033,T,F))
#   ) |> mapview::mapview(zcol = list(NULL, "issplit"))

#this still works but now has a different index value (from 73858 to 7131) due to `convert`
root_node_idn <- filter(as_tibble(sfn, "edges"), uid == root_edge_uid)$to

sfn <- sfn |> 
  tidygraph::activate("nodes") |> 
  mutate(
    idn = seq_along(SHAPE),
    #raindrop perspective paths from nodes to root
    d_root = tidygraph::node_distance_to(root_node_idn),
    d_root_ft = tidygraph::node_distance_to(root_node_idn, weights = length_ft),
    d_root_km = 0.0003 * d_root_ft,
    #looking 'up' perspective to perimeter from per-node 'me'
    in_n = tidygraph::local_size(order = tidygraph::graph_order(), mode = "in", mindist = 1)
    ) |> 
  arrange(d_root) |> 
#as in 'prep_sfn'
#associate the immediately upstream streamlength from edges to each 'to' node
#summed over both upstream tribs; also add 'below_culvert' logical:
# defining per-node state in {no-incoming, not, culvert, not+not, culvert+not, culvert+culvert}
  left_join(
    select(as_tibble(as_tibble(sfn, active = "edges")), uid, idn = to, length_ft, desc, zslope_est, zmax) |>
      summarise(
        node_below_culvert = if_else(any(str_detect(desc, "ulvert")), T, F),
        in_slope_max = max(zslope_est),
        in_ft = sum(length_ft),
        .by = "idn"),
    by = "idn"
  )
#wait to add 'above' neighborhood length until pruned by slope

#2421 NA node_below_culvert at upstream end of perimeter edges
as_tibble(as_tibble(sfn, "nodes")) |> count(node_below_culvert) 


#prune periphery by slope
#interesting options to consider. 
#more efficient to begin from root and breadth-first-search up, stopping at any reach over a gradient threshold?

plot.ecdf(as_tibble(sfn, "nodes")$in_slope_max)
slope_threshold <- -0.3

#some slopes >100%, but actually appear to be valid (e.g., ) 
list(
  edge = as_tibble(sfn, "edges") |> st_zm(),
  node = as_tibble(sfn, "nodes") |> st_zm()
  ) |> 
  #mapview::mapview(zcol = list("length_ft", "in_slope_max"))
  mapview::mapview(zcol = list("length_ft", "node_below_culvert"))

# #480 ft, 2455-1969=486ft drop, so actually the >100% estimate is not off...
# #and some individual point-to-points are well over that
# m <- filter(lines_3dhp, UniqueID==69444) |> pull(SHAPE) |> st_coordinates() 
# slopes::slope_matrix_weighted(m, lonlat = F)
# slopes::slope_matrix(m, lonlat = F)
# (p <- filter(lines_3dhp, UniqueID==69444) |> pull(SHAPE) |> st_cast("POINT") |> st_as_sf() |> rename(g = x) |>  
#   mutate(
#     g2 = lead(g),
#     grnd_dist = st_distance(g, g2, by_element = T),
#     e1 = st_coordinates(g)[,"Z"],
#     e2 = st_coordinates(g2)[,"Z"],
#     elev_drop = e1 - e2,
#     slp = elev_drop / grnd_dist
#     ))

# #ignoring routedness 'strands' some lower gradient components
# #and mostly trims leaf nodes with NA 'in_slope_max'
# x <- sfn |> 
#   tidygraph::activate("nodes") |> 
#   filter(in_slope_max > slope_threshold) 
# 
# list(
#   trim = as_tibble(x, "edges") |> st_zm(),
#   full = as_tibble(sfn, "edges") |> st_zm()
# ) |> mapview::mapview(color = c("orange","blue"))
# 
# #cleaner to filter out edges then just reduce to the component that includes root?
# y <- sfn |> 
#   tidygraph::activate("edges") |> 
#   filter(zslope_est > slope_threshold) 
# #yes. better. still gets rid of steep fringes and/but creates 767 components
# yy <- y |> tidygraph::convert(
#   tidygraph::to_largest_component, .clean = T
#   )
# yy |> filter(uid == root_edge_uid)
# 
# list(
#   yy = as_tibble(yy, "edges") |> st_zm(),
#   y = as_tibble(y, "edges") |> st_zm(),
#   full = as_tibble(sfn, "edges") |> st_zm()
# ) |> mapview::mapview(color = c("purple","green","blue"))

sfn <- sfn |>
  tidygraph::activate("edges") |>
  filter(zslope_est > slope_threshold) |> #greater than a negative value
  tidygraph::convert(tidygraph::to_largest_component, .clean = T) |> 
  tidygraph::activate("nodes") |> 
  mutate(
    idn = seq_along(SHAPE),
    in_ft_all = tidygraph::map_local_dbl(
      order = tidygraph::graph_order(), 
      mode = "in",
      .f = function(neighborhood, ...){
        sum(as_tibble(neighborhood, active = 'nodes')$in_ft, na.rm = T)
      }
    ),
    in_km_all = 0.0003 * in_ft_all,
    in_ft_pct = in_ft_all / max(in_ft_all)
  )


## ordinate culverts by root dist & upstream extent & slope & length
## slope_threshold filter creates some stranded perimeter 'below_culv' nodes
anti_join(
  as_tibble(as_tibble(sfn, "nodes")) |> filter(node_below_culvert)  #983
  ,
  as_tibble(as_tibble(sfn, active = "edges")) |> filter(str_detect(desc, "ulvert")) #955
  ,
  by = c("idn" = "to")
) #31

as_tibble(as_tibble(sfn, "nodes")) |> count(node_below_culvert)
#confirming that the 31 TRUE below_culvert nodes are at perimeter having had their upstream edges, at least one of which was culvert, cut off
list(
  edge = as_tibble(sfn, "edges") |> st_zm(),
  node = as_tibble(sfn, "nodes") |> st_zm(),
  # foo = anti_join(
  #   as_tibble(sfn, "nodes") |> filter(node_below_culvert),
  #   as_tibble(as_tibble(sfn, active = "edges")) |> filter(str_detect(desc, "ulvert")),
  #   by = c("idn" = "to")
  # ) |> st_zm()
  foo = as_tibble(sfn, "nodes") |> slice(415) |> st_zm()
) |>
  mapview::mapview(zcol = list("length_ft", "node_below_culvert", "d_root"))


#if topology attributes added back to edges, could also have edge/reach/culvert length
#not associated to nodes themselves 
#maybe better to overhaul/re-do treating edges as nodes? sfn |> tidygraph::convert(tidygraph::to_linegraph) 

sfn_culv <- as_tibble(sfn, active = "edges") |> 
  filter(str_detect(desc, "ulvert")) |> 
  left_join(
    as_tibble(as_tibble(sfn, "nodes")) |> 
      select(to = idn, 
             #node_below_culvert, 
             starts_with("d_"), 
             in_ft_all, in_km_all, in_ft_pct #starts_with("in_")
             )
    ,
    by = "to"
  )

sfn_culv |> st_zm() |> st_centroid() |> 
  mapview::mapview(zcol = "in_km_all")

#interesting/useful, and/but note the outlier high in_km is uid 64871 Lake Creek with Lk Cavanaugh and all its tribs above
#such that it could indeed be a critical culvert...or not, depending on expectations for through lake mvmt
plotly::ggplotly(
as_tibble(sfn_culv) |> 
  ggplot() + 
  #geom_point(aes(d_root_km, zslope_est, color = in_km_all, size = in_km_all)) +
  #geom_point(aes(d_root_km, in_km_all, color = abs(zslope_est), size = abs(zslope_est))) +
  geom_point(aes(d_root_km, in_km_all, color = abs(zslope_est), size = length_ft), alpha = 0.5) 
#+ scale_size_binned_area(max_size = 4) + scale_color_binned()
)

#need some 'marginal habitat gain' for optimization
#calc distance to next upstream culvert (or end of line if none)?
#could do all-to-all then determine nearest then use that to index out of all-all object
#this works, https://luukvdmeer.github.io/sfnetworks/reference/st_network_cost.html

#direction "in" defines distances from upstream nodes
#for first row (root node) equivalent to d_root_ft
d_in <- sfnetworks::st_network_cost(sfn, direction = 'in', weights = as_tibble(sfn, "edges")$length_ft)
d_in[1:10,1:10]
#per-row, the min distance is always 0 to self, the second (next largest in sort) is the nearest upstream
#so when subset to the "to-node" of culvert edges becomes the closest upstream culvert, possibly Inf if none above
#however need logic to use 'in_ft_all' if d_in==Inf since 'final' culverts may have lots of upstream length
apply(d_in[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])
which(sfn_culv$to==10) #index 415 is node 10: sfn_culv$to[415]
apply(d_in[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])[415] #9104 ~ 69+2793+267+975+2218+1839+76+197+668
apply(d_in[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])[which(sfn_culv$to==76)] #466 ~ 33.25+267+165.6
apply(d_in[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])[which(sfn_culv$to==121)] #207 ~ 66.3+140.8
apply(d_in[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])[which(sfn_culv$to==205)] #Inf as it should be with no higher culverts

#direction "out" defines distances downstream to nodes below
#first row (root/outlet) is all Inf other than 0 on diag
#the per-row min is again always 0 to self but now the next largest is distance to nearest downstream
#and when subset to culvert edge "to-node" means that Inf values are paths with no downstream culverts before the root
#which is quite interesting for highlighting 'already open' paths or equivalently the lowest culvert in any subbasin!
d_out <- sfnetworks::st_network_cost(sfn, direction = 'out', weights = as_tibble(sfn, "edges")$length_ft)
d_out[1:10,1:10]
apply(d_out[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])

#trying to get the subgraph of shortest paths from (lowest) culverts with none below to root
lowest_culv <- sfn_culv$to[which(apply(d_out[sfn_culv$to, sfn_culv$to], 1, \(x) sort(x)[2])==Inf)]
lowest_culv_paths <- sfnetworks::st_network_paths(sfn, from = 1, to = lowest_culv, mode = 'in')

sfn_lowest_culv <- sfn |> 
  tidygraph::activate("edges") |> 
  slice(unique(unlist(lowest_culv_paths$edge_paths))) |> 
  tidygraph::convert(
    tidygraph::to_largest_component
  ) 

ggplot() +
  geom_sf(data = as_tibble(sfn, "edges"), color = "lightblue") +
  geom_sf(data = as_tibble(sfn_lowest_culv, "edges"), color = "darkgreen")

list(
  edge = as_tibble(sfn, "edges") |> st_zm(),
  node = as_tibble(sfn, "nodes") |> st_zm(),
  lowculv = as_tibble(sfn_lowest_culv, "edges") |> st_zm()
) |>
  mapview::mapview(zcol = list("length_ft", "node_below_culvert", "length_ft"))

# #something more elegant? https://luukvdmeer.github.io/sfnetworks/articles/sfn05_routing.html
# #maybe riffing off https://stackoverflow.com/questions/76115139/create-subgraph-from-shortest-path-calculation-sfnetworks-r
# #could not get tidygraph::node_distance_from mapping correctly
# #this is correct but slow as-is, likely easy speed up at least from declaring weight vector standalone
# sfn |> 
#   mutate(
#     d_culv_up = map_dbl(
#       idn,
#       ~sort(sfnetworks::st_network_cost(sfn, from = sfn_culv$to, to = .x, direction = 'out', weights = as_tibble(sfn, "edges")$length_ft))[2]
#     )
#   ) 

#for now use the distance to the nearest upstream culvert or all upstream length if no culverts above
#non-culvert nodes need first position of sort: apply(d_in[, sfn_culv$to], 1, \(x) sort(x)[1])
sfn_culv <- sfn_culv |> 
  mutate(
    d_culv_up = apply(d_in[to, to], 1, \(x) sort(x)[2])
    ,
    hmarg = if_else(is.infinite(d_culv_up), in_ft_all, d_culv_up)
    ,
    psuedocost = abs(length_ft * as.vector(zslope_est))
  )

```

try for something derived from  `solve_opt` from https://github.com/sljardine/upstream/blob/main/upstream/R/fct_Suggest.R

```{r try_solve_opt}

# # objective function
# obj <- ROI::L_objective(v)
# # constraints
# bc <- ROI::L_constraint(L = brc, dir = "<=", rhs = budget) #budget
# hc <- ROI::L_constraint(L = diag(nb) - t(D), dir = rep("<=", nb), rhs = 1 - di) #hydrology
# # problem & solution
# # bc and hc
# prob <- ROI::OP(objective = obj,
#                 constraints = c(bc, hc),
#                 bounds = ROI::V_bound(li = 1 : nb, lb = rep.int(0, nb), ui = 1 : nb, ub = rep.int(1, nb)),
#                 types = rep.int("B", nb),
#                 maximum = TRUE)
# 
# sol <- ROI::ROI_solve(prob, "glpk",  control = list("verbose" = TRUE, "presolve" = TRUE))
# soln <- as.logical(ROI::solution(sol))

#start with unweighted hab: obj==1 s.t. points$hmarg -> h -> v, 
#with points sf df with: "hmarg_length, hmarg_area, hmarg_volume, cost, and wria_number"
#inspecting D.rda from repo, 'connectivity matrix' means binary but need to follow up with Sunny on exactly how generated
#here, 1 if culvert-i (row) has a non-infinite downstream distance to culvert-j (col)
#note need: install.packages("ROI.plugin.glpk")
library(ROI.plugin.glpk)
ROI::ROI_registered_solvers()
#so this runs but 'sol' is allNA and throws "PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION"
v <- sfn_culv$hmarg
D <- 1*!is.infinite(d_out[sfn_culv$to, sfn_culv$to])
brc <- sfn_culv$psuedocost
nb <- length(v)
di <- colSums(D)
budget <- 20 #plot(density(brc)); quantile(brc, p = seq(0,1, by = 0.1))
obj <- ROI::L_objective(v)
bc <- ROI::L_constraint(L = brc, dir = "<=", rhs = budget) #budget
hc <- ROI::L_constraint(L = diag(nb) - t(D), dir = rep("<=", nb), rhs = 1 - di) #hydrology
prob <- ROI::OP(objective = obj,
                constraints = c(bc, hc),
                bounds = ROI::V_bound(li = 1 : nb, lb = rep.int(0, nb), ui = 1 : nb, ub = rep.int(1, nb)),
                types = rep.int("B", nb),
                maximum = TRUE)
sol <- ROI::ROI_solve(prob, solver = "glpk",  control = list("verbose" = TRUE, "presolve" = TRUE))
soln <- as.logical(ROI::solution(sol))

```



fit for 'prob of being a real culvert'?

split the basin, train on fpb matches, test/validate in remaining?

much more to think about doing with `LinealGainMeasurement`, `SignificantReachCode`, etc.

run version of climate adapted, given slope and some DA/contributing channel length regression on width?


# examine Brinkerhoff et al. 2024

```{r first_look}
#from https://apps.nationalmap.gov/downloader/ via e.g. 'https://prd-tnm.s3.amazonaws.com/StagedProducts/Hydrography/NHDPlusHR/VPU/Current/GDB/NHDPLUS_H_1703_HU4_GDB.zip'

#1708 is Lwr Col; IMW basins are evident in HUC12s
st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1708_HU4_GDB.gdb/"), layer = "WBDHU12") |> ggplot() + geom_sf()

#1703 Yakima
#1708 Lwr Col, ~200K lines
#1711 Puget Sound, ~300K lines 
st_layers(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"))
sf_nhdphr_h12 <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "WBDHU12")
sf_nhdphr_h10 <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "WBDHU10")
sf_nhdphr_h08 <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "WBDHU8")

sf_nhdphr_flw <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "NHDFlowline") |> 
  left_join(
    read_csv(file.path(dir_data_common, "nhdphr/brinkerhoffetal2024_results_archive/results_1711.csv"))
    ,
    by = "NHDPlusID") |> 
  drop_na(percQEph_reach)  

# # drop_na based on check of object without it
# # ~9.5k flowlines missing percent ephemeral 
# # mostly shorelines/canals etc.
# sf_nhdphr_flw |> filter(is.na(percQEph_reach)) |> ggplot() + geom_sf()

#cut down to Deschutes R
aoi <- sf_nhdphr_h12 |> 
  #as_tibble() |> 
  filter(States == "WA") |> 
  select(HUC12, Name, AreaSqKm) |> 
  #filter(str_detect(Name, "Still"))
  filter(str_detect(Name, "Deschut")) #|> mapview::mapview()

aoi_flw <- sf_nhdphr_flw |> 
  #as_tibble() |> #glimpse()
  mutate(h8 = str_sub(ReachCode, 1, 8)) |> 
  filter(h8 == unique(str_sub(aoi$HUC12, 1,8))) |> 
  st_cast("LINESTRING") 

# 'percQEph_reach'

aoi_flw |> 
  as_tibble() |> 
  ggplot() + geom_boxplot(aes(perenniality, percQEph_reach), varwidth = T, outliers = F)

mapview::mapview(aoi) +
  mapview::mapview(aoi_flw, zcol = "percQEph_reach") 

#as suggested in paper & supp, examining median at smaller spatial scale
#drives towards bimodality
sf_nhdphr_h08 |> 
  left_join(
    sf_nhdphr_flw |> 
      as_tibble() |>  
      mutate(HUC8 = str_sub(ReachCode, 1, 8)) |> #distinct(HUC8) #21, matches poly count
      summarise(
        across(percQEph_reach, list(median = ~median(., na.rm = T))),
        .by = HUC8
      ) #|> print(n= Inf)
    , by = "HUC8"
  ) |> 
  ggplot() + geom_sf(aes(fill = percQEph_reach_median)) + wacolors::scale_fill_wa_c()

#to vis at H10 or H12 seems need to do spatial join...
#117 in polys but 40 derived from ReachCode and only 19 in common
inner_join(
  as_tibble(sf_nhdphr_h10) |> distinct(HUC10) |> arrange()
  ,
  as_tibble(sf_nhdphr_flw) |> mutate(HUC10 = str_sub(ReachCode, 1, 10)) |> distinct(HUC10) |> arrange()
  , by = "HUC10")

#the ~2.3K 'huc12' values derived from ReachCode line up poorly with the ~511 in the HUC12 poly
#only 87 in common; whatever the first 12 positions are, they're not the 'normal' HUC12
inner_join(
  as_tibble(sf_nhdphr_h12) |> distinct(HUC12)
  ,
  as_tibble(sf_nhdphr_flw) |> mutate(HUC12 = str_sub(ReachCode, 1, 12)) |> distinct(HUC12)
  , by = "HUC12")


#table of major rivers' percent ephemeral?
#also have NHDPlusFlowlineVAA
erom <- st_read(file.path(dir_data_common, "nhdphr/17/NHDPLUS_H_1711_HU4_GDB.gdb/"), layer = "NHDPlusEROMMA")
glimpse(erom)

#is this a fascinating plot? 
#suggests that below an average QAMA about ~1 unit, systems are either mostly E or not (bimodal 1/0)  
#but that the larger systems are 'well-blended'
#i think this is implied if not directly stated in the paper w/some of the remarks on the 'scale of assessment' and the choice of HUC4 
sf_nhdphr_flw |> 
  as_tibble() |>  
  drop_na(GNIS_Name) |> #down to ~51k 'named' flowlines
  #count(GNIS_Name) |> arrange(desc(n)) #1700 named reaches
  left_join(erom |> select(NHDPlusID, QAMA), by = "NHDPlusID") |> 
  summarise(
    across(
      c(QAMA, percQEph_reach), 
      list(median = ~median(., na.rm = T))),
    .by = "GNIS_Name"
  ) |> 
  arrange(desc(QAMA_median)) |> 
  #print(n= Inf)  
  ggplot() + 
  geom_point(aes(QAMA_median, percQEph_reach_median), size = 0.5, alpha = 0.5) +
  scale_x_log10()


sf_nhdphr_flw |> 
  #as_tibble() |> 
  select(GNIS_Name, NHDPlusID, perenniality, percQEph_reach) |> 
  drop_na(GNIS_Name) |> #down to ~51k 'named' flowlines
  left_join(erom |> select(NHDPlusID, QAMA), by = "NHDPlusID") |> 
  #pull(QAMA) |> plot.ecdf(); abline(v = 100)
  #filter(QAMA > 100) |> #100 is ~12.2K lines 
  filter(QAMA > 50) |> #100 is ~17.2K lines 
  # summarise(
  #   across(
  #     c(QAMA, percQEph_reach), 
  #     list(median = ~median(., na.rm = T))),
  #   .by = "GNIS_Name"
  # ) |> 
  # ggplot() + geom_point(aes(QAMA_median, percQEph_reach_median), size = 0.5, alpha = 0.5) + scale_x_log10()
  #slice_max(QAMA, n = 1, by = GNIS_Name) |> 
  ggplot() +
  geom_sf(data = sf_nhdphr_h10, fill = NA) +
  #geom_sf(aes(color = percQEph_reach)) + wacolors::scale_color_wa_c("vantage")
  geom_sf(aes(color = percQEph_reach, linewidth = log10(QAMA))) + wacolors::scale_color_wa_c("vantage")

```

```{r sf_dnr_hydro}
#1.3m features
sf_dnr_hydro <- st_read(file.path(dir_data_common, "wadnr_state_hydro/wchydro.gdb"), layer = "wchydro")

# sf_dnr_hydro[st_as_sfc(st_bbox(lines_3dhp)),] #31.4k
sf_dnr_hydro_3dhp <- sf_dnr_hydro[lines_3dhp_uch,] #23.9k

```

```{r nhdphr_dnr}
#NHDplus med res catchment for big jim creek, used above in default call to prep_sfn
sf_cat <- filter(sf_nhdp_wa_3dhp_fpb, comid == "24275471")


# list(
#   nhdphr = sf_nhdphr_flw_3dhp[sf_cat,] |> select(GNIS_Name, LengthKM, perenniality, percQEph_reach)  
#   ,
#   dnr = sf_dnr_hydro_3dhp[sf_cat,]
#   ,
#   `3dhp` = as_tibble(lsfn, "edges") |> st_zm()
# ) |> 
#   mapview::mapview(
#     zcol = list("perenniality", "WC_PERIOD_LABEL_NM", "desc"),
#     burst = T,
#     color = list("blue","green", "orange")
#     )


mapview::mapview(
  sf_nhdphr_flw_3dhp[sf_cat,] |> select(GNIS_Name, LengthKM, perenniality, percQEph_reach)
  ,
  zcol = "perenniality", burst = T, color = c("lightblue","darkblue") #ephemeral, non-ephem
  ) +
mapview::mapview(
  sf_dnr_hydro_3dhp[sf_cat,]
  ,
  zcol = "WC_PERIOD_LABEL_NM", burst = T, color = c("#831285","#C886E8") #Perennial, unknown
) +
mapview::mapview(
  as_tibble(lsfn, "edges") |> st_zm()
  ,
  zcol = "desc", burst = T, color = c("yellow","orange","#386276")
)


# #not yet run  
# sf_nhdphr_dnr <- st_join(
#   sf_nhdphr_flw_3dhp[sf_cat,] |> select(GNIS_Name, LengthKM, perenniality, percQEph_reach)  
#   ,
#   sf_dnr_hydro_3dhp[sf_cat,]
#   ,
#   join = st_intersects, 
#   left = T, #F for inner_
#   largest = T
#   )


union_3dhp <- st_zm(as_tibble(lsfn, "edges")) |> st_union()
union_dnr <- sf_dnr_hydro_3dhp[sf_cat,] |> st_union()
union_diff <- st_difference(
  union_3dhp,
  union_dnr |> st_buffer(dist = 30)
)

union_3dhp |> #163km
 union_dnr |> #54.7km
 union_diff |> #123km outside 30ft buffer
  st_length() |> units::set_units("km") |> as.numeric()


{ ggplot() +  geom_sf(data = union_dnr, color = "gold") + labs(title = "DNR Hydrography") } +
{ ggplot() + geom_sf(data = union_3dhp, color = "lightblue") + labs(title = "3DHP flowlines")} +
{
ggplot() +
  geom_sf(data = union_3dhp, color = "darkgreen") +
  geom_sf(data = union_diff, color = "lightblue") +
  labs(title = "Intersection (green) and difference (blue) between \n 3DHP flowlines and a 30ft buffer on DNR Hydro")
} + 
  plot_layout(ncol = 1)

ggsave(filename = "figures/f_bigjim_overlay.png", width = 7, height = 11, dpi = 150, bg = "white")

```

# testing SSNbler & SSN2

next thing here is modeling APPS, perhaps conditioning on HRCD?

```{r compare_ssnbler_vs_lsfn}

# SSNbler::copy_streams_to_temp()
# path <- paste0(tempdir(), "/streamsdata")
# MF_streams <- st_read(paste0(path, "/MF_streams.gpkg"))
# MF_obs <- st_read(paste0(path, "/MF_obs.gpkg"))
# MF_pred1km <- st_read(paste0(path, "/MF_pred1km.gpkg"))
# MF_CapeHorn <- st_read(paste0(path, "/MF_CapeHorn.gpkg"))
# MF_streams

path_ssn <- "~/T/DFW-Team WDFW Watershed Synthesis - General/hydrography/test_ssnbler"
ssn <- list()
#test with default lsfn from prep_sfn("24275471") of 1505 reaches
ssn$edges <-SSNbler::lines_to_lsn(
  streams = as_tibble(lsfn, "edges"),
  lsn_path = path_ssn,
  overwrite = T
)

ssn$edges <- SSNbler::updist_edges(
  ssn$edges,
  lsn_path = path_ssn,
  #calc_length = T
  length_col = "shape_length"
)

ssn$edges #looks to match d_root_ft at least to reasonable precision

#vignette example uses presupplied drainage area "h2oAreaKm2"
#this may be a 3DHP attribute, but for now testing with per-node upstream in-length 
#ran quickly...
ssn$edges <- SSNbler::afv_edges(
  ssn$edges,
  lsn_path = path_ssn,
  infl_col = "in_km",
  segpi_col = "in_km_PI",
  afv_col = "in_km_AFV"
)

ssn$edges |> 
  ggplot() + 
  #geom_density(aes(x = in_km_AFV))
  geom_density(aes(x = in_km_PI))

#not yet run...should pick a test set of obs points
#appears to require arg 'ssn_path' that copies/creates new gpkg
#edges.gpkg written from above is 2.5mb; avoidable/aliasable for larger data?
SSNbler::ssn_assemble(
  edges = ssn$edges,
  lsn_path = path_ssn,
  
)

```


# Flow regime via NWM on 3DHP

come back to examine relative to: https://gis.dnr.wa.gov/site2/rest/services/Public_Forest_Practices/WADNR_PUBLIC_FP_Water_Type/MapServer

```{r nwm_zarr_pull3, eval=FALSE}
#revised from flow_trees_apps.qmd, which built per-year rds objects into data_common/nmw_retro, moved into /v2_1

# #https://github.com/grimbough/Rarr
# if (!require("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# ## install Rarr
# BiocManager::install("Rarr")

# remotes::install_github(repo = "grimbough/Rarr")

library(tidyverse)
library(sf)
#for COMIDs all of WA
sf_nhdp_wa <- readRDS("~/T/DFW-Team WDFW Watershed Synthesis - flow_trees_heat/sf_nhdp_wa.rds")
sf_nhdp_wa_comid <- sf_nhdp_wa |> 
    as_tibble() |> 
    dplyr::filter(COMID > 0, areasqkm < 1500) |> 
    select(comid = COMID)
rm(sf_nhdp_wa)

library(Rarr)

#zarr_overview("https://noaa-nwm-retrospective-2-1-zarr-pds.s3.amazonaws.com/chrtout.zarr/streamflow/")
#structure changes slightly between 2.1 and 3.0
zarr_overview("https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/streamflow/")


# #may be worth truncating largest and smallest?
# as_tibble(sf_nhdp_wa) |> pull(areasqkm) |> summary()
# as_tibble(sf_nhdp_wa) |> filter(between(areasqkm, 0, 1000)) |>  ggplot() + stat_ecdf(aes(areasqkm))
# as_tibble(sf_nhdp_wa) |> arrange(desc(areasqkm))
# as_tibble(sf_nhdp_wa) |> arrange(areasqkm)
# as_tibble(sf_nhdp_wa) |> count(areasqkm < 0.001)
# #several very large, largest are cross border
# sf_nhdp_wa |> filter(areasqkm > 1500) |> select(-comid) |> mapview::mapview(zcol = "areasqkm")
# #various very small, scattered throughout
# sf_nhdp_wa |> filter(areasqkm < 0.001) |> select(-comid) |> sf::st_centroid() |> mapview::mapview(zcol = "areasqkm")

# #unclear what is happening with the negative COMIDs
# #somewhat randomly in the center of the state, some have NLCD riparian attrs
# #no negatives in NWM zarr feature_id
# as_tibble(sf_nhdp_wa) |> count(COMID < 0)
# sf_nhdp_wa |> filter(COMID < 0) |> select(-comid) |> mapview::mapview(zcol = "areasqkm")
# sf_nhdp_wa_nlcd19 |> filter(COMID < 0) |> as_tibble() |> print(n=100)


#zarr 'feature_id' *values* are COMIDs, but index position in feature_id != COMID 
#this reads the vector of all available NHD+ ids
#then subsets to those in the WA catchment polygons (excluding 5 in Canada and those with negative COMIDs)
#manually finding the correct url by hovering on html view; presumably a better way
#https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/index.html#CONUS/zarr/chrtout.zarr/feature_id/

# #submitted issue:  https://github.com/grimbough/Rarr/issues/10
# #old, worked, seems likely would now fail due to string rather than list datatype
# array_path <- "https://noaa-nwm-retrospective-2-1-zarr-pds.s3.amazonaws.com/chrtout.zarr/feature_id/"
# Rarr:::read_array_metadata(array_path)
# #new, fails
# array_path <- "https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/feature_id"  
# metadata <- Rarr:::read_array_metadata(array_path)
# #decompressor <- metadata$compressor$id # decompressor == "zstd", previously/above 'lz4'
# #fails: NULL, key is '$dtype', value is '<i8' which breaks get_chunk_size()
# datatype <- metadata$datatype
# #still wrong, returns string "<i8"
# datatype <- metadata$dtype 
# #should be list from
# datatype <- Rarr:::.parse_datatype(metadata$dtype)
# #which allows declaration
# buffer_size <- Rarr:::get_chunk_size(datatype, dimensions = metadata$chunks)

## horrible debugging
# #stepping through read_chunk
# i = 1
# chunk_id = required_chunks[i, ]
# dim_separator <- ifelse(is.null(metadata$dimension_separator),
#     yes = ".", no = metadata$dimension_separator
#   )
# chunk_id <- paste0(chunk_id, collapse = dim_separator)
# #chunk_file <- paste0(zarr_array_path, chunk_id) #looks like a bad path, maybe part of problem?
# chunk_file <- paste0(zarr_array_path,"/", chunk_id) #Yes
# #in the 'else' of isnull(s3)
# parsed_url <- Rarr:::parse_s3_path(chunk_file)
# #fails/returns NULL without coercing in the '/', fetches when included
# if(Rarr:::.s3_object_exists(s3_client, parsed_url$bucket, parsed_url$object)) {
#   compressed_chunk <- s3_client$get_object(Bucket = parsed_url$bucket,
#                                            Key = parsed_url$object)$Body
# }
# decompressed_chunk <- Rarr:::.decompress_chunk(compressed_chunk, metadata)
# #this is where the "e-320" appears
# converted_chunk <- Rarr:::.format_chunk(decompressed_chunk, metadata, alt_chunk_dim)
# 
# alt_chunk_dim <- unlist(metadata$chunks) #declared in extract_elements and passed
# datatype <- metadata$datatype
# actual_chunk_size <- length(decompressed_chunk) / datatype$nbytes
# if((datatype$base_type == "py_object") || (actual_chunk_size == prod(unlist(metadata$chunks)))) {
#     chunk_dim <- unlist(metadata$chunks)
#   } else if (actual_chunk_size == prod(alt_chunk_dim)) {
#     chunk_dim <- alt_chunk_dim
#   }
# output_type <- switch(datatype$base_type,
#                       "boolean" = 0L,
#                       "int" = 1L,
#                       "uint" = 1L,
#                       "float" = 2L
# )
# converted_chunk <- .Call("type_convert_chunk", decompressed_chunk,
#                              1L, #output_type, 
#                              datatype$nbytes, 
#                              datatype$is_signed,
#                              chunk_dim,
#                              PACKAGE = "Rarr"
#     )

# #trying to actually get the data, can prob due featureID part manually as fallback but not workable for flows
# #stepping throug read_zarr_array
# zarr_array_path <- "https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/feature_id"  
# s3_client <- Rarr:::.create_s3_client(path = zarr_array_path)
# metadata <- Rarr:::read_array_metadata(zarr_array_path)
# metadata$datatype
# metadata$datatype$base_type <- "float" #corrects buffer_size but not missing slash
# index <- vector(mode = "list", length = length(metadata$shape))
# index <- Rarr:::check_index(index = index, metadata = metadata)
# required_chunks <- as.matrix(Rarr:::find_chunks_needed(metadata, index))
# 
# #now step down into read_data
# #so deconstruct Rarr:::read_data
# chunk_idx <- mapply(\(x,y) { (x-1) %/% y }, index, metadata$chunks, SIMPLIFY = FALSE)
# #now step down into chunk_selections declaration
# #which means into .extract_elements and then into read_chunk and .format_chunk


my_format_chunk <- function (decompressed_chunk, metadata, alt_chunk_dim) 
{
    datatype <- metadata$datatype
    actual_chunk_size <- length(decompressed_chunk)/datatype$nbytes
    if ((datatype$base_type == "py_object") || (actual_chunk_size == 
        prod(unlist(metadata$chunks)))) {
        chunk_dim <- unlist(metadata$chunks)
    }
    else if (actual_chunk_size == prod(alt_chunk_dim)) {
        chunk_dim <- alt_chunk_dim
    }
    else {
        stop("Decompresed data doesn't match expected chunk size.")
    }
    if (metadata$order == "C") {
        chunk_dim <- rev(chunk_dim)
    }
    if (datatype$base_type == "string") {
        converted_chunk <- .format_string(decompressed_chunk, 
            datatype)
        dim(converted_chunk[[1]]) <- chunk_dim
    }
    else if (datatype$base_type == "unicode") {
        converted_chunk <- .format_unicode(decompressed_chunk, 
            datatype)
        dim(converted_chunk[[1]]) <- chunk_dim
    }
    else if (datatype$base_type == "py_object") {
        converted_chunk <- .format_object(decompressed_chunk, 
            metadata, datatype)
        dim(converted_chunk[[1]]) <- chunk_dim
    }
    else {
        output_type <- switch(datatype$base_type, boolean = 0L, 
#            int = 1L, uint = 1L, float = 2L)
            int = 1L, uint = 1L, float = 1L) #force float to 1L, not a general fix
        converted_chunk <- .Call("type_convert_chunk", decompressed_chunk, 
            output_type, datatype$nbytes, datatype$is_signed, 
           chunk_dim, PACKAGE = "Rarr")
    }
    if (metadata$order == "C") {
        converted_chunk[[1]] <- aperm(converted_chunk[[1]])
    }
    names(converted_chunk) <- c("chunk_data", "warning")
    return(converted_chunk)
}

#changes: add Rarr::: throughout; update chunk_file declaration to include '/' 
my_read_chunk <- function (zarr_array_path, chunk_id, metadata, s3_client = NULL, 
    alt_chunk_dim = NULL) 
{
    if (missing(metadata)) {
        metadata <- Rarr:::read_array_metadata(zarr_array_path, s3_client = s3_client)
    }
    dim_separator <- ifelse(is.null(metadata$dimension_separator), 
        yes = ".", no = metadata$dimension_separator)
    chunk_id <- paste(chunk_id, collapse = dim_separator)
    datatype <- metadata$datatype
#    chunk_file <- paste0(zarr_array_path, chunk_id)
    chunk_file <- paste0(zarr_array_path,"/", chunk_id)

    if (nzchar(Sys.getenv("RARR_DEBUG"))) {
        message(chunk_file)
    }
    if (is.null(s3_client)) {
        size <- file.size(chunk_file)
        if (file.exists(chunk_file)) {
            compressed_chunk <- readBin(con = chunk_file, what = "raw", 
                n = size)
        }
        else {
            compressed_chunk <- NULL
        }
    }
    else {
        parsed_url <- Rarr:::parse_s3_path(chunk_file)
        if (Rarr:::.s3_object_exists(s3_client, parsed_url$bucket, parsed_url$object)) {
            compressed_chunk <- s3_client$get_object(Bucket = parsed_url$bucket, 
                Key = parsed_url$object)$Body
        }
        else {
            compressed_chunk <- NULL
        }
    }
    if (!is.null(compressed_chunk)) {
        decompressed_chunk <- Rarr:::.decompress_chunk(compressed_chunk, 
            metadata)
#        converted_chunk <- Rarr:::.format_chunk(decompressed_chunk, 
        converted_chunk <- my_format_chunk(decompressed_chunk, 
            metadata, alt_chunk_dim)
    }
    else {
        converted_chunk <- list(chunk_data = array(metadata$fill_value, 
            dim = unlist(metadata$chunks)), warning = 0L)
    }
    return(converted_chunk)
}

#need to declare here as not even ':::' available
#also changing to:   chunk <- my_read_chunk(
my_extract_elements <- function(i, metadata, index, required_chunks, zarr_array_path, s3_client, chunk_idx) {
  ## find elements to select from the chunk and what in the output we replace
  index_in_result <- index_in_chunk <- list()
  alt_chunk_dim <- unlist(metadata$chunks)
  
  for (j in seq_len(ncol(required_chunks))) {
    index_in_result[[j]] <- which(chunk_idx[[j]] == required_chunks[i, j])
    ## are we requesting values outside the array due to overhanging chunks?
    outside_extent <- index_in_result[[j]] > metadata$shape[[j]]
    if (any(outside_extent))
      index_in_result[[j]] <- index_in_result[[j]][-outside_extent]
    if (any(index_in_result[[j]] == metadata$shape[[j]])) 
      alt_chunk_dim[j] <- length(index_in_result[[j]])
    
    index_in_chunk[[j]] <- ((index[[j]][index_in_result[[j]]] - 1) %% metadata$chunks[[j]]) + 1
  }
  
  ## read this chunk
  #chunk <- Rarr:::read_chunk(
  chunk <- my_read_chunk(
    zarr_array_path,
    chunk_id = required_chunks[i, ],
    metadata = metadata,
    s3_client = s3_client,
    alt_chunk_dim = alt_chunk_dim
  )

  warn <- chunk$warning[1]
  chunk_data <- chunk$chunk_data
  
  ## extract the required elements from the chunk
  selection <- R.utils::extract(chunk_data, indices = index_in_chunk, drop = FALSE)
  
  return(list(selection, index_in_result, warning = warn))
}

#change to FUN = .extract_elements
my_read_data <- function(required_chunks, zarr_array_path, s3_client, 
                      index, metadata) {
  warn <- 0L
  
  ## determine which chunk each of the requests indices belongs to
  chunk_idx <- mapply(\(x,y) { (x-1) %/% y }, index, metadata$chunks, SIMPLIFY = FALSE)
  
  ## hopefully we can eventually do this in parallel
  chunk_selections <- lapply(seq_len(nrow(required_chunks)), 
                             #FUN = .extract_elements,
                             FUN = my_extract_elements,
                             metadata = metadata, index = index,
                             required_chunks = required_chunks,
                             zarr_array_path = zarr_array_path,
                             s3_client = s3_client,
                             chunk_idx = chunk_idx)
  
  ## predefine our array to be populated from the read chunks
  output <- array(metadata$fill_value, dim = vapply(index, length, integer(1)))

  ## proceed in serial and update the output with each chunk selection in turn
  for (i in seq_along(chunk_selections)) {
    index_in_result <- chunk_selections[[i]][[2]]
    cmd <- Rarr:::.create_replace_call(x_name = "output", idx_name = "index_in_result",
                                idx_length = length(index_in_result), 
                                y_name = "chunk_selections[[i]][[1]]")
    eval(parse(text = cmd))
    warn <- max(warn, chunk_selections[[i]]$warning[1])
  }
  return(list(output = output, warn = warn))
}

my_read_zarr_array <- function (zarr_array_path, index, s3_client) 
{
    zarr_array_path <- Rarr:::.normalize_array_path(zarr_array_path)
    if (missing(s3_client)) 
        s3_client <- Rarr:::.create_s3_client(path = zarr_array_path)
    metadata <- Rarr:::read_array_metadata(zarr_array_path, s3_client = s3_client)
    #yikes
    metadata$datatype$base_type <- "float"
    
    if (missing(index)) {
        index <- vector(mode = "list", length = length(metadata$shape))
    }
    index <- Rarr:::check_index(index = index, metadata = metadata)
    required_chunks <- as.matrix(Rarr:::find_chunks_needed(metadata, index))
    res <- my_read_data(required_chunks, zarr_array_path, s3_client, index, metadata)
    if (isTRUE(res$warn > 0)) {
        warning("Integer overflow detected in at least one chunk.\n", 
            "Overflowing values have been replaced with NA", 
            call. = FALSE)
    }
    return(res$output)
}

#res <- Rarr:::read_data(required_chunks, zarr_array_path, s3_client, index, metadata)
#my_read_data(required_chunks, zarr_array_path, s3_client, index, metadata)$output

#Now runs...to 56K tibble
feat_id <- tibble(
  # #v2.1: comid = read_zarr_array("https://noaa-nwm-retrospective-2-1-zarr-pds.s3.amazonaws.com/chrtout.zarr/feature_id/")
  comid = my_read_zarr_array("https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/feature_id")
) |> 
  tibble::rowid_to_column(var = "i") |> 
  inner_join(sf_nhdp_wa_comid, by = "comid")

# {
#     "_ARRAY_DIMENSIONS": [
#         "time",
#         "feature_id"
#     ],
#     "add_offset": 0.0,
#     "coordinates": "latitude longitude",
#     "grid_mapping": "crs",
#     "long_name": "River Flow",
#     "missing_value": -999900,
#     "scale_factor": 0.009999999776482582,
#     "units": "m3 s-1"
# }
# Shape: 385704 x 2776734
# Chunk Shape: 672 x 30000
# No. of Chunks: 53382 (574 x 93)
# Data Type: int32

# q <- my_read_zarr_array(
#   "https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/streamflow/"
#   ,
#   index = list(1, feat_id$i)
# )


#for now hard coding expectation of passing year, month and hour(s)
#could add logic for other NULL handling to build time_dim_ind more flexibly
#defaults for single hour per day through full month (13:00 Aug 2020)
#may be faster to pull, for example, Apr-Sep rather than looping single months
#but then have [too] big objects and have to load and subset?
dir_data_common <- "~/T/DFW-Team WDFW Watershed Synthesis - data_common"

pull_nwm_zarr <- function(
    #s3_address = "https://noaa-nwm-retrospective-2-1-zarr-pds.s3.amazonaws.com/chrtout.zarr/streamflow/",
    s3_address =   "https://noaa-nwm-retrospective-3-0-pds.s3.amazonaws.com/CONUS/zarr/chrtout.zarr/streamflow/",
    y = 2022,
    m = 8,
    h = 13,
    feat = feat_id,
    dir_write = file.path(dir_data_common, "nwm_retro/v3_0")
) {
  sclfctr <- 0.009999999776482582
  cms2cfs <- 35.31468492
  
  #vector of all time available, v2.1 length was 367439
  time_dim <- seq(as.POSIXct("1979-02-01 01:00:00", tz = 'UTC'), length.out = 385704, by = "1 hour")
  #indices of desired subset
  time_i <- which(
    #year(time_dim) %in% 2023 & month(time_dim) %in% 8 & hour(time_dim) %in% 13
    year(time_dim) %in% y &
      month(time_dim) %in% m &
      hour(time_dim) %in% h
    )
  
  print(Sys.time())
  z <- my_read_zarr_array(s3_address, index = list(time_i, feat$i)) |> 
    as.data.frame() |> 
    set_names(feat$comid) |> 
    as_tibble() |> 
    mutate(
      #correct 'scale factor' brings to cms, then to cfs
      across(everything(), ~.* sclfctr * cms2cfs),
      across(everything(), ~if_else(.<0, NA_real_, .)),
      i = time_dim[time_i], #col of POSIXct
      m = month(i),
      yday = yday(i)
    )
  print(Sys.time())
  
  if(!is.null(dir_write)){
    if(identical(m, 1:12)) {
      f = file.path(dir_write, paste0("nwm_wa_",y, ".rds"))
    } else {
      f = file.path(dir_write, paste0("nwm_wa_",y,"_",paste0(str_pad(m,width = 2, pad = "0"),collapse = ""), ".rds"))
    }
    print(paste("data written to ", f))
    saveRDS(z, f)
  } else {
    return(z)
  }
}

# #365 days at 13:00 
# q <- pull_nwm_zarr(dir_write = NULL, y = 2022, m = 12)
# dim(as.matrix(q)); tail(names(q))
# q2 <- q |> select(-c(i:yday)) |> as.matrix()
# dim(q2)
# apply(q2,1, max, na.rm = T) #impossibly high values...>200k cfs
# rowMeans(q2)
# apply(q2,1, \(x) sum(is.na(x)))
# apply(q2,1, \(x) sum(x<0, na.rm = T)) #fixed, could drop these comids?
# #colnames(q2)[q2[1,]<0]
# colnames(q2)[q2[1,]<0]
# q2 <- q2[,!(is.na(q2[1,]))]
# rowMeans(q2)
# apply(q2,1, min)
# #okay looks like just a few weird/bad spots but the rescaling is correct
# plot.ecdf(q2[1,]) ; abline(v=c(10000,50000), col = 2:3)

# #default single hour (13:00) per day
# about 1m30s per y-m
# some change somewhere seems to have made less memory stable? 
# throws stubborn mem limit error after about 10 y-m
# gc() and googling did not show a better quick fix than just restarting session
# could script a meta-loop that spawns & provisions new session(s)
# either with Rscript or callr::r_bg, https://callr.r-lib.org/ 
# but this feels provisional and easy enough to babysit while doing other things
x <- expand_grid(y = 2021:2000, m = 1:12)
(x <- x |> filter(!(y == 2021 & m < 12)))
# for(i in seq_len(nrow(x))){
#   print(x[i,"y"]); print(x[i,"m"])
#   pull_nwm_zarr(y = x[i,"y"], m = x[i,"m"])
# }

#got annoyed and wrote something that starts fresh each y-m and seems more memory stable?
file.exists("~/T/DFW-Team WDFW Watershed Synthesis - General/hydrography/testing/nwm3_hacky_ym.R")

for(i in seq_len(nrow(x))){
  system2(
    command = "Rscript",
    args = c(
      #tilde expansion does NOT work
      shQuote("/Users/auerbdaa/T/DFW-Team WDFW Watershed Synthesis - General/hydrography/testing/nwm3_hacky_ym.R"),
      x[i,"y"], x[i,"m"]
    )
  )
  }

```

```{r nwm_on_3dhp}

#NWM 2022 testing values
q_3dhp <- q |> 
  select(m, yday, which(colnames(q) %in% as.character(unique(lines_3dhp$comid)))) |> 
  mutate(across(-c(m, yday), ~if_else(.<0 , NA_real_, .)))

q_3dhp_daily_summary <- q_3dhp |> 
  pivot_longer(cols = -c(m, yday), names_to = "comid", values_to = "cfs") |> 
  summarise(
    across(
      cfs, list(
        min = ~min(.), med = ~median(.), max = ~max(.),
        cv = ~sd(.) / mean(.)
        )
      ),
    .by = "comid"
    )

#add back to spatial objects
sf_nhdp_wa_aoi <- sf_nhdp_wa_aoi |> left_join(q_3dhp_daily_summary, by = "comid")
lines_3dhp <- lines_3dhp |> left_join(q_3dhp_daily_summary, by = "comid")

map(
  names(q_3dhp_daily_summary)[-1]
  ,
  ~sf_nhdp_wa_aoi |> 
    drop_na(cfs_min) |> 
    ggplot() + geom_sf(aes(color = eval(as.name(.x)), fill = eval(as.name(.x)))) +
    wacolors::scale_color_wa_c("diablo", reverse = T, aesthetics = c("color","fill")) +
    guides(color = guide_colorbar(.x), fill = guide_none())
  ) |> 
  wrap_plots(ncol = 2) +
  plot_annotation(subtitle = "NWM 3.0 2022 per med res NHDplus catchments")


lines_3dhp |> 
  ggplot() +
  geom_sf(aes(color = cfs_cv)) + 
  wacolors::scale_color_wa_c("stuart", midpoint = 1)


```