diff --git a/R/dam_locations.R b/R/dam_locations.R index 5c8abf5..f4c4a35 100644 --- a/R/dam_locations.R +++ b/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")) }