This is a follow-up question related to my [earlier post about using the `sf` package to split polylines using nearby points](https://gis.stackexchange.com/questions/288445/split-polyline-by-point-using-sf-package-in-r).

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](https://gis.stackexchange.com/questions/243994/how-to-calculate-distance-from-point-to-linestring-in-r-using-sf-library-and-g) and [Ana](https://gis.stackexchange.com/questions/269746/find-nearest-line-segment-to-each-point-in-r); 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](https://postgis.net/docs/ST_Snap.html) and [`sf` documentation](https://r-spatial.github.io/sf/reference/geos_binary_ops.html) 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?