COMID errors from get_flowline_index

[zcjohnson3]

On occasion, the get_flowline_index function returns an incorrect COMID assignment (and therefore an erroneous REACH_meas value) for coordinates that are clearly in a nearby COMID and on or near the NHD flowline. Below is a reproducible example that returns COMID 23809098 & REACH_meas 52.7, putting the site just downstream of the North Fork Dam on the Clackamas River in Oregon. However, these coordinates are actually just upstream of the dam in COMID 23809100 (REACH_meas ~0-10).

nhdplusTools::get_flowline_index('download_nhdplusv2', sf::st_sfc(sf::st_point(c(-122.278972,45.243732)), crs=4326))

[dblodgett-usgs]

Hi @zcjohnson3 --

This is a hard one. It may just be a limitation of the method, but I think I at least have a work around for you.

If you were to look at the distance from your point to the geometry nodes, it is equidistant to the outlet of the flowline it is along and the flowline that it is upstream of. Notice that the portion of the line that your point appears to touch is straight -- this indicates that it is a segment between two geometry nodes. nhdplusTools is matching to the geometry node locations -- not the projected line connecting them (by default)

If you set precision -- it will add geometry nodes and get you somewhat better results.

point <- sf::st_sfc(sf::st_point(c(-122.278972,45.243732)), crs=4326)

index_1 <- nhdplusTools::get_flowline_index(
  'download_nhdplusv2', 
  point
)
#> Warning in st_buffer.sfc(points, search_radius): st_buffer does not correctly
#> buffer longitude/latitude data

index_2 <- nhdplusTools::get_flowline_index(
  'download_nhdplusv2', 
  point, 
  precision = 10
)
#> Warning in st_buffer.sfc(points, search_radius): st_buffer does not correctly
#> buffer longitude/latitude data
#> although coordinates are longitude/latitude, st_intersects assumes that they
#> are planar
#> Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is FALSE

index_1  
#>   id    COMID      REACHCODE REACH_meas      offset
#> 1  1 23809098 17090011000083    52.7117 0.001746542
index_2
#>   id    COMID      REACHCODE REACH_meas       offset
#> 1  1 23809100 17090011000083    55.5661 0.0001060655

^{Created on 2026-03-13 with reprex v2.1.1}

I don't really think there's a fix here unfortunately. Do you think there's anything else we could do to make this more likely to succeed given the data?

[zcjohnson3]

Thanks, @dblodgett-usgs . I had overlooked the precision input for that function. Does the function run slower if you give it a precision value? Since I'm looking to do this for many sites, I'm wondering if it would be best to try it first with no precision input but if the resulting COMID does not match what I expect then run again with precision, OR if I should just always run it with a precision input to guard against any erroneous COMID results. Also, any recommendations for precision values (or is 10 a good default)?

[dblodgett-usgs]

Hi -- It may slow down, but more than anything it can use a lot more memory. But if you keep the precision relatively course (50m is probably fine) it's not going to add much except in places like you are having this problem. Maybe try doing it everywhere with a 30m or 50m tolerance and see if performance is ok?