This is a similar question to my recent questions on how to find the nearest point along a polyline and how to split a polyline by an intersecting point using only the sf package.

My question has mainly to do with the task in the first question I asked but addresses an unresolved problem: given 378 or more points, find the nearest points that intersect a line. I know what functions to use in sf, the problem is that many of the nearest points are still some small distance (e.g., 10^-11 units away) from the line.

Using sf_0.6-4 the st_nearest_feature and st_nearest_points functions will locate the nearest point along a polyline. However, due to the high precision of the coordinates the nearest point is still a very small distance from the reach and does not intersect it. According to the sf author setting the precision may help, but I've had no success in my own attempts to do intersects using a precision of 1000 (three digits following the decimal).

I've also tried using st_snap with a tolerance of 1e-09 but it does not consistently snap the point; in fact this tool should be unnecessary if the precision is used correctly (which evidently I am not doing). Regardless of messing around with precision settings or snapping the points, the points do not consistently intersect the line.

Does anyone know how this can be done? I've included all sites and nearby rivers and a script below demonstrating my basic approach with sf, along with examples of how I've used st_set_precision and st_snap. I've uploaded the input data here.

library(sf) # requires version 0.6.4

pts <- st_read("./inputs", "sites", stringsAsFactors = F)
rivers <- st_read("./inputs", "rivers_buffers", stringsAsFactors = F)

pts.near <- lapply(pts$SiteId, function(i){
  # Select site, site buffer, and reaches within the site buffer
  pt <- filter(pts, SiteId==i) # get individual site
  pt.buffer <- filter(pt.buffers, SiteId==i)

  rivers.buffer <- st_intersection(rivers, pt.buffer)

  # # Set precision
  # p <- 1000
  # pt <- st_set_precision(pt, precision=p)
  # rivers.buffer <- st_set_precision(rivers.buffer, precision=p)

  # Select nearest reach
  reach.nearest <- st_nearest_feature(pt, rivers.buffer)
  reach.nearest <- rivers.buffer[reach.nearest,]

  # # Locate nearest point along nearest reach (Returns line)
  pt.nearest <- st_nearest_points(pt, reach.nearest)
  pt.nearest <- st_cast(pt.nearest, "POINT")[2]

  # Snap point
  pt.nearest <- st_snap(pt.nearest, reach.nearest, tolerance=1e-09)

  pt.dt <- data.frame(SiteId = i, 
                      ReachId = reach.nearest$FID,
                      Intersects = st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1],
                      Distance = st_distance(pt.nearest, reach.nearest),
                      X = st_coordinates(pt.nearest)[1],
                      Y = st_coordinates(pt.nearest)[2]

result <- do.call(rbind, pts.near)

PS; if I can't solve this in R I will commit heresy and use ArcGIS.

  • Instead of snapping points to the line, why don't you select points that you know will be on the line by using st_line_sample? gis.stackexchange.com/questions/251101/…
    – jsta
    Commented Jul 10, 2018 at 14:43
  • 1
    Three reasons: (1) sampling along the line gives only an approximation of the true nearest point (and the error would be affected by the sample size and length of the line), (2) what I want to do is a fairly straightforward task in a GIS (e.g., ArcGIS) so I'd like to do it in R, and (3) I can always extend my current approach and buffer the nearest point by ~1mm. I can then intersect it with the line (I want to eventually split each line by an intersecting point).
    – bcaradima
    Commented Jul 10, 2018 at 14:56


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.