Equivalent function to rgeos::gProject() in sf

[jfbourdon]

gProject() returns the distance along geometry to point nearest the specified point and as far I know, there is no direct equivalent in sf. The function sf::st_nearest_points() however does part of the job, by returning the point on the line nearest to the specified point. The issue is that getting the distance like gProject involves a lot of R code because we cannot split a line with a point due to floating point precision issues. Thus, does a native equivalent is planned?

Below is how I am currently replicating gProject():

# Create example `line` and `point`
line <- sf::st_linestring(matrix(c(1,2,10,6,12,20), ncol = 2)) |> sf::st_sfc()
point <- sf::st_point(c(9,14)) |> sf::st_sfc()


# Get distance along the line to the point nearest to the specified point using `rgeos`
rgeos::gProject(sf::as_Spatial(line), sf::as_Spatial(point))  # 12.44672


# Get distance along the line to the point nearest to the specified point using `sf`
# https://stackoverflow.com/a/70266714/8512568
st_ends_heading <- function(line)
{
  M <- sf::st_coordinates(line)
  i <- c(2, nrow(M) - 1)
  j <- c(1, -1)
  
  headings <- mapply(i, j, FUN = function(i, j) {
    Ax <- M[i-j,1]
    Ay <- M[i-j,2]
    Bx <- M[i,1]
    By <- M[i,2]
    unname(atan2(Ay-By, Ax-Bx))
  })
  
  return(headings)
}

st_extend_line <- function(line, distance, end = "BOTH")
{
  if (!(end %in% c("BOTH", "HEAD", "TAIL")) | length(end) != 1) stop("'end' must be 'BOTH', 'HEAD' or 'TAIL'")
  
  M <- sf::st_coordinates(line)[,1:2]
  keep <- !(end == c("TAIL", "HEAD"))
  
  ends <- c(1, nrow(M))[keep]
  headings <- st_ends_heading(line)[keep]
  distances <- if (length(distance) == 1) rep(distance, 2) else rev(distance[1:2])
  
  M[ends,] <- M[ends,] + distances[keep] * c(cos(headings), sin(headings))
  newline <- sf::st_linestring(M)
  
  # If input is sfc_LINESTRING and not sfg_LINESTRING
  if (is.list(line)) newline <- sf::st_sfc(newline, crs = sf::st_crs(line))
  
  return(newline)
}

st_project <- function(line, point)
{
  bridge <- sf::st_nearest_points(line, point)
  blade <- st_extend_line(bridge, 0.01)
  split <- lwgeom::st_split(line, blade) |> sf::st_collection_extract("LINESTRING")
  dist_from_head <- sf::st_length(split[1]) |> as.numeric()
  
  return(dist_from_head)
}


st_project(line, point)  # 12.44672

[edzer]

does a native equivalent is planned?

No, but I'd be happy to consider a PR (probably under a different name).

[rsbivand]

For guidance, see: https://r-forge.r-project.org/scm/viewvc.php/pkg/R/rgeos_linearref.R?root=rgeos and https://r-forge.r-project.org/scm/viewvc.php/pkg/src/rgeos_linearref.c?root=rgeos.

However, see https://paleolimbot.github.io/geos/reference/geos_project.html for:

> geos::geos_project(sf::st_as_text(line), sf::st_as_text(point))
[1] 12.44672

or equivalently:

> geos::geos_project(geos::as_geos_geometry(line), geos::as_geos_geometry(point))
[1] 12.44672

[edzer]

Thanks, that looks relatively trivial to implement.

[edzer]

So this is implemented in a branch (for now):

library(sf)
# Linking to GEOS 3.11.1, GDAL 3.6.4, PROJ 9.1.1; sf_use_s2() is TRUE
st_interpolate_line(st_as_sfc("LINESTRING (0 0, 1 1)"), 1)
# Geometry set for 1 feature 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 0.7071068 ymin: 0.7071068 xmax: 0.7071068 ymax: 0.7071068
# CRS:           NA
# POINT (0.7071068 0.7071068)
st_interpolate_line(st_as_sfc("LINESTRING (0 0, 1 1)"), 1, TRUE)
# Geometry set for 1 feature 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 1 ymin: 1 xmax: 1 ymax: 1
# CRS:           NA
# POINT (1 1)
st_project_point(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc("POINT (5 5)"))
# [1] 7.071068
st_project_point(st_as_sfc("LINESTRING (0 0, 10 10)"), st_as_sfc("POINT (5 5)"), TRUE)
# [1] 0.5

[edzer]

Maybe we want to follow PostGIS' naming here? https://postgis.net/docs/ST_LineInterpolatePoint.html and https://postgis.net/docs/ST_LineLocatePoint.html

[edzer]

Since we already have st_line_merge and st_line_sample, I suggest we go for st_line_interpolate and st_line_locate for interpolate and project, respectively.

[jfbourdon]

Thanks! It will be of great help!
When I saw that C++ would be involved, I set aside my hope of making a PR myself.

[edzer]

I settled on st_line_interpolate and st_line_project