st_line_interpolate fails in WSG84

[mpadge]

Describe the bug

Hi @edzer, just discovered exactly what i need via the st_line_interpolate() function - thanks! Unfortunately, interpolation seems to fail after intermediate call to st_project_point(), in this case from CRS 4326 / WSG84. (Use case is one where inaccuracy from interpolating spherical coordinates is acceptable.) Reprex follows:

To Reproduce

library (sf)
#> Linking to GEOS 3.13.1, GDAL 3.11.3, PROJ 9.6.0; sf_use_s2() is TRUE
packageVersion ("sf")
#> [1] '1.0.21'
sf_use_s2 ()
#> [1] TRUE
sf::sf_extSoftVersion ()
#>           GEOS           GDAL         proj.4 GDAL_with_GEOS     USE_PROJ_H 
#>       "3.13.1"       "3.11.3"        "9.6.0"         "true"         "true" 
#>           PROJ 
#>        "9.6.0"

This works as expected:

l1 <- st_as_sfc("LINESTRING (10.1 50.1, 10.2 50.2)")
d <- as.numeric (st_length (l1))
st_line_interpolate(l1, seq (0, d, length.out = 5))
#> Geometry set for 5 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 10.1 ymin: 50.1 xmax: 10.2 ymax: 50.2
#> CRS:           NA
#> POINT (10.1 50.1)
#> POINT (10.125 50.125)
#> POINT (10.15 50.15)
#> POINT (10.175 50.175)
#> POINT (10.2 50.2)

But same data in 4326 then do not interpolate:

l1 <- st_as_sfc("LINESTRING (10.1 50.1, 10.2 50.2)", crs = 4326)
d <- as.numeric (st_length (l1))
l2 <- st_line_interpolate(l1, seq (0, d, length.out = 5))
#> although coordinates are longitude/latitude, st_project_point assumes that they
#> are planar
l2
#> Geometry set for 5 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 10.1 ymin: 50.1 xmax: 10.2 ymax: 50.2
#> Geodetic CRS:  WGS 84
#> POINT (10.1 50.1)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
duplicated (l2)
#> [1] FALSE FALSE  TRUE  TRUE  TRUE

Only start and end points are returned, with no interpolated, intermediate points at all. Passing "dist" as units object also fails:

l2 <- st_line_interpolate(l1, seq (units::set_units(0, "m"), d, length.out = 5))
#> although coordinates are longitude/latitude, st_project_point assumes that they
#> are planar
l2
#> Geometry set for 5 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 10.1 ymin: 50.1 xmax: 10.2 ymax: 50.2
#> Geodetic CRS:  WGS 84
#> POINT (10.1 50.1)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
#> POINT (10.2 50.2)
duplicated (l2)
#> [1] FALSE FALSE  TRUE  TRUE  TRUE

Re-projecting works fine:

l1 <- st_as_sfc("LINESTRING (10.1 50.1, 10.2 50.2)", crs = 4326) |>
    st_transform ("+proj=merc +a=6378137 +b=6378137")
d <- as.numeric (st_length (l1))
st_line_interpolate(l1, seq (0, d, length.out = 5)) |>
    st_transform (4326)
#> Geometry set for 5 features 
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: 10.1 ymin: 50.1 xmax: 10.2 ymax: 50.2
#> Geodetic CRS:  WGS 84
#> POINT (10.1 50.1)
#> POINT (10.125 50.12502)
#> POINT (10.15 50.15003)
#> POINT (10.175 50.17502)
#> POINT (10.2 50.2)

^{Created on 2025-09-03 with reprex v2.1.1}

Any help on interpolating within 4326 / WSG84 would be appreciated! Thanks

[edzer]

dist values need to be in units of the coordinates; for unprojected lines this is now required and errors in the above case. What (now) works is

library(sf)
# Linking to GEOS 3.12.2, GDAL 3.11.3, PROJ 9.4.1; sf_use_s2() is TRUE
l1 <- st_as_sfc("LINESTRING (10.1 50.1, 10.2 50.2)", crs = 4326)
u = units::set_units(seq(0, sqrt(2)/10, length.out = 5), degrees)
st_line_interpolate(l1, u)
# although coordinates are longitude/latitude, st_project_point assumes that they
# are planar
# Geometry set for 5 features 
# Geometry type: POINT
# Dimension:     XY
# Bounding box:  xmin: 10.1 ymin: 50.1 xmax: 10.2 ymax: 50.2
# Geodetic CRS:  WGS 84
# POINT (10.1 50.1)
# POINT (10.125 50.125)
# POINT (10.15 50.15)
# POINT (10.175 50.175)
# POINT (10.2 50.2)

[edzer]

Thanks for raising this, @mpadge !

[edzer]

Now covering the case of wrong units, as in

library(sf)
# Linking to GEOS 3.12.2, GDAL 3.11.3, PROJ 9.4.1; sf_use_s2() is TRUE
l1 <- st_as_sfc("LINESTRING (10.1 50.1, 10.2 50.2)", crs = 'OGC:CRS84')
d <- st_length (l1)
di = seq(d * 0, d, length.out = 5)
di
# Units: [m]
# [1]     0.000  3301.622  6603.244  9904.866 13206.488
st_line_interpolate(l1, di)
# although coordinates are longitude/latitude, st_project_point assumes that they
# are planar
# Error: cannot convert m into °
# Execution halted

[mpadge]

Thanks @edzer!