Find coordinates x distance along linestring

Question

I would like to extract the coordinates of a point at a known distance along the linestring, starting from one end of the linestring.

For example:

library(sf)

path <- st_as_sfc('LINESTRING(10 20, 11 21, 12 21, 13 22)')
start_point <- st_as_sfc('POINT(10 20)')
nodes <- st_as_sfc('MULTIPOINT(10 20, 11 21, 12 21, 13 22)')

plot(st_geometry(path))
plot(st_geometry(nodes), add = T, pch = 16, col = 'grey')
plot(st_geometry(start_point), add = T, pch = 16, col = 'red')

In the example code/image we have a linestring (nodes in grey) and the starting point (in this example the linestring starting coordinates) in red.

Example of desired output:

distance_line <- st_as_sfc('LINESTRING(10 20, 11 21, 11.5 21)')
point_wanted <- st_as_sfc('POINT(11.5 21)')  

plot(st_geometry(distance_line), col = 'green', lwd = 4, add = T)
plot(st_geometry(point_wanted), add = T, pch = 16, col = 'blue')

Ultimately, I'm wanting to extract the coordinates (e.g. using st_coordinates) of the point at distance X along the linestring from the starting point. This feels like a common enough desire, so apologies if I have missed an obvious solution.

The only method I can see is to sample, using sf::st_line_sample, at a high resolution and extract the nearest value. This seems inefficient as I have many thousand linestrings with each only needing one distance coordinates. Ideally, the proposed method would be sf compatible.

Updated with more realistic data

path <- st_as_sf(data.frame(X = c(444618, 444640, 444661), Y = c(216561, 216556, 216550), L1 = 1), coords = c('X', 'Y'), crs = 27700) %>% 
      group_by(L1) %>%
      summarise(do_union = F) %>% 
      st_cast('LINESTRING')

nodes <- st_as_sf(data.frame(X = c(444618, 444640, 444661), Y = c(216561, 216556, 216550), L1 = 1), coords = c('X', 'Y'), crs = 27700) 

Testing method proposed by @agila:

    st_distance(nodes)[1,]
    Units: [m]
    [1]  0.00000 22.56103 44.38468

Testing with point 2 and 3.
pt1 <- path %>% st_startpoint()

desired_distance <- units::set_units(22.56103, "m")
ratio <- desired_distance / st_length(path)
pt2 <- st_linesubstring(path, from = 0, to = ratio) %>% st_endpoint()

desired_distance <- units::set_units(44.38468, "m")
ratio <- desired_distance / st_length(path)
pt3 <- st_linesubstring(path, from = 0, to = ratio) %>% st_endpoint()

    (st_distance(pt1, pt2))
    Units: [m]
         [,1]
[1,] 22.56103
    (st_distance(pt1, pt3))
Units: [m]
         [,1]
[1,] 44.36801

I don't know why the accuracy of this method appears to scale with the distance but this error is acceptable for my task.

Answer 1

I want to propose the following solution. As you can see, it has some drawbacks, but I think it may fix your problem (depending on the required spatial accuracy) First, load some packages

library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(lwgeom)
#> Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.9.0, PROJ 7.2.1

Then, create the linestring object

path <- st_as_sfc('LINESTRING(10 20, 11 21, 12 21, 13 22)', crs = 4326)

Calculate its length

st_length(path)
#> 407726.3 [m]

If we want to estimate the point at distance 200000 [m] from the starting point of the linestring, then we can use st_linesubstring():

desired_distance <- units::set_units(200000, "m")
ratio <- desired_distance / st_length(path)
(pt <- st_linesubstring(path, from = 0, to = ratio) %>% st_endpoint())
#> Warning in st_linesubstring.sfc(path, from = 0, to = ratio): st_linesubstring
#> does not follow a geodesic; you may want to use st_geod_segmentize first
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 11.46373 ymin: 21 xmax: 11.46373 ymax: 21
#> Geodetic CRS:  WGS 84
#> POINT (11.46373 21)

I’m not 100% sure about that warning message (and you may want to wait for other answers or better explanations), but we can fix it by converting the input object to a projected CRS. For example:

path2 <- st_transform(path, 32632)
(pt2 <- st_linesubstring(path2, from = 0, to = ratio) %>% st_endpoint() %>% st_transform(4326))
#> Geometry set for 1 feature 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 11.46165 ymin: 21.00073 xmax: 11.46165 ymax: 21.00073
#> Geodetic CRS:  WGS 84
#> POINT (11.46165 21.00073)

The two points are not identical but quite close. They both lay around midway in the linestring object.

st_distance(pt, pt2)
#> Units: [m]
#>          [,1]
#> [1,] 230.2247

Plot

par(mar = rep(0, 4))
plot(path, reset = FALSE)
plot(pt, add = TRUE, pch = 16, col = "darkgreen", cex = 3)
plot(pt2, add = TRUE, pch = 16, col = "darkred", cex = 2)

^{Created on 2021-06-16 by the reprex package (v2.0.0)}