add precise locations and fix up mainstem additions

R/dam_locations.R

@@ -158,6 +158,69 @@ get_dam_locations <- function(dams, nid) {
            feature_data_source)
 }
 
+utility_link_dams <- function(to_link, dams, nhdpv2_fline, vaa, search_radius_m, max_matches_in_radius) {
+  new_hl <- nhdplusTools::get_flowline_index(nhdpv2_fline, 
+                                             to_link, 
+                                             search_radius = units::set_units(
+                                               search_radius_m, "m"),
+                                             max_matches = max_matches_in_radius)
+
+  sf::st_drop_geometry(select(to_link, provider_id)) %>%
+    mutate(id = seq_len(nrow(.))) %>%
+    left_join(new_hl, by = "id") %>%
+    left_join(select(sf::st_drop_geometry(dams), 
+                     provider_id, drainage_area_sqkm), 
+              by = "provider_id") %>%
+    left_join(select(sf::st_drop_geometry(nhdpv2_fline),
+                     COMID, drainage_area_sqkm_nhdpv2 = TotDASqKM), 
+              by = "COMID") %>%
+    left_join(select(vaa, COMID = comid, hydroseq), by = "COMID") %>%
+    mutate(da_diff = abs(drainage_area_sqkm - drainage_area_sqkm_nhdpv2) / drainage_area_sqkm)
+}
+
+get_dam_network_location <- function(dams, update_index, nhdpv2_fline, vaa, 
+                                     search_radius_m, max_matches_in_radius, index_type) {
+  linked_dams_base <- select(dams[update_index, ], provider_id, nhdpv2_COMID)
+
+  linked_dams <- linked_dams_base %>%
+    left_join(select(sf::st_drop_geometry(nhdpv2_fline), COMID, REACHCODE, FromMeas), 
+              by = c("nhdpv2_COMID" = "COMID"))
+
+  # get precise network location
+
+  fline_sites <- dplyr::filter(nhdpv2_fline, COMID %in% linked_dams$nhdpv2_COMID) |>
+    sf::st_transform(sf::st_crs(linked_dams)) |>
+    sf::st_cast("LINESTRING")
+
+  linked_dams <- utility_link_dams(linked_dams, dams, fline_sites, vaa, 
+                                   search_radius_m, max_matches_in_radius)
+
+  linked_dams_dedup <- linked_dams %>%
+    left_join(select(sf::st_drop_geometry(dams), provider_id, nhdpv2_COMID), by = "provider_id") %>%
+    group_by(provider_id) %>%
+    filter(COMID == nhdpv2_COMID) %>%
+    ungroup()
+
+  if(!any(duplicated(linked_dams_dedup$provider_id))) {
+    linked_dams_base <- linked_dams_base %>%
+      left_join(select(linked_dams_dedup, provider_id, REACHCODE, REACH_meas, drainage_area_sqkm_nhdpv2),
+                by = "provider_id")
+  } else {
+    stop()
+  }
+
+  if(!nrow(linked_dams_base) == sum(update_index)) stop()
+
+  dams$nhdpv2_REACHCODE[update_index] <- 
+    linked_dams_base$REACHCODE
+  dams$nhdpv2_REACH_measure[update_index] <- 
+    linked_dams_base$REACH_meas
+
+  dams$index_type[update_index] <- index_type 
+
+  dams
+}
+
 #' get dam hydrolocations
 #' @param dams sf data.frame table of dams with best estimate of comid and drainage area
 #' @param nhdpv2_fline sf data.frame of flowlines to index to
@@ -208,17 +271,23 @@ get_dam_hydrolocations <- function(dams, nhdpv2_fline, vaa,
 
   if(any(update_index)) {
 
-    linked_dams <- select(dams[update_index, ], provider_id, nhdpv2_COMID) %>%
-      left_join(select(sf::st_drop_geometry(nhdpv2_fline), COMID, REACHCODE, FromMeas), 
-                by = c("nhdpv2_COMID" = "COMID"))
-
-    #TODO: could figure out a precise location here -- but not that critical.
-    dams$nhdpv2_REACHCODE[update_index] <- 
-      linked_dams$REACHCODE
-    dams$nhdpv2_REACH_measure[update_index] <- 
-      linked_dams$FromMeas
+    dams <- get_dam_network_location(dams, update_index, nhdpv2_fline, vaa, 
+                                     search_radius_m, max_matches_in_radius, 
+                                     "nawqa_on_network_da_match")
+  }
 
-    dams$index_type[update_index] <- "nawqa_on_network_da_match"  
+  # these are where we have a small catchment worth indexing to despite a da mismatch
+  update_index <- is.na(dams$nhdpv2_REACH_measure) & # no reachcode measure yet
+    !is.na(dams$nhdpv2_COMID) & # has a comid
+    (!is.na(dams$norm_diffda) & # has a normalized area diff
+       abs(dams$norm_diffda) > da_diff_thresh & # diff is larger than thresh
+       dams$drainage_area_sqkm_nhdpv2 < 10) # but this is near a headwater
+
+  if(any(update_index)) {
+
+    dams <- get_dam_network_location(dams, update_index, nhdpv2_fline, vaa, 
+                                     search_radius_m, max_matches_in_radius, 
+                                     "nawqa_<10sqkm_da_mismatch")
   }
 
   # these are where we don't have a network drainage area match but NAWQA assigned a COMID
@@ -236,23 +305,8 @@ get_dam_hydrolocations <- function(dams, nhdpv2_fline, vaa,
 
   no_location <- sf::st_transform(no_location, 5070)
 
-  new_hl <- nhdplusTools::get_flowline_index(nhdpv2_fline, 
-                                             no_location, 
-                                             search_radius = units::set_units(
-                                               search_radius_m, "m"),
-                                             max_matches = max_matches_in_radius)
-
-  linked_dams <- sf::st_drop_geometry(select(no_location, provider_id)) %>%
-    mutate(id = seq_len(nrow(.))) %>%
-    left_join(new_hl, by = "id") %>%
-    left_join(select(sf::st_drop_geometry(dams), 
-                     provider_id, drainage_area_sqkm), 
-              by = "provider_id") %>%
-    left_join(select(sf::st_drop_geometry(nhdpv2_fline),
-                     COMID, drainage_area_sqkm_nhdpv2 = TotDASqKM), 
-              by = "COMID") %>%
-    left_join(select(vaa, COMID = comid, hydroseq), by = "COMID") %>%
-    mutate(da_diff = abs(drainage_area_sqkm - drainage_area_sqkm_nhdpv2) / drainage_area_sqkm) %>%
+  linked_dams <- utility_link_dams(no_location, dams, nhdpv2_fline, vaa, 
+                                   search_radius_m, max_matches_in_radius) %>%
     filter(is.na(da_diff) | da_diff < da_diff_thresh)
 
   linked_dams_dedup <- bind_rows(
@@ -291,16 +345,25 @@ get_dam_hydrolocations <- function(dams, nhdpv2_fline, vaa,
 }
 
 add_mainstems <- function(dam_hydrologic_locations, mainstems, vaa) {
-  mainstems <- mainstems[,c("id", "uri"), drop = TRUE]
-  mainstems$id <- as.integer(mainstems$id)
-  vaa <- right_join(vaa, mainstems, by = c("levelpathi" = "id"))
-
-  vaa <- vaa[,c("comid", "uri")]
-
-  names(vaa) <- c("comid", "mainstem_uri")
-
-  vaa <- vaa[!is.na(vaa$comid),]
+  mainstems <- mainstems[,c("head_nhdpv2_COMID", "uri"), drop = TRUE]
+
+  mainstems$head_nhdpv2_COMID <- as.integer(gsub("https://geoconnex.us/nhdplusv2/comid/", "", 
+                                                 mainstems$head_nhdpv2_COMID))
+
+  # stolen from ref_gages
+  mainstem_lookup <- group_by(vaa, levelpathi) |>
+    filter(hydroseq == max(hydroseq)) |>
+    ungroup() |>
+    select(head_nhdpv2_COMID = comid, levelpathi) |>
+    distinct() |>
+    left_join(mainstems, by = "head_nhdpv2_COMID") |>
+    filter(!is.na(uri)) |>
+    select(-head_nhdpv2_COMID) |>
+    right_join(select(vaa, comid, levelpathi), 
+               by = "levelpathi") |>
+    select(-levelpathi, comid, mainstem_uri = uri)
+
+  dplyr::left_join(dam_hydrologic_locations, mainstem_lookup, 
+                   by = c("nhdpv2_COMID" = "comid"))
 
-  left_join(dam_hydrologic_locations, vaa, 
-            by = c("nhdpv2_COMID" = "comid"))
 }