# How to calculate distance from POINT to LINESTRING in R using `sf` library and get all POINT features from LINESTRING where distances were calculated?

I'm trying to calculate distances from a set of `POINT` features (1, 2, 3, 4, 5, 6 & 7) to a `LINESTRING` feature using `sf` library in `R`. I could calculate all the distances using the function `st_distance` but I can't get which `POINT` feature from `LINESTRING` feature were used by the function to calculate those distances:

1. First example: if the minimum distance between `POINT` "2" and the `LINESTRING` is the distance between POINT "2" and the POINT "A" in the `LINESTRING`, how can I get POINT "A"?
2. Second example: if the minimum distance between `POINT` "4" and the `LINESTRING` is the distance between POINT "4" and the POINT "X" that is located between POINT "B" and POINT "C" in the `LINESTRING`, how can I get POINT "X"?

Note: I know it's possible to achieve this using `geosphere::dist2Line` function but at the moment of this question `geosphere` package don't have compatibility with simple features objects. So I'm trying to solve this issue using `sf` without using `Spatial` objects from `sp`.

Below is a reproducible code example and a figure of the problem:

``````# Load Libraries ----------------------------------------------------------

library('sf')

# Test data ---------------------------------------------------------------

points.df <- data.frame(
'x' = c(-53.50000, -54.15489, -54.48560, -52.00000, -52.57810, -49.22097, -48.00000),
'y' = c(-38.54859, -39.00000, -38.80000, -38.49485, -38.00000, -38.50000, -37.74859)
)

line.df <- data.frame(
'x' = c(-52.53557, -52.00000, -50.00000, -48.00000, -46.40190),
'y' = c(-41.00000, -40.60742, -40.08149, -40.82503, -39.00000)
)

# Create 'sf' objects -----------------------------------------------------

points.sf <- st_as_sf(points.df, coords = c("x", "y"))

st_crs(points.sf) <- st_crs(4326) # assign crs
points.sf <- st_transform(points.sf, crs = 32721) # transform

line.sf <- st_sf(id = 'L1', st_sfc(st_linestring(as.matrix(line.df), dim = "XY")))
st_crs(line.sf) <- st_crs(4326) # assign crs
line.sf <- st_transform(line.sf, crs = 32721) # transform

# Plots -------------------------------------------------------------------

xmin <- min(st_bbox(points.sf), st_bbox(line.sf))
ymin <- min(st_bbox(points.sf), st_bbox(line.sf))
xmax <- max(st_bbox(points.sf), st_bbox(line.sf))
ymax <- max(st_bbox(points.sf), st_bbox(line.sf))

plot(points.sf, pch = 19, col = "#53A8BD", xlab = "Longitude", ylab = "Latitude",
xlim = c(xmin,xmax), ylim = c(ymin,ymax), graticule = st_crs(4326), axes = TRUE)
plot(line.sf, col = "#C72259", add = TRUE)
text(st_coordinates(points.sf), as.character(1:7), pos = 3)
text(st_coordinates(line.sf), LETTERS[1:5], pos = 1)
`````` ``````# Distances ---------------------------------------------------------------

distances <- st_distance(x = points.sf, y = line.sf)
print(distances)

Units: m
[,1]
[1,] 262727.8
[2,] 256710.2
[3,] 292476.9
[4,] 223153.4
[5,] 291143.6
[6,] 188868.4
[7,] 198670.4
``````
• Can you clarify: you want the distances from some points to the nearest points of a LINESTRING? Because its possible the minimum distance from a point to a linestring is between points. – Spacedman Jun 15 '17 at 20:18
• @Spacedman I need the distances from some points to the nearest LINESTRING: It could be an existent POINT from the LINESTRING or maybe a non existent point of the LINESTRING (eg: middle point between two POINT of a segment). And also, know the points from LINESTRING where the distances were calculated. You can check geosphere::dist2Line R documentation example. – Guzmán Jun 15 '17 at 20:24
• I would try using `st_line_sample` to generate points along your linestring, then calculate distances using `st_distance`. – jsta Aug 6 '17 at 15:53

While `sf` package don't have a built-in function or `geosphere` is not compatible with `sf` objects I would use a wrapper around `geosphere::dist2Line` function: just getting the matrix of coordinates instead using the entire `sf` object.

I also tried @jsta answer based on sampling the line and I compared the differences between both approaches. Since I'm working with a big dataset of points, I prefer not to add more points.

# Wrapping geosphere::dist2Line function approach

``````# Transform sf objects to EPSG:4326 ---------------------------------------

pointsWgs84.sf <- st_transform(points.sf, crs = 4326)
lineWgs84.sf <- st_transform(line.sf, crs = 4326)

# Calculate distances -----------------------------------------------------

dist <- geosphere::dist2Line(p = st_coordinates(pointsWgs84.sf), line = st_coordinates(lineWgs84.sf)[,1:2])

# Create 'sf' object ------------------------------------------------------

dist.sf <- st_as_sf(as.data.frame(dist), coords = c("lon", "lat"))
dist.sf <- st_set_crs(x = dist.sf, value = 4326)

# Plot --------------------------------------------------------------------

xmin2 <- min(st_bbox(pointsWgs84.sf), st_bbox(lineWgs84.sf))
ymin2 <- min(st_bbox(pointsWgs84.sf), st_bbox(lineWgs84.sf))
xmax2 <- max(st_bbox(pointsWgs84.sf), st_bbox(lineWgs84.sf))
ymax2 <- max(st_bbox(pointsWgs84.sf), st_bbox(lineWgs84.sf))

plot(pointsWgs84.sf, pch = 19, col = "#53A8BD", xlab = "Longitude", ylab = "Latitude",
xlim = c(xmin2, xmax2), ylim = c(ymin2, ymax2), graticule = st_crs(4326), axes = TRUE)
plot(lineWgs84.sf, col = "#C72259", add = TRUE)
plot(dist.sf, col = "#6C5593", add = TRUE, pch = 19, cex = 0.75)
text(st_coordinates(pointsWgs84.sf), as.character(1:7), pos = 3, col = "#53A8BD")
text(st_coordinates(lineWgs84.sf), LETTERS[1:5], pos = 1, col = "#C72259")
text(st_coordinates(dist.sf), as.character(1:7), pos = 2, cex = 0.75, col = "#6C5593")
`````` # Sampling approach based on jsta answer

``````# i = which point, n = number of points to add to segment

samplingApproach <- function(i, n) {
line.sf_sample <- st_line_sample(x = line.sf, n = n)
line.sf_sample <- st_cast(line.sf_sample, "POINT")
closest <- line.sf_sample[which.min(st_distance(line.sf_sample, points.sf[i,]))]
closest <- st_transform(closest, crs = 4326)
return(closest)
}

# Differences between approaches with different 'n'
n2 <- lapply(X = seq(100, 6000, 50), FUN = function(x)
st_distance(samplingApproach(i = 2, n = x), dist.sf[2,]))

# Plot differences --------------------------------------------------------

plot(y = n2,
x = seq(100, 6000, 50),
ylim = c(0, 1000),
xlab = "Density (n)",
ylab = "Difference (meters)", type = "p", cex = 0.5, pch = 19,
col = "#6C5593") lines(x = c(0,6000), y = c(0, 0), col = "#EB266D", lwd = 2) lines(smooth.spline(y = n2, x = seq(100, 6000, 50), spar = 0.75), lwd = 2, col = "#272822")
`````` • It seems like `st_coordinates(pointsWgs84.sf)` also needs column subsetting (`st_coordinates(pointsWgs84.sf) [,1:2]`). Also, this approach only works with WGS84 data.. – Ratnanil Jun 18 '18 at 7:20

Try this:

``````line.sf_sample <- st_line_sample(line.sf, 60)
line.sf_sample <- st_cast(line.sf_sample, "POINT")
points.sf <- st_cast(points.sf, "POINT")

closest <- list()
for(i in seq_len(nrow(points.sf))){
closest[[i]] <- line.sf_sample[which.min(
st_distance(line.sf_sample, points.sf[i,]))]
}

plot(points.sf, pch = 19, col = "#53A8BD",
xlab = "Longitude", ylab = "Latitude",
xlim = c(xmin,xmax), ylim = c(ymin,ymax),
graticule = st_crs(4326), axes = TRUE)
plot(line.sf, col = "#C72259", add = TRUE)
text(st_coordinates(points.sf), as.character(1:7), pos = 3)
text(st_coordinates(line.sf), LETTERS[1:5], pos = 1)
plot(closest[], add = TRUE, col = "green", pch = 19)
`````` You can adjust the `density` argument on `st_line_sample` based on your tolerance for error.

• Thanks for your answer! Since I'm working with a big dataset, I prefer not to add more points by sampling and then have to calculate more distances between pairs of points. Also, the exact point (nearest distance) will depend on the number of points added (tolerance). – Guzmán Aug 7 '17 at 19:45