5

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)[1], st_bbox(line.sf)[1])
ymin <- min(st_bbox(points.sf)[2], st_bbox(line.sf)[2])
xmax <- max(st_bbox(points.sf)[3], st_bbox(line.sf)[3])
ymax <- max(st_bbox(points.sf)[4], st_bbox(line.sf)[4])

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)

example

# 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
  • 2
    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
6

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)[1], st_bbox(lineWgs84.sf)[1])
ymin2 <- min(st_bbox(pointsWgs84.sf)[2], st_bbox(lineWgs84.sf)[2])
xmax2 <- max(st_bbox(pointsWgs84.sf)[3], st_bbox(lineWgs84.sf)[3])
ymax2 <- max(st_bbox(pointsWgs84.sf)[4], st_bbox(lineWgs84.sf)[4])

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")

plotDist2Line

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")

differences

  • 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
1

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[[2]], add = TRUE, col = "green", pch = 19)

enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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