Solve geodesic problems on the WGS84 ellipsoid using exact algorithms. These functions provide high-precision solutions for:

  • Finding destination points given start, azimuth, and distance (direct problem)

  • Finding distance and azimuths between two points (inverse problem)

  • Generating points along geodesic paths

  • Computing distance matrices

geodesic_direct(x, azi, s)

geodesic_inverse(x, y)

geodesic_path(x, y, n = 100L)

geodesic_line(x, azi, distances)

geodesic_distance(x, y)

geodesic_distance_matrix(x, y = NULL)

Arguments

x

A two-column matrix or data frame of starting coordinates (longitude, latitude) in decimal degrees.

azi

Numeric vector of azimuths (bearings) in degrees, measured clockwise from north.

s

Numeric vector of distances in meters.

y

A two-column matrix or data frame of ending coordinates (longitude, latitude) in decimal degrees.

n

Integer number of points to generate along the path (including start and end points).

distances

Numeric vector of distances from the starting point in meters.

Value

  • geodesic_direct(): Data frame with columns:

    • lon1, lat1: Starting coordinates

    • azi1: Starting azimuth (degrees)

    • s12: Distance (meters)

    • lon2, lat2: Destination coordinates

    • azi2: Azimuth at destination (degrees)

    • m12: Reduced length (meters)

    • M12, M21: Geodesic scale factors

    • S12: Area under geodesic (square meters)

  • geodesic_inverse(): Data frame with columns:

    • lon1, lat1: Starting coordinates

    • lon2, lat2: Ending coordinates

    • s12: Distance (meters)

    • azi1: Azimuth at start (degrees)

    • azi2: Azimuth at end (degrees)

    • m12: Reduced length (meters)

    • M12, M21: Geodesic scale factors

    • S12: Area under geodesic (square meters)

  • geodesic_path(): Data frame with columns:

    • lon, lat: Coordinates along the path

    • azi: Azimuth at each point (degrees)

    • s: Distance from start (meters)

  • geodesic_line(): Data frame with columns:

    • lon, lat: Coordinates at specified distances

    • azi: Azimuth at each point (degrees)

    • s: Distance from start (meters)

  • geodesic_distance(): Numeric vector of distances in meters (pairwise).

  • geodesic_distance_matrix(): Matrix of distances in meters.

Details

These functions use the exact geodesic algorithms from GeographicLib, which provide full double-precision accuracy for all points on the WGS84 ellipsoid.

The direct problem finds the destination given a starting point, initial azimuth (bearing), and distance. This is useful for navigation and creating buffers.

The inverse problem finds the shortest path (geodesic) between two points and returns the distance and azimuths at both endpoints.

The azimuth is measured in degrees from north, with positive values clockwise (east) and negative values counter-clockwise (west). The range is -180° to 180° (e.g., 90° = east, -90° = west, 180° or -180° = south).

Examples

# Direct problem: Where do you end up starting from London,
# heading east for 1000 km?
geodesic_direct(c(-0.1, 51.5), azi = 90, s = 1000000)
#>   lon1 lat1 azi1   s12     lon2     lat2     azi2    m12       M12       M21
#> 1 -0.1 51.5   90 1e+06 14.12014 50.62607 101.0838 995914 0.9877522 0.9877514
#>            S12
#> 1 7.838198e+12

# Inverse problem: Distance from London to New York
geodesic_inverse(c(-0.1, 51.5), c(-74, 40.7))
#>   lon1 lat1 lon2 lat2     s12      azi1      azi2     m12       M12       M21
#> 1 -0.1 51.5  -74 40.7 5587820 -71.62462 -128.7635 4900877 0.6407216 0.6404073
#>             S12
#> 1 -4.040644e+13

# Generate a great circle path
path <- geodesic_path(c(-0.1, 51.5), c(-74, 40.7), n = 100)
head(path)
#>          lon      lat       azi         s
#> 1 -0.1000000 51.50000 -71.62462      0.00
#> 2 -0.8740577 51.65738 -72.23107  56442.62
#> 3 -1.6534078 51.80961 -72.84297 112885.24
#> 4 -2.4379364 51.95665 -73.46021 169327.86
#> 5 -3.2275218 52.09842 -74.08265 225770.49
#> 6 -4.0220347 52.23486 -74.71015 282213.11

# Multiple distances along a bearing
geodesic_line(c(-0.1, 51.5), azi = 45, distances = c(100, 500, 1000) * 1000)
#>          lon      lat      azi     s
#> 1  0.9326282 52.13103 45.81169 1e+05
#> 2  5.3673926 54.55656 49.37084 5e+05
#> 3 11.6730192 57.32605 54.59801 1e+06