st_nearest_points do not intersect nearest line feature #790
Description
I have a sites (pt, POINT) and river network (lines, LINESTRING) dataset. Often, a site is located a few meters away from the nearest line feature. For a given site I would like to locate the nearest point along the nearest line that intersects that line using only the sf package (which is a fantastic package BTW!). I have posted this problem on Stack Exchange but have not found a satisfactory and reliable solution since posting. @barryrowlingson and others replied on SE, and @barryrowlingson kindly posted an issue (#788) requesting such a feature on sf that has already been implemented in 0.6.4.
I installed sf_0.6-4 and made use of st_nearest_feature() and st_nearest_points(), but found that the LINESTRING output from st_nearest_points sometimes does not intersect the nearest line. I'm not sure why this is happening and have tried using the st_snap function as well without success. I would appreciate any help in understanding why the nearest point here is still a small distance away from intersecting the nearest line. I've posted a reproducible example below since I'm uncertain if this is a package issue or not; if it isn't then closing the issue would be OK with me.
library(sf)
library(tidyverse)
pt <- structure(list(SiteId = "NAWA_9", EZG_NR = 102382, h1 = 187641,
h2 = 188276, geometry = structure(list(structure(c(605975,
220845), class = c("XY", "POINT", "sfg"))), class = c("sfc_POINT",
"sfc"), precision = 0, bbox = structure(c(605975, 220845,
605975, 220845), .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",
"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_), class = "factor", .Label = c("constant",
"aggregate", "identity"), .Names = c("SiteId", "EZG_NR", "h1",
"h2")))
lines <- structure(list(OBJECTID_1 = c(50285, 50459, 50460, 50893, 51233,
51826), FNODE_ = c(67310, 67495, 67197, 67935, 68287, 68934),
TNODE_ = c(67197, 65833, 67495, 67495, 67935, 67935), LPOLY_ = c(0,
0, 0, 0, 0, 0), RPOLY_ = c(0, 0, 0, 0, 0, 0), LENGTH = c(450.060414716,
1714.92888098, 349.503218869, 381.106076374, 605.268729687,
835.277620554), OBJECTID = c(1555924, 1556346, 1556073, 17889262,
17889263, 17889261), OBJECTVAL = c("Bach", "Bach", "Bach_U",
"Bach", "Bach", "Bach"), NAME = c(NA, "Limpach", NA, "Limpach",
NA, "Limpach"), UNTERIRDIS = c(NA, NA, "unbestimmt", NA,
NA, NA), OBJECTORIG = c("LK25", "LK25", "LK25", "GN25", "GN25",
"GN25"), YEAROFCHAN = c(1998, 2000, 2000, 2005, 2005, 2005
), GEWISSNR = c(12343, 796, 12343, 796, 380438, 796), LAUFNR = c(0,
0, 0, 0, 0, 0), LINST = c("CH", "CH", "CH", "CH", "CH", "CH"
), BACHNR = c(NA_character_, NA_character_, NA_character_,
NA_character_, NA_character_, NA_character_), GWLNR = c("CH0123430000",
"CH0007960000", "CH0123430000", "CH0007960000", "CH3804380000",
"CH0007960000"), Binary = c(1L, 1L, 1L, 1L, 1L, 1L), Shape_Leng = c(450.060414716,
1714.92888098, 349.503218869, 381.106025792, 605.268728778,
835.27768502), ReachId = c(50285L, 50459L, 50460L, 50893L,
51233L, 51826L), SiteId = c("NAWA_9", "NAWA_9", "NAWA_9",
"NAWA_9", "NAWA_9", "NAWA_9"), EZG_NR = c(102382, 102382,
102382, 102382, 102382, 102382), h1 = c(187641, 187641, 187641,
187641, 187641, 187641), h2 = c(188276, 188276, 188276, 188276,
188276, 188276), geometry = structure(list(structure(c(605489.669045818,
605491.875, 605509.375, 605519.375, 605529.375, 605546.875,
605566.875, 605584.375, 605596.875, 605609.375, 605635.625,
220964.545921764, 220965.59375, 220966.90625, 220969.40625,
220969.40625, 220973.09375, 220975.59375, 220976.203125,
220973.703125, 220968.703125, 220950.59375), .Dim = c(11L,
2L), class = c("XY", "LINESTRING", "sfg")), structure(c(605916.125,
605931.625, 605985.375, 606039.1875, 606092.625, 606113.125,
606126.625, 606173.125, 606266.375, 606311.483253449, 220742.09375,
220766.5, 220852.796875, 220939.203125, 221025.09375, 221058.796875,
221073.90625, 221109.09375, 221179.5, 221214.752345892), .Dim = c(10L,
2L), class = c("XY", "LINESTRING", "sfg")), structure(c(605635.625,
605916.125, 220950.59375, 220742.09375), .Dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(605756.541259766, 605767.375,
605800.6875, 605836.3125, 605854.875, 605877.375, 605916.125,
220397.253540039, 220423.09375, 220508.296875, 220594.59375,
220641.5, 220680.59375, 220742.09375), .Dim = c(7L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(605762.325220285, 605756.541259766,
220392.666224025, 220397.253540039), .Dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg")), structure(c(605755.910450333, 605756.541259766,
220395.72591234, 220397.253540039), .Dim = c(2L, 2L), class = c("XY",
"LINESTRING", "sfg"))), n_empty = 0L, 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"), class = c("sfc_LINESTRING",
"sfc"), precision = 0, bbox = structure(c(605489.669045818,
220392.666224025, 606311.483253449, 221214.752345892), .Names = c("xmin",
"ymin", "xmax", "ymax"), class = "bbox"))), .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", "EZG_NR", "h1", "h2", "geometry"), row.names = c("50285",
"50459", "50460", "50893", "51233", "51826"), 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_, 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("OBJECTID_1", "FNODE_",
"TNODE_", "LPOLY_", "RPOLY_", "LENGTH", "OBJECTID", "OBJECTVAL",
"NAME", "UNTERIRDIS", "OBJECTORIG", "YEAROFCHAN", "GEWISSNR",
"LAUFNR", "LINST", "BACHNR", "GWLNR", "Binary", "Shape_Leng",
"ReachId", "SiteId", "EZG_NR", "h1", "h2")))
# Select nearest reach
reach.nearest <- st_nearest_feature(pt, lines)
reach.nearest <- lines[reach.nearest,]
st_distance(pt, reach.nearest)
# Locate nearest point along nearest reach (Returns line)
pt.nearest <- st_nearest_points(pt, reach.nearest)
cat("Line intersects:", st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
st_distance(pt.nearest, reach.nearest)
pt.nearest <- st_cast(pt.nearest, "POINT")[2]
cat("Point (st_cast) intersects:", st_intersects(pt.nearest, reach.nearest, sparse=FALSE)[1,1], "|")
pt.snap <- st_snap(pt.nearest, reach.nearest, tolerance = 1e-9)
cat("Point (st_snap) intersects:", st_intersects(pt.snap, reach.nearest, sparse=FALSE)[1,1], "|\n")
st_distance(pt.snap, reach.nearest)