Find nearest point along polyline using 'sf' package in R

This is a follow-up question related to my earlier post about using the sf package to split polylines using nearby points.

I have 450 points with XY coordinates and a large river network (~214,000 reaches). Many of the points do not intersect the river network, so for each point I would like to find the nearest location and distance along a reach, and to also snap the point to that location.

My current implementation uses a custom wrapper function calling geosphere::dist2line() (largely based on the following questions by Guzman and Ana; see my first question) to do the following (I've include one simple feature POINT (site) and LINESTRING or polyline (reach) to provide a reproducible example):

library(sf)
site <- structure(list(SiteId = "PF_00183", SampId = "PF_00183 _ 2014-10-15", 
    MonitoringProgram = "ProgettoFiumi", X1 = 688524.841170469, 

    Y1 = 172082.839248429, X2 = 688405.761917025, Y2 = 172008.150777719, 
    Width = NA_real_, TransectDistance = 140.563993461242, Area = NA_real_, 
    EZG_NR = 54757, h1 = 163815, h2 = 168422, geometry = structure(list(
        structure(c(688524.8412, 172082.8392), class = c("XY", 
        "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(688524.8412, 
    172082.8392, 688524.8412, 172082.8392), .Names = c("xmin", 
    "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
        epsg = NA_integer_, proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("SiteId", 
"SampId", "MonitoringProgram", "X1", "Y1", "X2", "Y2", "Width", 

"TransectDistance", "Area", "EZG_NR", "h1", "h2", "geometry"), row.names = c(NA, 
-1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = c("SiteId", "SampId", "MonitoringProgram", 
"X1", "Y1", "X2", "Y2", "Width", "TransectDistance", "Area", 
"EZG_NR", "h1", "h2")))

reach <- structure(list(FNODE_ = 148558, TNODE_ = 148142, LENGTH = 343.00790968, 

    ReachId = 125010L, geometry = structure(list(structure(c(688400.125, 
    688394.375, 688385.875, 688375.6875, 688363.625, 688359.8125, 
    688358.8125, 688361.125, 688375.6875, 688389.375, 688411.1875, 
    688432.875, 688445.1875, 688448.125, 688445.875, 688446.625, 
    688448.125, 688455.375, 688465, 171876.40625, 171890.40625, 
    171906.203125, 171916.703125, 171929.09375, 171944.703125, 
    171952.59375, 171959.90625, 171980.296875, 171997.09375, 
    172017.796875, 172041.90625, 172063.5, 172086.703125, 172109.90625, 
    172129, 172138.5, 172153.5, 172171.203125), .Dim = c(19L, 
    2L), class = c("XY", "LINESTRING", "sfg"))), class = c("sfc_LINESTRING", 

    "sfc"), precision = 0, bbox = structure(c(688358.8125, 171876.40625, 
    688465, 172171.203125), .Names = c("xmin", "ymin", "xmax", 
    "ymax"), class = "bbox"), crs = structure(list(epsg = NA_integer_, 
        proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("FNODE_", 
"TNODE_", "LENGTH", "ReachId", "geometry"), row.names = 125010L, class = c("sf", 
"data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_), .Names = c("FNODE_", 
"TNODE_", "LENGTH", "ReachId"), .Label = c("constant", "aggregate", 
"identity"), class = "factor"))


# Transform to CRS WGS1984 (EPSG: 4326) and create 'sp' inputs for geosphere::dist2line()
site.sp <- st_transform(site, crs = 4326)
site.sp <- as_Spatial(st_geometry(site.sp))

reach.sp <- st_transform(reach, crs = 4326)
reach.sp <- as_Spatial(st_zm(st_geometry(reach.sp)))

# Calculate the location and distance of the nearest point along the polyline
d <- geosphere::dist2Line(site.sp, reach.sp)

d <- as.data.frame(d)

# Create simple feature with CRS WGS1984 (unprojected)
site.near <- st_as_sf(d, coords=c("lon", "lat"), crs=4326)

# Transform CRS from WGS1984 (EPSG: 4326) to CH1903/LV03
site.near <- st_transform(site.near, crs=st_crs(site))

# Check if the nearest point intersects the reach (FALSE)
st_intersects(site.near, reach, sparse = F)

The problem with this implementation is that (1) the nearest point along the reach does not intersect the reach polyline, and (2) it switches between the geosphere and sf packages, sp and sf object classes, and transforms between CRSs. In short, it doesn't work and its too complicated.

For these reasons, I would prefer to use only on the sf package to accomplish this entire task. I have some intuition that st_snap() would be useful here, but it always returns the same point as the input (the PostGIS and sf documentation pages don't really explain what the tolerance argument does).

EDIT As suggested by Spacedman's comment I used rgeos::gNearestPoints. It provided a the nearest point intersecting the reach for the site/reach I provided above, however for another site it provides a nearest point that is still a small distance from the reach:

> st_distance(site.near, reach)
Units: m
            [,1]
[1,] 3.81693e-09

I provide the site and reach below:

site <- structure(list(SiteId = "NAWA_2", SampId = "NAWA_2_2012-08-13", 
    MonitoringProgram = "NAWA", X1 = 613474, Y1 = 267088, X2 = 613603, 
    Y2 = 266936, Width = 21.87, TransectDistance = 199.361480732864, 
    Area = 4360.03558362773, EZG_NR = 91795, h1 = 74792, h2 = 79777, 
    geometry = structure(list(structure(c(613474, 267088), class = c("XY", 
    "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(613474, 
    267088, 613474, 267088), .Names = c("xmin", "ymin", "xmax", 
    "ymax"), class = "bbox"), crs = structure(list(epsg = NA_integer_, 
        proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("SiteId", 

"SampId", "MonitoringProgram", "X1", "Y1", "X2", "Y2", "Width", 
"TransectDistance", "Area", "EZG_NR", "h1", "h2", "geometry"), row.names = c(NA, 
-1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = c("SiteId", "SampId", "MonitoringProgram", 
"X1", "Y1", "X2", "Y2", "Width", "TransectDistance", "Area", 
"EZG_NR", "h1", "h2")))


reach <- structure(list(OBJECTID_1 = 5381, FNODE_ = 15469, TNODE_ = 13378, 
    LPOLY_ = 0, RPOLY_ = 0, LENGTH = 3248.0468712, OBJECTID = 1517084, 
    OBJECTVAL = "Fluss", NAME = "Birs", UNTERIRDIS = NA_character_, 
    OBJECTORIG = "GN25", YEAROFCHAN = 2005, GEWISSNR = 443, LAUFNR = 0, 
    LINST = "CH", BACHNR = NA_character_, GWLNR = "CH0004430000", 
    Binary = 1L, Shape_Leng = 3248.0468712, ReachId = 5381L, 
    SiteId = "PF_00449", SampId = "PF_00449 _ 2014-09-25", MonitoringProgram = "ProgettoFiumi", 
    X1 = 613481.575730715, Y1 = 267350.715242873, X2 = NA_real_, 
    Y2 = NA_real_, Width = NA_real_, TransectDistance = NA_real_, 
    Area = NA_real_, EZG_NR = 91795, h1 = 74792, h2 = 79777, 

    SiteId.1 = "NAWA_2", SampId.1 = "NAWA_2_2012-08-13", MonitoringProgram.1 = "NAWA", 
    X1.1 = 613474, Y1.1 = 267088, X2.1 = 613603, Y2.1 = 266936, 
    Width.1 = 21.87, TransectDistance.1 = 199.361480732864, Area.1 = 4360.03558362773, 
    EZG_NR.1 = 91795, h1.1 = 74792, h2.1 = 79777, geometry = structure(list(
        structure(c(613680.112152525, 613672, 613574.625, 613541.625, 
        613535.1875, 613507.125, 613490.1875, 613476.1875, 613472.375, 
        613471.5, 613472.375, 613500.625, 613500.625, 613485.625, 
        613473.8125, 613432.1875, 266891.941045502, 266898.3125, 
        266968.59375, 267000, 267006.09375, 267039.8125, 267070.6875, 
        267106.3125, 267128.8125, 267151.3125, 267177.59375, 

        267405.3125, 267431.5, 267461.5, 267487.1875, 267578.3125
        ), .Dim = c(16L, 2L), class = c("XY", "LINESTRING", "sfg"
        ))), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(613432.1875, 
    266891.941045502, 613680.112152525, 267578.3125), .Names = c("xmin", 
    "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
        epsg = NA_integer_, proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("OBJECTID_1", 
"FNODE_", "TNODE_", "LPOLY_", "RPOLY_", "LENGTH", "OBJECTID", 
"OBJECTVAL", "NAME", "UNTERIRDIS", "OBJECTORIG", "YEAROFCHAN", 
"GEWISSNR", "LAUFNR", "LINST", "BACHNR", "GWLNR", "Binary", "Shape_Leng", 

"ReachId", "SiteId", "SampId", "MonitoringProgram", "X1", "Y1", 
"X2", "Y2", "Width", "TransectDistance", "Area", "EZG_NR", "h1", 
"h2", "SiteId.1", "SampId.1", "MonitoringProgram.1", "X1.1", 
"Y1.1", "X2.1", "Y2.1", "Width.1", "TransectDistance.1", "Area.1", 
"EZG_NR.1", "h1.1", "h2.1", "geometry"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 

NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_
), .Names = c("OBJECTID_1", "FNODE_", "TNODE_", "LPOLY_", "RPOLY_", 
"LENGTH", "OBJECTID", "OBJECTVAL", "NAME", "UNTERIRDIS", "OBJECTORIG", 
"YEAROFCHAN", "GEWISSNR", "LAUFNR", "LINST", "BACHNR", "GWLNR", 
"Binary", "Shape_Leng", "ReachId", "SiteId", "SampId", "MonitoringProgram", 
"X1", "Y1", "X2", "Y2", "Width", "TransectDistance", "Area", 
"EZG_NR", "h1", "h2", "SiteId.1", "SampId.1", "MonitoringProgram.1", 

"X1.1", "Y1.1", "X2.1", "Y2.1", "Width.1", "TransectDistance.1", 
"Area.1", "EZG_NR.1", "h1.1", "h2.1"), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = "5381.5", class = c("sf", 
"data.frame"))

Obtaining the nearest point along the reach gives:

site.near <- st_as_sf(rgeos::gNearestPoints(as(site,"Spatial"), as(st_zm(reach),"Spatial"))[2,])
> st_intersects(site.near, reach, sparse=FALSE)
      [,1]
[1,] FALSE

> st_distance(site.near, reach)
Units: m
            [,1]
[1,] 3.81693e-09

Using st_snap does not seem to work:

> site.snap <- st_snap(site.near, reach, tol=1e-9)
> st_intersects(site.near, reach, sparse=FALSE)
      [,1]
[1,] FALSE

Attempting to use maptools::snapPointsToLines also fails to intersect the nearest point to the reach:

> site.near <- st_as_sf(maptools::snapPointsToLines(as(site,"Spatial"), as(st_zm(reach),"Spatial")))
> st_intersects(site.near, reach, sparse=FALSE)
      [,1]
[1,] FALSE
> st_distance(site.near, reach)
Units: m
             [,1]
[1,] 3.269482e-11

Could this be an issue with numerical precision in R?

Answer

This is a longer post containing what I've found so far: installing sf_0.6-4 and making use of the new st_nearest_feature and st_nearest_points functions will locate the nearest point along a polyline using only sf. Here is the implementation with sf_0.6-4:

library(sf)

site <- structure(list(SiteId = "NAWA_2", SampId = "NAWA_2_2012-08-13", 
    MonitoringProgram = "NAWA", X1 = 613474, Y1 = 267088, X2 = 613603, 
    Y2 = 266936, Width = 21.87, TransectDistance = 199.361480732864, 
    Area = 4360.03558362773, EZG_NR = 91795, h1 = 74792, h2 = 79777, 

    geometry = structure(list(structure(c(613474, 267088), class = c("XY", 
    "POINT", "sfg"))), class = c("sfc_POINT", "sfc"), precision = 0, bbox = structure(c(613474, 
    267088, 613474, 267088), .Names = c("xmin", "ymin", "xmax", 
    "ymax"), class = "bbox"), crs = structure(list(epsg = NA_integer_, 
        proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("SiteId", 
"SampId", "MonitoringProgram", "X1", "Y1", "X2", "Y2", "Width", 
"TransectDistance", "Area", "EZG_NR", "h1", "h2", "geometry"), row.names = c(NA, 
-1L), class = c("sf", "data.frame"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 

NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_), class = "factor", .Label = c("constant", 
"aggregate", "identity"), .Names = c("SiteId", "SampId", "MonitoringProgram", 
"X1", "Y1", "X2", "Y2", "Width", "TransectDistance", "Area", 
"EZG_NR", "h1", "h2")))

reach <- structure(list(OBJECTID_1 = 5381, FNODE_ = 15469, TNODE_ = 13378, 
    LPOLY_ = 0, RPOLY_ = 0, LENGTH = 3248.0468712, OBJECTID = 1517084, 
    OBJECTVAL = "Fluss", NAME = "Birs", UNTERIRDIS = NA_character_, 
    OBJECTORIG = "GN25", YEAROFCHAN = 2005, GEWISSNR = 443, LAUFNR = 0, 

    LINST = "CH", BACHNR = NA_character_, GWLNR = "CH0004430000", 
    Binary = 1L, Shape_Leng = 3248.0468712, ReachId = 5381L, 
    SiteId = "PF_00449", SampId = "PF_00449 _ 2014-09-25", MonitoringProgram = "ProgettoFiumi", 
    X1 = 613481.575730715, Y1 = 267350.715242873, X2 = NA_real_, 
    Y2 = NA_real_, Width = NA_real_, TransectDistance = NA_real_, 
    Area = NA_real_, EZG_NR = 91795, h1 = 74792, h2 = 79777, 
    SiteId.1 = "NAWA_2", SampId.1 = "NAWA_2_2012-08-13", MonitoringProgram.1 = "NAWA", 
    X1.1 = 613474, Y1.1 = 267088, X2.1 = 613603, Y2.1 = 266936, 
    Width.1 = 21.87, TransectDistance.1 = 199.361480732864, Area.1 = 4360.03558362773, 
    EZG_NR.1 = 91795, h1.1 = 74792, h2.1 = 79777, geometry = structure(list(

        structure(c(613680.112152525, 613672, 613574.625, 613541.625, 
        613535.1875, 613507.125, 613490.1875, 613476.1875, 613472.375, 
        613471.5, 613472.375, 613500.625, 613500.625, 613485.625, 
        613473.8125, 613432.1875, 266891.941045502, 266898.3125, 
        266968.59375, 267000, 267006.09375, 267039.8125, 267070.6875, 
        267106.3125, 267128.8125, 267151.3125, 267177.59375, 
        267405.3125, 267431.5, 267461.5, 267487.1875, 267578.3125
        ), .Dim = c(16L, 2L), class = c("XY", "LINESTRING", "sfg"
        ))), class = c("sfc_LINESTRING", "sfc"), precision = 0, bbox = structure(c(613432.1875, 
    266891.941045502, 613680.112152525, 267578.3125), .Names = c("xmin", 

    "ymin", "xmax", "ymax"), class = "bbox"), crs = structure(list(
        epsg = NA_integer_, proj4string = "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +units=m +no_defs"), .Names = c("epsg", 
    "proj4string"), class = "crs"), n_empty = 0L)), .Names = c("OBJECTID_1", 
"FNODE_", "TNODE_", "LPOLY_", "RPOLY_", "LENGTH", "OBJECTID", 
"OBJECTVAL", "NAME", "UNTERIRDIS", "OBJECTORIG", "YEAROFCHAN", 
"GEWISSNR", "LAUFNR", "LINST", "BACHNR", "GWLNR", "Binary", "Shape_Leng", 
"ReachId", "SiteId", "SampId", "MonitoringProgram", "X1", "Y1", 
"X2", "Y2", "Width", "TransectDistance", "Area", "EZG_NR", "h1", 
"h2", "SiteId.1", "SampId.1", "MonitoringProgram.1", "X1.1", 
"Y1.1", "X2.1", "Y2.1", "Width.1", "TransectDistance.1", "Area.1", 

"EZG_NR.1", "h1.1", "h2.1", "geometry"), sf_column = "geometry", agr = structure(c(NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_, 
NA_integer_, NA_integer_, NA_integer_, NA_integer_, NA_integer_

), .Names = c("OBJECTID_1", "FNODE_", "TNODE_", "LPOLY_", "RPOLY_", 
"LENGTH", "OBJECTID", "OBJECTVAL", "NAME", "UNTERIRDIS", "OBJECTORIG", 
"YEAROFCHAN", "GEWISSNR", "LAUFNR", "LINST", "BACHNR", "GWLNR", 
"Binary", "Shape_Leng", "ReachId", "SiteId", "SampId", "MonitoringProgram", 
"X1", "Y1", "X2", "Y2", "Width", "TransectDistance", "Area", 
"EZG_NR", "h1", "h2", "SiteId.1", "SampId.1", "MonitoringProgram.1", 
"X1.1", "Y1.1", "X2.1", "Y2.1", "Width.1", "TransectDistance.1", 
"Area.1", "EZG_NR.1", "h1.1", "h2.1"), .Label = c("constant", 
"aggregate", "identity"), class = "factor"), row.names = "5381.5", class = c("sf", 
"data.frame"))


reach.nearest <- st_nearest_feature(site, reach)
reach.nearest <- reach[reach.nearest,]

# Locate nearest point along nearest reach (Returns line)
site.nearest <- st_nearest_points(site, reach.nearest)
cat("Line intersects:", st_intersects(site.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
st_distance(site.nearest, reach.nearest)

# st_nearest_points returns lines, so use st_cast

site.nearest <- st_cast(site.nearest, "POINT")[2]

# Test whether the nearest point intersects the reach
cat("Point intersects:", st_intersects(site.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
st_distance(site.nearest, reach.nearest)

The problem is that st_nearest_points returns a line connecting the site to the nearest point along the reach. The line has a distance of 0m and does intersect the reach. However, extracting the point with st_cast does not intersect the reach:

> # st_nearest_points returns lines, so use st_cast
> site.nearest <- st_cast(site.nearest, "POINT")[2]
> 

> # Test whether the nearest point intersects the reach
> cat("Point intersects:", st_intersects(site.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
Point intersects: FALSE |> st_distance(site.nearest, reach.nearest)
Units: m
             [,1]
[1,] 3.269482e-11

This minuscule distance is the only thing preventing me from intersecting the point, and its driving me crazy. The underlying reason for this the use of floating point coordinates and their precision. See here for an ongoing issue on GitHub. According to the documentation of sf (see vignette 1, precision section):

One of the attributes of a geometry list-column (sfc) is the precision: a double number that, when non-zero, causes some rounding during conversion to WKB, which might help certain geometrical operations succeed that would otherwise fail due to floating point representation. The model is that of GEOS, which copies from the Java Topology Suite (JTS), and works like this:

• if precision is zero (default, unspecified), nothing is modified


• negative values convert to float (4-byte real) precision


• positivevalues convert to round(x*precision)/precision.

For the precision model, see also here, where it is written that: “… to specify 3 decimal places of precision, use a scale factor of 1000. To specify -3 decimal places of precision (i.e. rounding to the nearest 1000), use a scale factor of 0.001.”

Setting the precision as an sf attribute will not change the coordinates of the sf object until they it is used in GEOS for doing geometric operations (source). From the doc a precision of 1000 (3 decimal places of precision) seems reasonable, but returns this:

> sp <- function(sf, precision=1000){
+   st_set_precision(sf, precision)

+ }
> 
> site <- sp(site, 1000)
> reach <- sp(reach, 1000)
> 
> reach.nearest <- st_nearest_feature(site, reach)
> reach.nearest <- reach[reach.nearest,]
> 
> # Locate nearest point along nearest reach (Returns line)
> site.nearest <- st_nearest_points(site, reach.nearest)

> cat("Line intersects:", st_intersects(site.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
Line intersects: FALSE |> st_distance(site.nearest, reach.nearest)
Units: m
             [,1]
[1,] 0.0001012353
> 
> # st_nearest_points returns lines, so use st_cast
> site.nearest <- st_cast(site.nearest, "POINT")[2]
> 
> # Test whether the nearest point intersects the reach

> cat("Point intersects:", st_intersects(site.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
Point intersects: FALSE |> st_distance(site.nearest, reach.nearest)
Units: m
             [,1]
[1,] 0.0001012353

At this point I'm stuck with understanding why the precision argument does not work, but this post at least addresses the original question of using only the sf package. I may open a new question to fully address the problem of precision.

TLDR: use new geometric functions in sf version 0.6.4 to find nearest point, distances of nearest point is still 3.2e-11m away from river due to precision of floating point coordinates. Set precision accordingly with st_set_precision.