geodesic_intersect: Geodesic intersections

In geographiclib: Access to 'GeographicLib'

Geodesic intersections

Description

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

Usage

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)

Arguments

x	Coordinates for geodesic X: a vector of c(lon, lat), a matrix with columns ⁠[lon, lat]⁠, or a list with lon and lat components.
azi_x	Azimuth(s) for geodesic X in degrees.
y	Coordinates for geodesic Y (same format as x).
azi_y	Azimuth(s) for geodesic Y in degrees.
x1, x2	Start and end coordinates for segment X.
y1, y2	Start and end coordinates for segment Y.
maxdist	Maximum distance (in meters) for finding all intersections.

Details

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). \Sexpr[results=rd]{tools:::Rd_expr_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.

Value

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).

See Also

geodesic_inverse() for computing azimuths between points

Examples

# 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)

# 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)
)
# 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)
)
# 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)
# Find all intersections within 1,000,000 meters
geodesic_intersect_all(c(0, 0), 45, c(1, 0), 315, maxdist = 1e6)

geographiclib documentation built on March 4, 2026, 9:07 a.m.