lwgeom icon indicating copy to clipboard operation
lwgeom copied to clipboard

Feat/pairwise azimuth

Open robitalec opened this issue 1 year ago • 0 comments

Adds an optional pairwise azimuth implementation for st_geod_azimuth.

See #96.

Reprex (updated 2024-12-06)

# Setup points
library(sf)
#> Linking to GEOS 3.12.2, GDAL 3.9.2, PROJ 9.4.1; sf_use_s2() is TRUE
p <- st_sfc(st_point(c(7,8)), st_point(c(50,20)), crs = 4326)
p2 <- st_sfc(st_point(c(1,10)), st_point(c(1,20)), crs = 4326)

# geosphere bearing with one and two sets of points
library(geosphere)
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, were retired in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
geosphere_seq <- bearing(st_coordinates(p))
geosphere_pairwise <- bearing(p1 = st_coordinates(p), p2 = st_coordinates(p2))

# lwgeom st_geod_azimuth with one and two sets of points
library(lwgeom)
#> Linking to liblwgeom 3.0.0beta1 r16016, GEOS 3.12.2, PROJ 9.4.1
#> 
#> Attaching package: 'lwgeom'
#> The following object is masked from 'package:sf':
#> 
#>     st_perimeter
lwgeom_seq <- st_geod_azimuth(x = p)
lwgeom_pairwise <- st_geod_azimuth(x = p, y = p2)

# Check
library(CircStats)
#> Loading required package: MASS
#> Loading required package: boot
rad(geosphere_seq)
#> [1] 1.209948       NA
lwgeom_seq
#> 1.209948 [rad]

rad(geosphere_pairwise)
#> [1] -1.239384 -1.416121
lwgeom_pairwise
#> Units: [rad]
#> [1] -1.239384 -1.416121

# Test 
p3 <- st_sfc(st_point(c(10,20)), st_point(c(15,25)), st_point(c(20,30)), crs = 4326)
st_geod_azimuth(p, y = p3)
#> Error: length of x and y must match

Created on 2024-12-06 with reprex v2.1.1

robitalec avatar Dec 05 '24 15:12 robitalec