Find the intersection of two geodesics on the WGS84 ellipsoid. Several methods are available:
geodesic_intersect() - Find the closest intersection of two geodesics
defined by starting points and azimuths
geodesic_intersect_segment() - Find the intersection of two geodesic
segments defined by their endpoints
geodesic_intersect_next() - Find the next closest intersection from a
known intersection point
geodesic_intersect_all() - Find all intersections within a maximum
distance
geodesic_intersect(x, azi_x, y, azi_y)
geodesic_intersect_segment(x1, x2, y1, y2)
geodesic_intersect_next(x, azi_x, azi_y)
geodesic_intersect_all(x, azi_x, y, azi_y, maxdist)Coordinates for geodesic X: a vector of c(lon, lat), a matrix
with columns [lon, lat], or a list with lon and lat components.
Azimuth(s) for geodesic X in degrees.
Coordinates for geodesic Y (same format as x).
Azimuth(s) for geodesic Y in degrees.
Start and end coordinates for segment X.
Start and end coordinates for segment Y.
Maximum distance (in meters) for finding all intersections.
A data frame with columns:
x - Displacement along geodesic X from its starting point (meters)
y - Displacement along geodesic Y from its starting point (meters)
coincidence - Indicator: 0 = normal intersection, +1 = geodesics are
parallel and coincident, -1 = geodesics are antiparallel and coincident
lat - Latitude of intersection point (degrees)
lon - Longitude of intersection point (degrees)
For geodesic_intersect_segment(), an additional column segmode indicates
whether the intersection lies within both segments (0), or which segment(s)
the intersection lies outside of.
For geodesic_intersect_all(), returns a list of data frames (one per input
pair of geodesics).
The intersection point is found using the algorithm described in:
C. F. F. Karney, "Geodesic intersections", J. Surveying Eng. 150(3), 04024005:1-9 (2024). doi:10.1061/JSUED2.SUENG-1483
The "closest" intersection minimizes the L1 distance, defined as
|x| + |y| where x and y are the displacements along the two geodesics.
For segment intersection, segmode encodes whether the intersection lies
within the segments:
segmode = 0 means the intersection lies within both segments
Non-zero values indicate the intersection lies outside one or both segments
The coincidence indicator is useful for detecting when geodesics are parallel or antiparallel at their intersection.
geodesic_inverse() for computing azimuths between points
# Two geodesics from different starting points
# Geodesic X: starts at (0, 0), azimuth 45 degrees
# Geodesic Y: starts at (1, 0), azimuth 315 degrees
geodesic_intersect(c(0, 0), 45, c(1, 0), 315)
#> x y coincidence lat lon
#> 1 78712.76 78712.76 0 0.5033503 0.5
# Vectorized: multiple pairs of geodesics
geodesic_intersect(
cbind(c(0, 0, 0), c(0, 0, 0)),
c(30, 45, 60),
cbind(c(1, 1, 1), c(0, 0, 0)),
c(330, 315, 300)
)
#> x y coincidence lat lon
#> 1 111310.96 111310.96 0 0.8717833 0.5
#> 2 78712.76 78712.76 0 0.5033503 0.5
#> 3 64269.79 64269.79 0 0.2906144 0.5
# Intersection of two geodesic segments
# Segment X: (0, -1) to (0, 1)
# Segment Y: (-1, 0) to (1, 0)
geodesic_intersect_segment(
c(0, -1), c(0, 1),
c(-1, 0), c(1, 0)
)
#> x y segmode coincidence lat lon
#> 1 110574.4 111319.5 0 0 0 0
# Find the next intersection from a known intersection point
# Two geodesics crossing at (0, 0) with azimuths 45 and 315
geodesic_intersect_next(c(0, 0), 45, 315)
#> x y coincidence lat lon
#> 1 -19987140 19987140 0 -4.978254e-15 -179.5734
# Find all intersections within 1,000,000 meters
geodesic_intersect_all(c(0, 0), 45, c(1, 0), 315, maxdist = 1e6)
#> x y coincidence lat lon
#> 1 78712.76 78712.76 0 0.5033503 0.5