# 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?

• I think until something like `st_nearest` is implemented, it might be best to convert to `sp` and use `rgeos`: `cutpt = st_as_sf(rgeos::gNearestPoints(as(site,"Spatial"), as(reach,"Spatial"))[2,])` - the resulting point has zero distance to the line and `st_split` will split it. Jul 5, 2018 at 21:25

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`.

• Have you found a resolution to this? I am in exactly the same boat, and it is also driving me crazy. Feb 20, 2019 at 16:18
• @mbcardamia Yeah, I had seen your GitHub issue. My ultimate workaround was to not use an R solution altogether and use ArcPy (which I have access to), but I wish I could keep all of my analysis in R. I understand that floating point precision is a bear, but it does seem crazy that there's no workaround for a use case like this. Feb 26, 2019 at 14:36
• To everyone finding this thread in 2022, I have a pretty dumb but effective solution! If blah = st_nearest_points(x, y) followed by blahblah = st_cast(blah, "POINT") is giving you problems with st_intersects(blahblah, y) because of floating point math/really small distances between blahblah and y, just throw st_buffer() around y inside st_nearest_points() with a really small negative value, ala st_nearest_points(x, st_buffer(y, -0.001)). This will force the "nearest point" to move slightly further into the y feature so that blahblah will then safely intersect it! May 4, 2022 at 21:36

Here's an outline method. Maybe code later but this should get you started:

`st_distance(site, reach)` gets the distance from the site to the reach. So a buffer with this radius should touch the reach at that point (or set of points) and so you could then use `st_intersection` to get the intersection of the buffer and that must be the nearest point.

Numerical precision bites you though, and you may find the buffer doesn't intersect the reach, so add a small tolerance value to the buffer radius. So...

``````tol = 1
units(tol)="m"
r = st_cast(st_buffer(site, tol+st_distance(site,reach[1,1])),"LINESTRING")
``````

I cast `r` to a LINESTRING so that its not a POLYGON any more. Then when you intersect you get a LINESTRING and not a chunk of a circle:

``````> nearest = st_intersection(r, reach)
> nearest\$geom
Geometry set for 1 feature
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: 688447.1 ymin: 172078.9 xmax: 688447.6 ymax: 172092.1
epsg (SRID):    NA
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
MULTIPOINT (688447.1 172078.9, 688447.6 172092.1)
``````

Those two points are very close to the true nearest point, and `st_cast(nearest, "POINT")` gets you the two points.

However those points don't seem to `st_split` the line, which is what you want to do, so maybe we need to get `st_snap` working.

• I've opened an issue in `sf` to implement geos' gNearestPoints: github.com/r-spatial/sf/issues/788 Jul 6, 2018 at 11:19
• I think such a function would be very useful, but the "ideal" implementation would be more like the Near tool in ArcGIS, which can mix geometries, and returns distances and near locations. Incidentally, you refer to the st_snap not working. Is it a known issue in `sf`? Because it doesn't always work for all sites near a reach in this case. Jul 6, 2018 at 11:28
• I think I understand st_snap - it doesn't snap to the nearest point on the second geom, but to a point on the geom on a grid of tolerance. I think... Jul 6, 2018 at 13:31
• gNearestPoints can mix geometries, and once you have the points gDistance will get you the distance. Jul 6, 2018 at 13:32

One approach is to convert the line into points and find the nearest point. Depending on trade-off between accuracy and processing time, you can select the optimal density of points to create.

Here, I create a point every 1m along the line. Use mapview to check.

``````library(mapview)
points <- st_line_sample(reach, density = 1/1) %>%
st_sf() %>%
st_cast('POINT')
nearest <- points[which.min(st_distance(site, points)),]
mapview(nearest) + site + reach
``````

Unfortunately even when points are created along a line using `st_line_sample()`, they do not intersect the line. I am guessing this is due to rounding. You might have to buffer the point slightly to intersect with the reach.

Going back to your original question (splitting the line at the nearest point), you could then recreate the two segments from the sampled points, like so:

``````points <- st_line_sample(reach, density = 1/1) %>%
st_sf() %>%
st_cast('POINT') %>%
mutate(group = 1)

which.point <- which.min(st_distance(site, points))

nearest <- points[which.point,]

segment1 <- points[1:which.point,] %>%
group_by(group) %>%
summarise(do_union = FALSE) %>%
st_cast("LINESTRING") %>%
ungroup

segment2 <- points[which.point:nrow(points),] %>%
group_by(group) %>%
summarise(do_union = FALSE) %>%
st_cast("LINESTRING") %>%
ungroup
``````

mapview(segment1, color = 'red') + segment2 + nearest + site

Finally calculate length using st_length:

``````segment1 %<>% mutate(Length = st_length(.))
``````
• Thank you for your suggestion. Sampling the polyline and finding the nearest point could work, however it only approximates the nearest location along the reach. Wouldn't st_cast(reach\$geometry, "POINTS") find the exact closest point to each site? Jul 5, 2018 at 17:10
• st_cast(reach\$geometry, "POINTS") will simply convert the linestring to its vertices. So if there are long segments with no vertices the nearest point may actually be very inaccurate. Jul 5, 2018 at 17:24
• see my edits addressing the issue of splitting the line at nearest point. Jul 5, 2018 at 19:36

I have a similar problem: I need to snap a point on a linestring. Even though I don't split the linestrings, I hope my solution would be useful for you.

My approach is this: Due to the rounding errors, the nearest point needs not to be precisely on a linestring. My solution is to add a new point to the linestring that is precisely the nearest point. (Splitting the line here should be trivial then.)

My code is this:

``````library(sf)
library(dplyr)
library(purrr)
library(tibble)

np <- st_cast(st_nearest_points(point, nf), "POINT")[2]
nf_points <- st_cast(st_geometry(nf), "POINT")
if (any(st_distance(nf_points, np) < tolerance))
return(NULL)
ppg <- st_coordinates(nf_points)
res <- map(1:(nrow(ppg) - 1),
~ st_set_crs(st_geometry(st_linestring(ppg[.:(. + 1),])),
map_dbl(~ units::drop_units(st_distance(.,  np)))
row <- which(res == min(res))[1]
new_ppg <- rbind(ppg[1:row, ], st_coordinates(np), ppg[row:nrow(ppg), ])
new_geometry <- new_ppg %>% st_linestring() %>% st_geometry()
list(idx = idx_of_nf, geometry = new_geometry)
}

for (k in seq_len(nrow(points))) {
point <- points[k, ]
if (!is.null(ng))
}