Feature request: compute azimuth between two sets of points

[robitalec]

I am wondering if it would be possible to provide a second set of points to st_geod_azimuth like eg. geosphere::bearing and sf::st_distance. My use case is calculating the azimuth between pairs of points in two sf objects, instead of sequential azimuth in one sf object.

Reprex:

# Setup points
library(sf)
p <- st_sfc(st_point(c(7,52)), st_point(c(8,53)), crs = 4326)
p2 <- st_sfc(st_point(c(1,10)), st_point(c(2,12)), crs = 4326)

# geosphere bearing with one and two sets of points
library(geosphere)
geosphere_seq <- bearing(st_coordinates(p))
geosphere_pairwise <- bearing(st_coordinates(p), st_coordinates(p2))

# lwgeom st_geod_azimuth with one set of points (and example of pairwise interface)
library(lwgeom)
lwgeom_seq <- st_geod_azimuth(p)
# Feature request:
# lwgeom_pairwise <- st_geod_azimuth(p, y = p2)

# Check
library(CircStats)
rad(geosphere_seq)
#> [1] 0.5410385        NA
lwgeom_seq
#> 0.5410385 [rad]

rad(geosphere_pairwise)
#> [1] -2.987452 -2.985470
# Feature request:
# lwgeom_pairwise

# Example interface from {sf}
st_distance(p)
#> Units: [m]
#>          [,1]     [,2]
#> [1,]      0.0 130175.6
#> [2,] 130175.6      0.0
st_distance(p, p2, by_element = TRUE)
#> Units: [m]
#> [1] 4701732 4590227

^{Created on 2024-12-04 with reprex v2.1.0}

Edit: draft changes are in #97

Let me know what you think, thanks!