diff --git a/workspace/R/NHD_navigate.R b/workspace/R/NHD_navigate.R
index d0de0145e21cfe668ce8e02951b31f77249287a6..9deee9b020b0a075d48be9cf1241b82595d715cf 100644
--- a/workspace/R/NHD_navigate.R
+++ b/workspace/R/NHD_navigate.R
@@ -409,8 +409,10 @@ movePOI_NA_DA <- function(poi_fix, nhdDF){
# 2/3 - correpsonding segments (both raw and dissolved)
POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
+ nexus <- filter(POIs, nexus == TRUE)
+
POIs <- POIs %>% mutate_if(is.numeric, function(x) ifelse(is.infinite(x), NA, x)) %>%
- filter(nexus != TRUE)
+ filter(nexus != 1)
# downstream POIs and their drainage area ratios and incremental drainage areas
segs_down <- Seg %>%
@@ -530,7 +532,8 @@ POI_move_down<-function(out_gpkg, POIs, Seg, Seg2, exp, DA_diff){
# Bring new POI data over to new data
fin_POIs <- filter(POIs, !COMID %in% c(merged_POIs$oldPOI, merged_POIs$COMID)) %>%
- rbind(merged_POIs %>% select(-c(DirDA, oldPOI)))
+ rbind(merged_POIs %>% select(-c(DirDA, oldPOI))) %>%
+ rbind(nexus)
#inner_join(st_drop_geometry(merged_POIs) %>% select(COMID, Type_Gages_A), by = "COMID") %>%
#mutate(Type_Gages = ifelse(!is.na(Type_Gages_A), Type_Gages_A, Type_Gages)) %>% select(-Type_Gages_A)
@@ -832,6 +835,7 @@ HR_Area_coms <- function(nhd, WBs, data_paths, crs){
# Download NHD HR HUC04 if we dont' have it, other wise load and
# Bind NHDHR waterbodies into one layer
hr_wb <- do.call("rbind", lapply(unique(huc04), function(x){
+ print(x)
if(!file.exists(file.path(data_paths$nhdplus_dir, vpu,
paste0("NHDPLUS_H_", x, "_HU4_GDB.gdb")))){
@@ -848,7 +852,8 @@ HR_Area_coms <- function(nhd, WBs, data_paths, crs){
st_zm(.) %>%
st_as_sf() %>%
mutate(source = "NHDHR") %>%
- st_transform(crs)
+ st_transform(crs) #%>%
+ #st_compatibalize(wb)
# Bind or create new object
if(exists("wb")){
@@ -897,12 +902,14 @@ WB_event <- function(WBs, nhd_wb, type){
outlet_fl <- filter(nhd_wb, COMID %in% filter(WBs_table, outlet == "outlet")$comid)
# Get the downstream flowline for continuity
- ds_fl <- filter(nhd_wb, DnHydroseq %in% outlet_fl$Hydroseq,
+ ds_fl <- filter(nhd_wb, Hydroseq %in% outlet_fl$DnHydroseq,
LevelPathI %in% outlet_fl$LevelPathI) %>%
rbind(outlet_fl) %>%
- group_by(WBAREACOMI) %>%
+ arrange(desc(Hydroseq)) %>%
+ group_by(LevelPathI) %>%
# union together
- summarize(do_union = T)
+ summarize(do_union = T) %>%
+ st_cast("LINESTRING")
WBs_area <- filter(WBs_layer, source == "NHDv2Area")
WBs_HR <- filter(WBs_layer, source == "NHDHR")
@@ -911,13 +918,29 @@ WB_event <- function(WBs, nhd_wb, type){
if (nrow(WBs_area) > 0){
# Intersect unioned FL with waterbody and get intersecting point
- outlet_pnts <- sf::st_intersection(ds_fl, WBs_area) %>%
- # Cast to point
+ outlet_pnts_int <- sf::st_intersection(ds_fl, WBs_area) #%>%
+ # # Cast to point
+ # st_cast("POINT")
+ # group_by(LevelPathI) %>%
+ # # Get furthest downstream point of intersection
+ # filter(row_number() == max(row_number(), na.rm = T)) %>%
+ # ungroup()
+
+ outlet_point <- outlet_pnts_int[st_geometry_type(outlet_pnts_int) %in% c("POINT"),]
+ if(nrow(outlet_point) > 0){
+ outlet_fl <- outlet_pnts_int[st_geometry_type(outlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
+ st_cast("LINESTRING") %>%
+ st_cast("POINT")
+ outlet_pnts_int <- rbind(outlet_point, outlet_fl)
+ } else {
+ outlet_pnts_int <- st_cast(outlet_pnts_int, "MULTILINESTRING")
+ }
+
+ outlet_pnts <- outlet_pnts_int %>%
+ group_by(LevelPathI) %>%
st_cast("POINT") %>%
- group_by(WBAREACOMI) %>%
- # Get furthest downstream point of intersection
- filter(row_number() == max(row_number(), na.rm = T)) %>%
- ungroup()
+ mutate(id = row_number()) %>%# Was comid
+ filter(id == max(id))
# Derive event from point to use for splitting within hyRefactor
wb_events <- get_flowline_index(nhd_wb, outlet_pnts) %>%
@@ -1002,29 +1025,50 @@ WB_event <- function(WBs, nhd_wb, type){
inlet_FL <- nhd_wb[lengths(st_intersects(start_pts, WBs_layer)) == 0,] %>%
rbind(filter(nhd_wb, Hydroseq %in% .$DnHydroseq,
- LevelPathI %in% .$LevelPathI))
+ LevelPathI %in% .$LevelPathI)) %>%
+ arrange(desc(Hydroseq)) %>%
+ unique()
inlet_ls <- inlet_FL %>%
group_by(LevelPathI) %>%
st_cast("POINT") %>%
summarize(do_union = F) %>%
- st_cast("LINESTRING")
+ st_cast("LINESTRING") %>%
+ ungroup()
- inlet_pnts <- sf::st_intersection(inlet_ls, WBs_layer)
-
- # the intersection can return linestring features.
- if(sf::st_geometry_type(inlet_pnts, by_geometry = FALSE) == "LINESTRING") {
- inlet_pnts <- inlet_pnts %>%
- st_cast("POINT") %>%
- st_as_sf() %>%
- group_by(COMID) %>%
- filter(row_number() == 1) %>%
- ungroup()
- } else {
- inlet_pnts <- inlet_pnts %>%
- st_cast("POINT") %>%
- st_as_sf()
+ #inlet_pnts <- sf::st_intersection(inlet_ls, st_cast(WBs_layer, "MULTILINESTRING", group_or_split = FALSE)) #%>%
+
+ inlet_pnts_int <- sf::st_intersection(inlet_ls, WBs_layer) #%>%
+ #st_cast("LINESTRING")
+ int_point <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("POINT"),]
+ if(nrow(int_point) > 0){
+ inlet_fl <- inlet_pnts_int[st_geometry_type(inlet_pnts_int) %in% c("LINESTRING", "MULTILINESTRING"),] %>%
+ st_cast("LINESTRING")
+ inlet_pnts_int <- rbind(int_point, inlet_fl)
}
+
+ inlet_pnts <- inlet_pnts_int %>%
+ group_by(LevelPathI) %>%
+ st_cast("POINT") %>%
+ mutate(id = row_number()) %>%# Was comid
+ filter(id == min(id))
+
+ # # the intersection can return linestring features.
+ # if(sf::st_geometry_type(inlet_pnts, by_geometry = FALSE) %in% c("LINESTRING", "MULTILINESTRING")) {
+ # inlet_pnts2 <- inlet_pnts %>%
+ # #st_cast("LINESTRING") %>%
+ # st_cast("POINT") %>%
+ # group_by(COMID) %>%
+ # mutate(id = row_number()) %>%# Was comid
+ # filter(id == 1)
+ # #filter(row_number() == 1) %>%
+ # #filter(row_number() == min(row_number())) %>%
+ # ungroup()
+ # } else {
+ # inlet_pnts <- inlet_pnts %>%
+ # st_cast("POINT") %>%
+ # st_as_sf()
+ # }
wb_events <- get_flowline_index(nhd_wb, inlet_pnts) %>%
inner_join(select(st_drop_geometry(nhd_wb), COMID, WBAREACOMI, LevelPathI), by = "COMID") %>%
@@ -1062,7 +1106,7 @@ gage_POI_creation <- function(tmp_POIs, gages_info, nhd, combine_meters, reach_m
# Reset gage flag if even is created
gage_POIs_nonevent <- filter(gage_POIs, !Type_Gages %in% events$Type_Gages) %>%
- addType(., tmp_POIs, "Gages")
+ addType(., tmp_POIs, "Gages", nexus = FALSE, bind = TRUE)
tmp_POIs <- data.table::rbindlist(list(gage_POIs_nonevent,
select(events, COMID, Type_Gages, nexus)), fill = TRUE) %>%
@@ -1206,12 +1250,16 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
# Drive the inlet events
nhd_inlet <- filter(nhd_WBs, !COMID %in% wb_outlet_events$COMID)
+
+ if(unique(nhd$VPUID) == "12"){
+ nhd_inlet <- rbind(nhd_inlet, filter(nhd_WBs, COMID == 1304709))
+ }
+
if (nrow(nhd_inlet) > 0){
# Get inlet events and bind with outlets
wb_inlet_events <- WB_event(WBs_VPU_areaHR, nhd_inlet, "inlet") %>%
- st_compatibalize(wbin_POIs) %>%
mutate(nexus = TRUE) %>%
- inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq, ToMeas), by = "COMID")
+ inner_join(select(st_drop_geometry(nhd), COMID, Hydroseq, ToMeas), by = "COMID")
# Determine which events are close enough to the end of an upstream flowline to just use that geometry for the POI
wbin_fl_upstream <- filter(nhd, DnHydroseq %in% filter(nhd, COMID %in% wb_inlet_events$COMID)$Hydroseq) %>%
@@ -1227,8 +1275,8 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
mutate(dsCOMID = COMID, COMID = usCOMID)
if(nrow(wb_inlet_POIs) > 0) {
- wbin_POIs <- bind_rows(POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI),
- nhd, "WBIn"), wbin_POIs)
+ wbin_POIs <- bind_rows(wbin_POIs, POI_creation(select(st_drop_geometry(wb_inlet_POIs), dsCOMID, Type_WBIn = WBAREACOMI),
+ nhd, "WBIn"))
wb_inlet_events <- filter(wb_inlet_events, !COMID %in% wb_inlet_POIs$dsCOMID)
}
@@ -1237,6 +1285,7 @@ wbin_POIcreation <- function(nhd, WBs_VPU, data_paths, crs,
if(nrow(wb_inlet_events) > 0){
wbin_POIs <- addType(wbin_POIs, tmp_POIs, "WBIn")
+ wb_inlet_events <- st_compatibalize(wb_inlet_events, wbin_POIs)
wbin_POIs_fin <- data.table::rbindlist(list(wbin_POIs,
select(wb_inlet_events, COMID,
Type_WBIn = WBAREACOMI, nexus)),