Calculating distance of point on line from the start of line

Question

I have a line and points placed on this line. I am looking for a solution that allows to determine the distance of each point from the beginning of the line using R (preferably sf packet).

This problem and its very good solution for QGIS is described at Points layer distance from the start of line layer in QGIS I use this solutions and works well. However, I cannot implement it in R. I cannot use the qgisprocess package for this purpose https://github.com/paleolimbot/qgisprocess This package has an algorithm for generating geometry using QGIS expression ("native:geometrybyexpression"), but the above-mentioned the QGIS expression does not generate new geometry, it only updates the attribute table.

I have added the source code below with an example line and points for calculating the distance from the beginning of the line.

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
linestring_matrix <-  rbind(c(0, 0), c(100, 0), c(200, 100), c(400, 200), c(500, 0))
linia <- st_sfc(st_linestring(linestring_matrix))
punkty <- st_sfc(st_point(c(50,0)), st_point(c(150,50)), st_point(c(175,75)),
  st_point(c(250,125)), st_point(c(350,175)), st_point( c(475,50)))
plot(linia)
plot(punkty, add=TRUE)

^{Created on 2022-08-09 by the reprex package (v2.0.1)}

I added the expected result in the form of labels at the points (generated in QGIS):

score 1 · Accepted Answer · answered Aug 09 '22 at 15:48

Here's a method. Create lots of points along the line with a known separation, then find the nearest point on the line to your sample points, then you know the distance along the line.

If the line length is L then sampling L points gets you the precision of 1 unit. Setting scale to 2 here gets you precision of half a unit:

scale = 2
pts = st_cast(st_sample(linia,round(st_length(linia)*scale),type="regular"), "POINT")

The nearest of those points to the sample points is:

st_nearest_feature(punkty, pts)
# [1]  100  342  412  595  818 1266

that's an index, but we know each point is regularly sampled at half a distance unit along the line. So divide by scale:

st_nearest_feature(punkty, pts)/2
# [1]  50.0 171.0 206.0 297.5 409.0 633.0

which matches the QGIS labels, pretty much.

Assuming the coordinates are metres, you can get the distance to the nearest cm (takes a bit longer):

> scale = 100
> pts = st_cast(st_sample(linia,round(st_length(linia)*scale),type="regular"), "POINT")
> st_nearest_feature(punkty, pts)
[1]  5001 17072 20607 29733 40913 63274
> st_nearest_feature(punkty, pts)/scale
[1]  50.01 170.72 206.07 297.33 409.13 632.74

Thank you for the idea of solving the problem. I have improved the algorithm for greater accuracy. Spliting a line with st_cast andst_sample does not take into account the trailing, non-integer part of line (eg with a line length of 5.9, each segment will be 1.18 units, which is a bug). There is a QGIS pointsalonglines algorithm that can be used in R, which marks points along the line from the beginning of the line at a defined distance. Corrected code in the answer below (cannot be pasted in the comment). — wacekk, Aug 10 '22 at 13:06
I think if you generate points based on the line as-is, and the line reversed, then work out line distance from each end using the above algorithm, and take like an average or something, you can probably work around that... — Spacedman, Aug 10 '22 at 14:56

score 0 · Answer 2 · answered Aug 10 '22 at 13:16

Here is an improved code to solve my problem (base on Spacedman code):

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(tmap) # for rendering map
library(qgisprocess) # for using QGIS algorithms in R, installed from https://github.com/paleolimbot/qgisprocess
#> (...)
library(tidyverse)
linia <- st_sfc(st_linestring(rbind(c(0, 0), c(100, 0), c(200, 100), c(400, 200), c(500, 0))))
punkty <- st_sfc(st_point(c(50,0)), st_point(c(150,50)), st_point(c(175,75)), st_point(c(250,125)),
                 st_point(c(350,175)), st_point(c(475,50)))
scale = 10
pts = st_as_sf(qgis_run_algorithm("native:pointsalonglines", INPUT = st_sf(linia), DISTANCE = 1/scale))
#> (...)
st_crs(pts) <- NA
nearest_pts <- st_nearest_feature(punkty, pts)
distances <- as.character(round(pts[nearest_pts,]$distance), 1)
punkty2 <- st_sf(tibble(distances), punkty)
tm_shape(linia)+
  tm_lines()+
  tm_shape(punkty2)+
  tm_dots(size = 1)+
  tm_text("distances", ymod = 1)
#> (...)

^{Created on 2022-08-10 by the reprex package (v2.0.1)}

The solution provides greater accuracy, but is slower due to the use of the algorithm from the qgisprocess package.

score -1 · Answer 3 · answered Aug 09 '22 at 14:41

There is a function in the sf package for that: st_distance

linestring_matrix <-  rbind(c(0, 0), c(100, 0), c(200, 100), c(400, 200), c(500, 0))
linia <- st_sfc(st_linestring(linestring_matrix))
punkty <- st_sfc(st_point(c(50,0)), st_point(c(150,50)), st_point(c(175,75)),
                 st_point(c(250,125)), st_point(c(350,175)), st_point( c(475,50)))
plot(linia)
plot(punkty, add=TRUE)
plot(punkty[1], add = T, col = 'red') #first point

minDist <- st_distance(punkty, by_element = F) # distance matrix
nearestTo1 <- which.min(minDist[2:length(punkty),1])
index <- as.numeric(names(nearestTo1))
plot(punkty[index], add = T, col = 'blue') #nearest point

The st_distance function calculates the distance between points (in this case, Euclidean), not a distance along a line (such as a road or river). I added the expected result in the original post. — wacekk, Aug 09 '22 at 15:11

Calculating distance of point on line from the start of line

3 Answers3