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?

  • 1
    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. – Spacedman Jul 5 at 21:25

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.

  • 1
    I've opened an issue in sf to implement geos' gNearestPoints: github.com/r-spatial/sf/issues/788 – Spacedman Jul 6 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. – mbcaradima Jul 6 at 11:28
  • 1
    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... – Spacedman Jul 6 at 13:31
  • gNearestPoints can mix geometries, and once you have the points gDistance will get you the distance. – Spacedman Jul 6 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

enter image description here

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

enter image description here

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? – mbcaradima Jul 5 at 17:10
  • 1
    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. – sebdalgarno Jul 5 at 17:24
  • see my edits addressing the issue of splitting the line at nearest point. – sebdalgarno Jul 5 at 19:36
up vote 0 down vote accepted

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.

strangely enough i've encountered the same problem but trying to locate bus stop along a bus route. Had to make extensive use of st_as_sfc() and st_precision() to get st_split() to work.

Your Answer

 

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.