Snap a point to the closest point on a line segmen

Question

I would like to snap a point to the closest point on a road segment using sf::st_snap. However, the function seems to return the wrong result, it's snapping my point to a point at the start of the road segment. Any ideas on how to fix this?

Reproducible example provided below, including a comparision of the results I get when using sf::st_snap vs maptools::snapPointsToLines

Using sf::st_snap

# Max distance
  cut_dist = 200 # meters

# snap points to closest road
  new_point <- sf::st_snap(point1, roads, tolerance = cut_dist)

# view points on the map
  mapView(point1, color="red") + mapView( st_buffer(point1, dist = cut_dist)) + mapView(new_point) + mapView(roads) 

# Distance between pont1 and new_point
  st_distance( point1, new_point)
> 1591 meters # note this is above the set maximun distance

Using maptools::snapPointsToLines (the result I would expect)

# convert sf to sp
  point_sp <-  as_Spatial(point1)
  roads_sp <-  as_Spatial(roads)

# snap points
  new_point_sp <- snapPointsToLines(point_sp, roads_sp, maxDist = cut_dist)

# view points on the map
  mapView(point1, color="red") + mapView( st_buffer(point1, dist = cut_dist)) + mapView(new_point_sp) + mapView(roads) 

# Distance between pont1 and new_point
  spDistsN1( point_sp, new_point_sp)
>  116 meters

Data and libraries

library(sf)
library(mapview)
library(maptools)
library(sp)

point1 <- structure(list(idhex = 9L, geometry = structure(list(structure(c(665606.970079183, 
          6525003.41418009), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT", 
          "sfc"), precision = 0, bbox = structure(c(xmin = 665606.970079183, 
          ymin = 6525003.41418009, xmax = 665606.970079183, ymax = 6525003.41418009
          ), class = "bbox"), crs = structure(list(epsg = 32633L, proj4string = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"), class = "crs"), n_empty = 0L)), sf_column = "geometry", agr = structure(c(idhex = NA_integer_), .Label = c("constant", 
                                                                                                                                                         "aggregate", "identity"), class = "factor"), row.names = 2L, class = c("sf", 
                                                                                                                                                                                                                                                                                                                                                                                                                                                              "data.table", "data.frame"))


   roads <- structure(list(id = 139885, osm_id = 250886593, geometry = structure(list(
        structure(c(665387.589147313, 665367.867159783, 665363.008143169, 
        665363.051860059, 665363.308104069, 665366.519781353, 665368.635421323, 
        665370.846894641, 665370.829724196, 665367.910645335, 665361.777524054, 
        665355.967776345, 665351.649946698, 665343.44353825, 665334.917779131, 
        665313.306069501, 665309.001351385, 665310.66019677, 665313.528620709, 
        665341.742306731, 665351.854389331, 665354.981775569, 665360.254611229, 
        665365.006104512, 665379.034588542, 665394.435039616, 665409.282519288, 
        665410.676785182, 665448.890797438, 665458.917562631, 665471.042094353, 
        665485.485001236, 665495.899212422, 665504.535684257, 665509.674854913, 
        665506.145837246, 665483.727146874, 665481.426949686, 665462.311063365, 
        665445.215460573, 665450.424049663, 665450.837897892, 665491.036360788, 
        665491.419140717, 665469.507518623, 665458.677850808, 665455.926197775, 
        665462.873809047, 665460.283684337, 665426.046702616, 665396.279686035, 
        665368.373253059, 665357.878521323, 665304.347529357, 665221.04051616, 
        665170.777462125, 665144.670345016, 665106.030568334, 665073.2789218, 
        665018.208956171, 664947.693178271, 664921.708297412, 664861.659061389, 
        664797.900403384, 664745.001666066, 664730.200174759, 664717.892651619, 
        664706.473711845, 664697.750102392, 664688.215719591, 664681.544531593, 
        664672.960647368, 664665.064067202, 664636.446517023, 664622.930521655, 
        664518.065243846, 664442.725560545, 664423.048166559, 664411.132259582, 
        664407.05972929, 664398.364646172, 664391.348502443, 664382.558239303, 
        664372.012526058, 664354.354954718, 664332.995014599, 664311.609706282, 
        664271.102641808, 664228.816287751, 664150.088321471, 664069.895400484, 
        6526138.02793883, 6526135.40749336, 6526130.11578605, 6526111.34403368, 
        6526087.4978365, 6526054.13445288, 6526022.49962268, 6525982.74830288, 
        6525959.40435839, 6525944.55197219, 6525918.33886077, 6525894.18611795, 
        6525874.55473851, 6525840.53410542, 6525813.96628006, 6525767.42907088, 
        6525745.21917638, 6525733.51582599, 6525713.24841331, 6525627.57847652, 
        6525608.06984863, 6525568.30170735, 6525550.71644271, 6525539.76231607, 
        6525491.25651378, 6525446.12690364, 6525433.36256694, 6525431.23562504, 
        6525372.98235432, 6525354.13376808, 6525331.3288195, 6525309.59511696, 
        6525293.92174422, 6525270.21980161, 6525256.11455612, 6525228.35885783, 
        6525217.10943051, 6525215.95489587, 6525195.91355696, 6525158.79257025, 
        6525134.01851773, 6525131.70940566, 6525050.96446632, 6524950.68358502, 
        6524851.23226232, 6524806.24052727, 6524749.34394609, 6524714.63617193, 
        6524660.07336072, 6524612.21010524, 6524583.84484865, 6524562.03540982, 
        6524557.38094998, 6524533.67136837, 6524510.74454804, 6524495.56823805, 
        6524486.9387399, 6524475.63373441, 6524465.4404841, 6524468.04929815, 
        6524475.95178632, 6524478.86036788, 6524470.76472937, 6524447.96214429, 
        6524448.06967557, 6524443.4855897, 6524435.86812114, 6524425.93373791, 
        6524417.67487537, 6524409.79262886, 6524399.64960133, 6524378.79085156, 
        6524360.33496349, 6524303.24355601, 6524302.70486651, 6524293.01335665, 
        6524290.81442892, 6524298.30279414, 6524309.46697681, 6524313.27442914, 
        6524337.22831533, 6524364.43083297, 6524376.27944935, 6524382.92319852, 
        6524389.6474774, 6524406.74565716, 6524430.82326744, 6524462.46041311, 
        6524492.20009833, 6524544.74318075, 6524591.10483188), .Dim = c(91L, 
        2L), class = c("XY", "LINESTRING", "sfg"))), class = c("sfc_LINESTRING", 
        "sfc"), precision = 0, bbox = structure(c(xmin = 664069.895400484, 
        ymin = 6524290.81442892, xmax = 665509.674854913, ymax = 6526138.02793883
        ), class = "bbox"), crs = structure(list(epsg = 32633L, proj4string = "+proj=utm +zone=33 +datum=WGS84 +units=m +no_defs"), class = "crs"), n_empty = 0L)), row.names = 139885L, class = c("sf", 
        "data.frame"), sf_column = "geometry", agr = structure(c(id = NA_integer_, 
        osm_id = NA_integer_), .Label = c("constant", "aggregate", "identity"
        ), class = "factor"))

Answer 1

I don't know why st_snap returns the start/end point of the linestring. Maybe that is something for an issue at https://github.com/r-spatial/sf.

Here's a workaround that seems to identify the correct point. Note that st_nearest_points was only introduced recently, so you may need to install sf from github.

nrst = st_nearest_points(point1, roads)
new_point2 = st_cast(nrst, "POINT")[2]

mapView(point1, color="red") + st_buffer(point1, dist = cut_dist) + new_point2 + roads

Wrapping this is a function to return a new geometry set with snapped points below a certain max_dist:

library(sf)
library(mapview)

st_snap_points = function(x, y, max_dist = 1000) {

  if (inherits(x, "sf")) n = nrow(x)
  if (inherits(x, "sfc")) n = length(x)

  out = do.call(c,
                lapply(seq(n), function(i) {
                  nrst = st_nearest_points(st_geometry(x)[i], y)
                  nrst_len = st_length(nrst)
                  nrst_mn = which.min(nrst_len)
                  if (as.vector(nrst_len[nrst_mn]) > max_dist) return(st_geometry(x)[i])
                  return(st_cast(nrst[nrst_mn], "POINT")[2])
                })
  )
  return(out)
}

brew = st_transform(breweries, st_crs(trails))

tst = st_snap_points(brew, trails, 500)

mapview(breweries) + mapview(tst, col.regions = "red") + trails