diff --git a/R/get_drainage_area_estimates.R b/R/get_drainage_area_estimates.R index a6f4ad19..afc93647 100644 --- a/R/get_drainage_area_estimates.R +++ b/R/get_drainage_area_estimates.R @@ -23,7 +23,15 @@ #' the \code{ncontrb_a} (non-contributing acres) attribute on HUC12 features. #' #' @param start list with \code{featureSource} and \code{featureID} compatible -#' with \code{\link{get_nldi_feature}}. +#' with \code{\link{get_nldi_feature}}. May also include optional +#' \code{reachcode} and \code{measure} fields; when present these override the +#' NLDI-supplied values used to position the outlet-catchment split point, +#' without affecting upstream network navigation. NLDI can be bypassed +#' entirely by passing only \code{reachcode} and \code{measure} (no +#' \code{featureSource}/\code{featureID}); in that case the start COMID is +#' resolved from the reachcode via the VAA when \code{local_navigation = TRUE}, +#' otherwise via the NHD OGC API, and the start point geometry is computed +#' from the flowline at the supplied measure. #' @param catchments logical. If TRUE, fetch and return NHDPlusV2 catchment #' polygons for the full upstream network. Default FALSE. #' @param nhdplushr logical. If TRUE (the default), compute a drainage area @@ -357,10 +365,43 @@ resolve_start_feature <- function(start, vaa = NULL, message("Found ", length(outlet_comids), " outlets for waterbody") start_feature <- wb + } else if(is.list(start) && !is.null(start$reachcode) && + !is.null(start$measure) && + is.null(start$featureSource) && is.null(start$featureID)) { + message("Resolving start COMID from reachcode/measure ", + "(bypassing NLDI)...") + measure <- as.numeric(start$measure) + comid <- resolve_comid_from_reachcode(start$reachcode, measure, vaa) + message(" Resolved COMID = ", comid) + + fl <- get_nhdplus(comid = comid, realization = "flowline") + if(is.null(fl) || nrow(fl) == 0) + stop("Could not fetch flowline geometry for COMID ", comid) + indexes <- data.frame( + COMID = comid, + REACHCODE = start$reachcode, + REACHCODE_measure = measure, + offset = 0 + ) + start_point <- get_hydro_location(indexes, dplyr::rename(fl, id = "comid")) + + start_feature <- sf::st_sf( + comid = comid, + reachcode = start$reachcode, + measure = measure, + geometry = sf::st_sfc(start_point, crs = sf::st_crs(fl)) + ) + outlet_comids <- comid + } else { message("Resolving start feature via NLDI...") - start_feature <- get_nldi_feature(start) + start_feature <- get_nldi_feature(start[c("featureSource", "featureID")]) outlet_comids <- start_feature$comid + + if(!is.null(start$reachcode)) + start_feature$reachcode <- start$reachcode + if(!is.null(start$measure)) + start_feature$measure <- as.numeric(start$measure) } if(length(outlet_comids) == 0 || all(is.na(outlet_comids))) @@ -370,6 +411,55 @@ resolve_start_feature <- function(start, vaa = NULL, list(start_feature = start_feature, outlet_comids = outlet_comids) } +#' Resolve a COMID from a reachcode + measure +#' +#' Used by \code{resolve_start_feature} when the caller bypasses NLDI by passing +#' a reachcode/measure directly in \code{start}. A reachcode can contain +#' multiple flowlines; the one whose \code{[frommeas, tomeas]} interval covers +#' \code{measure} is selected. If \code{measure} is outside every interval, the +#' flowline whose interval boundary is closest is chosen. +#' +#' @param reachcode character. NHDPlusV2 reachcode (14-digit string). +#' @param measure numeric. Reach measure 0-100. +#' @param vaa data.frame or NULL. NHDPlusV2 VAA table with \code{comid}, +#' \code{reachcode}, \code{frommeas}, \code{tomeas} columns. When NULL, the +#' USGS NHD OGC API is queried by reachcode. +#' @return integer COMID. +#' @noRd +resolve_comid_from_reachcode <- function(reachcode, measure, vaa = NULL) { + measure <- as.numeric(measure) + + if(!is.null(vaa)) { + rows <- vaa[vaa$reachcode == reachcode, + c("comid", "frommeas", "tomeas")] + if(nrow(rows) == 0) + stop("No flowlines found in VAA for reachcode ", reachcode) + } else { + base <- get("usgs_water_root", envir = nhdplusTools_env) + out <- get_oafeat( + base = base, AOI = NULL, type = "nhdflowline_network", + filter = paste0("reachcode%20=%20%27", reachcode, "%27"), + properties = c("comid", "frommeas", "tomeas"), + skip_geometry = TRUE, status = FALSE + ) + if(is.null(out) || nrow(out) == 0) + stop("No flowlines found via OGC API for reachcode ", reachcode) + rows <- data.frame( + comid = as.integer(out$comid), + frommeas = as.numeric(out$frommeas), + tomeas = as.numeric(out$tomeas) + ) + } + + in_bounds <- measure >= rows$frommeas & measure <= rows$tomeas + if(any(in_bounds)) + return(as.integer(rows$comid[which(in_bounds)[1]])) + + dist <- pmin(abs(measure - rows$frommeas), + abs(measure - rows$tomeas)) + as.integer(rows$comid[which.min(dist)]) +} + #' Find the HUC12 containing a point #' #' Used to detect closed-basin / endorheic point starts whose containing HUC12 diff --git a/man/get_drainage_area_estimates.Rd b/man/get_drainage_area_estimates.Rd index 57b7768a..d88b4734 100644 --- a/man/get_drainage_area_estimates.Rd +++ b/man/get_drainage_area_estimates.Rd @@ -19,7 +19,15 @@ get_drainage_area_estimates( } \arguments{ \item{start}{list with \code{featureSource} and \code{featureID} compatible -with \code{\link{get_nldi_feature}}.} +with \code{\link{get_nldi_feature}}. May also include optional +\code{reachcode} and \code{measure} fields; when present these override the +NLDI-supplied values used to position the outlet-catchment split point, +without affecting upstream network navigation. NLDI can be bypassed +entirely by passing only \code{reachcode} and \code{measure} (no +\code{featureSource}/\code{featureID}); in that case the start COMID is +resolved from the reachcode via the VAA when \code{local_navigation = TRUE}, +otherwise via the NHD OGC API, and the start point geometry is computed +from the flowline at the supplied measure.} \item{catchments}{logical. If TRUE, fetch and return NHDPlusV2 catchment polygons for the full upstream network. Default FALSE.} diff --git a/tests/testthat/test_get_drainage_area_estimates.R b/tests/testthat/test_get_drainage_area_estimates.R index 9222fab2..2c029d12 100644 --- a/tests/testthat/test_get_drainage_area_estimates.R +++ b/tests/testthat/test_get_drainage_area_estimates.R @@ -316,6 +316,86 @@ test_that("get_drainage_area_estimates Black Earth Creek smoke test", { expect_s3_class(result$split_catchment, "sf") }) +test_that("resolve_comid_from_reachcode picks flowline by measure interval", { + vaa <- data.frame( + comid = c(101L, 102L, 103L, 999L), + reachcode = c("14070003001234", "14070003001234", + "14070003001234", "14070003009999"), + frommeas = c(0, 50, 80, 0), + tomeas = c(50, 80, 100, 100) + ) + + # measure in second interval + expect_equal( + nhdplusTools:::resolve_comid_from_reachcode("14070003001234", 65, vaa), + 102L + ) + # measure on boundary -> first matching interval + expect_equal( + nhdplusTools:::resolve_comid_from_reachcode("14070003001234", 50, vaa), + 101L + ) + # measure out of bounds (>100) -> closest interval boundary (comid 103, tomeas=100) + expect_equal( + nhdplusTools:::resolve_comid_from_reachcode("14070003001234", 110, vaa), + 103L + ) + # missing reachcode -> error + expect_error( + nhdplusTools:::resolve_comid_from_reachcode("00000000000000", 50, vaa), + "No flowlines found in VAA" + ) +}) + +test_that("get_drainage_area_estimates accepts reachcode/measure in start", { + skip_on_cran() + + # source the actual reachcode/measure from a one-time NLDI call + ref <- get_nldi_feature( + list(featureSource = "nwissite", featureID = "USGS-05406500") + ) + ref_measure <- as.numeric(ref$measure) + + # 1. override path: NLDI still called for COMID, reachcode/measure overridden + override_measure <- ref_measure - 5 + start_ovr <- list( + featureSource = "nwissite", featureID = "USGS-05406500", + reachcode = ref$reachcode, measure = override_measure + ) + res_ovr <- suppressWarnings(suppressMessages( + get_drainage_area_estimates(start_ovr, nhdplushr = FALSE) + )) + expect_equal(res_ovr$start_feature$measure, override_measure) + expect_equal(res_ovr$start_feature$reachcode, ref$reachcode) + expect_true(res_ovr$da_huc12_sqkm > 0) + expect_true(res_ovr$network_da_sqkm > 0) + + # 2. no-NLDI path: only reachcode/measure provided; COMID resolved via OGC, + # start point derived from the flowline at the supplied measure + start_no_nldi <- list( + reachcode = ref$reachcode, measure = ref_measure + ) + res_no_nldi <- suppressWarnings(suppressMessages( + get_drainage_area_estimates(start_no_nldi, nhdplushr = FALSE) + )) + expect_equal(as.integer(res_no_nldi$start_feature$comid), + as.integer(ref$comid)) + expect_equal(res_no_nldi$start_feature$reachcode, ref$reachcode) + expect_equal(res_no_nldi$start_feature$measure, ref_measure) + expect_true(res_no_nldi$da_huc12_sqkm > 0) + expect_equal(res_no_nldi$network_da_sqkm, res_ovr$network_da_sqkm) + + # derived start point should be a real POINT near the NLDI-reported location + geom <- sf::st_geometry(res_no_nldi$start_feature) + expect_s3_class(geom, "sfc_POINT") + expect_false(sf::st_is_empty(geom)[1]) + dist_m <- as.numeric(sf::st_distance( + sf::st_transform(geom, 5070), + sf::st_transform(sf::st_geometry(ref), 5070) + )) + expect_lt(dist_m, 500) +}) + test_that("get_drainage_area_estimates local_navigation smoke test", { skip_on_cran()