Geodesic Calculations
In geographiclib: Access to 'GeographicLib'

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(geographiclib)

Example Locations
The Inverse Problem: Distance and Bearing
The Direct Problem: Finding Destinations
Geodesic Paths
Distance Matrices
Rhumb Lines (Loxodromes)
Polygon Area
Exact vs Fast Geodesic
Geodesic Intersections
Nearest Neighbor Search
See Also

Example Locations

# Major cities including Southern Hemisphere
cities <- cbind(
  lon = c(151.21, 147.32, -0.13, -74.01, 139.69, -43.17, -68.30, 166.67),
  lat = c(-33.87, -42.88, 51.51, 40.71, 35.69, -22.91, -54.80, -77.85)
)
rownames(cities) <- c("Sydney", "Hobart", "London", "New York", 
                       "Tokyo", "Rio", "Ushuaia", "McMurdo")

# Antarctic stations
antarctic <- cbind(
  lon = c(166.67, 77.85, -68.13, 39.58, 0),
  lat = c(-77.85, -67.60, -67.57, -69.41, -90)
)
rownames(antarctic) <- c("McMurdo", "Davis", "Palmer", "Progress", "South Pole")

The Inverse Problem: Distance and Bearing

Given two points, find the distance and azimuths between them.

Basic Distance Calculation

# Sydney to London
geodesic_inverse(c(151.21, -33.87), c(-0.13, 51.51))

The output includes: - s12: Distance in meters - azi1: Azimuth at start (bearing to depart on) - azi2: Azimuth at end (bearing of arrival) - m12, M12, M21: Geodesic scale factors - S12: Area under the geodesic

Understanding Azimuth

Azimuth is measured in degrees from north (-180° to 180°): - 0° = North - 90° = East - ±180° = South - -90° = West

# Various directions from Sydney
destinations <- rbind(
  "Due North" = c(151.21, -23),     # Same longitude, further north
  "Due East" = c(170, -33.87),      # Same latitude, further east
  "Due South" = c(151.21, -50),     # Same longitude, further south
  "Due West" = c(130, -33.87)       # Same latitude, further west
)

sydney <- c(151.21, -33.87)
results <- geodesic_inverse(
  cbind(rep(sydney[1], 4), rep(sydney[2], 4)),
  destinations
)

data.frame(
  destination = rownames(destinations),
  azimuth = round(results$azi1, 1),
  distance_km = round(results$s12 / 1000)
)

Distances Between Cities

# Distances from Sydney to other cities
sydney_idx <- 1
results <- geodesic_inverse(
  cbind(rep(cities[sydney_idx, 1], nrow(cities) - 1), 
        rep(cities[sydney_idx, 2], nrow(cities) - 1)),
  cities[-sydney_idx, ]
)

data.frame(
  from = "Sydney",
  to = rownames(cities)[-sydney_idx],
  distance_km = round(results$s12 / 1000),
  bearing = round(results$azi1, 1)
)

Antarctic Distances

# Distances from South Pole to Antarctic stations
pole <- c(0, -90)
results <- geodesic_inverse(
  cbind(rep(pole[1], nrow(antarctic) - 1),
        rep(pole[2], nrow(antarctic) - 1)),
  antarctic[-5, ]  # Exclude South Pole itself
)

data.frame(
  station = rownames(antarctic)[-5],
  distance_from_pole_km = round(results$s12 / 1000),
  bearing = round(results$azi1, 1)
)

The Direct Problem: Finding Destinations

Given a starting point, bearing, and distance, find the destination.

Basic Direct Calculation

# Travel 1000 km due north from Sydney
geodesic_direct(c(151.21, -33.87), azi = 0, s = 1000000)

Creating a "Circle" of Destinations

# Points 1000 km from Sydney in all directions
bearings <- seq(0, 330, by = 30)
results <- geodesic_direct(c(151.21, -33.87), azi = bearings, s = 1000000)

data.frame(
  bearing = bearings,
  direction = c("N", "NNE", "ENE", "E", "ESE", "SSE", 
                "S", "SSW", "WSW", "W", "WNW", "NNW"),
  lon = round(results$lon2, 2),
  lat = round(results$lat2, 2)
)

Antarctic Traverse Planning

# Plan traverse from McMurdo heading inland
mcmurdo <- c(166.67, -77.85)

# Head south-southwest at 195° bearing
waypoints <- geodesic_direct(
  mcmurdo, 
  azi = 195, 
  s = c(0, 100000, 200000, 300000, 400000, 500000)
)

data.frame(
  distance_km = waypoints$s12 / 1000,
  lon = round(waypoints$lon2, 2),
  lat = round(waypoints$lat2, 2)
)

Geodesic Paths

Generate points along the shortest path between two locations.

Sydney to London Great Circle

path <- geodesic_path(c(151.21, -33.87), c(-0.13, 51.51), n = 12)
path

# The path crosses into the Northern Hemisphere via Asia

Antarctic Circle Path

# Path around Antarctica at 70°S
start <- c(0, -70)
end <- c(180, -70)  # Go halfway around

# This will follow a geodesic, not a parallel of latitude!
path <- geodesic_path(start, end, n = 10)
path

# Note: lat varies slightly - geodesics don't follow parallels

Trans-Antarctic Path

# McMurdo to Palmer Station (across the continent)
path <- geodesic_path(c(166.67, -77.85), c(-68.13, -67.57), n = 10)
path

Distance Matrices

Compute all pairwise distances between sets of points.

City Distance Matrix

# Distance matrix for all cities (km)
dist_matrix <- geodesic_distance_matrix(cities) / 1000
round(dist_matrix)

Antarctic Station Distances

# Distances between Antarctic stations
dist_matrix <- geodesic_distance_matrix(antarctic) / 1000
round(dist_matrix)

Cross-Matrix: Cities to Antarctic Stations

# Distance from each city to each Antarctic station
dist_matrix <- geodesic_distance_matrix(cities, antarctic) / 1000
rownames(dist_matrix) <- rownames(cities)
colnames(dist_matrix) <- rownames(antarctic)
round(dist_matrix)

Rhumb Lines (Loxodromes)

Rhumb lines maintain constant bearing. They're longer than geodesics but easier to navigate.

Geodesic vs Rhumb

# Sydney to London: geodesic vs rhumb
start <- c(151.21, -33.87)
end <- c(-0.13, 51.51)

geo_result <- geodesic_inverse(start, end)
rhumb_result <- rhumb_inverse(start, end)

data.frame(
  method = c("Geodesic", "Rhumb"),
  distance_km = round(c(geo_result$s12, rhumb_result$s12) / 1000),
  starting_bearing = round(c(geo_result$azi1, rhumb_result$azi12), 1),
  extra_distance_km = c(0, round((rhumb_result$s12 - geo_result$s12) / 1000))
)

East-West Travel

For east-west travel along a parallel, rhumb = geodesic:

# Due east from Sydney (along parallel)
start <- c(151.21, -33.87)
end <- c(170, -33.87)

geo <- geodesic_inverse(start, end)
rhumb <- rhumb_inverse(start, end)

data.frame(
  method = c("Geodesic", "Rhumb"),
  distance_km = round(c(geo$s12, rhumb$s12) / 1000, 3),
  bearing = round(c(geo$azi1, rhumb$azi12), 3)
)
# Nearly identical for E-W travel

Rhumb Path

# Rhumb path vs geodesic path: Sydney to Cape Town
start <- c(151.21, -33.87)
end <- c(18.42, -33.92)

geo_path <- geodesic_path(start, end, n = 6)
rhumb_path_result <- rhumb_path(start, end, n = 6)

data.frame(
  type = rep(c("Geodesic", "Rhumb"), each = 6),
  point = rep(1:6, 2),
  lon = c(geo_path$lon, rhumb_path_result$lon),
  lat = round(c(geo_path$lat, rhumb_path_result$lat), 2)
)
# Note: rhumb line maintains more constant latitude

Polygon Area

Calculate the area of polygons on the ellipsoid.

Antarctic Ice Shelf Area

# Approximate Ross Ice Shelf boundary
ross_ice_shelf <- cbind(
  lon = c(158, 170, -175, -160, -150, -158, -170, 170, 158),
  lat = c(-77, -78, -78.5, -79, -78, -77.5, -77, -77, -77)
)

result <- polygon_area(ross_ice_shelf)
result

# Area in km²
abs(result$area) / 1e6

Multiple Polygons

# Compare areas of different regions
pts <- cbind(
  lon = c(
    # Tasmania (approximate)
    144, 148, 148, 144, 144,
    # South Island NZ (approximate)
    166, 174, 174, 166, 166
  ),
  lat = c(
    -40, -40, -44, -44, -40,
    -41, -41, -47, -47, -41
  )
)

result <- polygon_area(pts, id = c(rep(1, 5), rep(2, 5)))
data.frame(
  region = c("Tasmania", "South Island NZ"),
  area_km2 = round(abs(result$area) / 1e6)
)

Winding Direction Matters

# Counter-clockwise vs clockwise
ccw <- cbind(lon = c(0, 1, 1, 0), lat = c(0, 0, 1, 1))
cw <- cbind(lon = c(0, 0, 1, 1), lat = c(0, 1, 1, 0))

data.frame(
  winding = c("Counter-clockwise", "Clockwise"),
  area_km2 = c(polygon_area(ccw)$area, polygon_area(cw)$area) / 1e6
)
# Areas have opposite signs

Exact vs Fast Geodesic

geographiclib provides both exact and fast (series approximation) versions:

# Compare accuracy and speed
start <- c(151.21, -33.87)
end <- c(-0.13, 51.51)

exact <- geodesic_inverse(start, end)
fast <- geodesic_inverse_fast(start, end)

data.frame(
  version = c("Exact", "Fast"),
  distance_m = c(exact$s12, fast$s12),
  difference_mm = c(0, (fast$s12 - exact$s12) * 1000)
)
# Difference is typically nanometers - negligible for most uses

Performance Comparison

# Generate 10,000 random point pairs
n <- 10000
pts1 <- cbind(runif(n, -180, 180), runif(n, -90, 90))
pts2 <- cbind(runif(n, -180, 180), runif(n, -90, 90))

system.time(geodesic_inverse(pts1, pts2))
#>  user  system elapsed 
#>  0.05    0.00    0.05

system.time(geodesic_inverse_fast(pts1, pts2))
#>  user  system elapsed 
#>  0.04    0.00    0.04

Geodesic Intersections

Find where two geodesics cross on the ellipsoid. This is useful for navigation, surveying, and geometric calculations.

Closest Intersection

Given two geodesics defined by starting points and azimuths, find their closest intersection:

# Two geodesics starting from different points
# Geodesic X: from (0, 0) heading NE at 45°
# Geodesic Y: from (1, 0) heading NW at 315°
result <- geodesic_intersect(c(0, 0), 45, c(1, 0), 315)
result

# The intersection is found ahead on both geodesics
data.frame(
  geodesic = c("X", "Y"),
  start = c("(0, 0)", "(1, 0)"),
  azimuth = c(45, 315),
  displacement_km = round(c(result$x, result$y) / 1000, 2)
)

The output includes: - x, y: Displacements along each geodesic from the starting point (meters) - coincidence: 0 = normal intersection, +1 = parallel/coincident, -1 = antiparallel - lat, lon: Coordinates of the intersection point

Segment Intersection

Find if and where two geodesic segments (defined by endpoints) intersect:

# Two crossing segments
# Segment X: from (0, -0.5) to (0, 0.5) - roughly N-S
# Segment Y: from (-0.5, 0) to (0.5, 0) - roughly E-W
result <- geodesic_intersect_segment(
  c(0, -0.5), c(0, 0.5),   # Segment X endpoints
  c(-0.5, 0), c(0.5, 0)    # Segment Y endpoints
)
result

# segmode = 0 means intersection is within both segments
if (result$segmode == 0) {
  cat("Segments intersect at:", result$lat, result$lon, "\n")
} else {
  cat("Segments do not intersect (extended intersection at:", 
      result$lat, result$lon, ")\n")
}

The segmode value indicates whether the intersection lies within the segments: - segmode = 0: Intersection is within both segments - Non-zero: Intersection is outside one or both segments

Non-Intersecting Segments

# Two parallel segments that don't intersect
result <- geodesic_intersect_segment(
  c(0, 0), c(0, 1),    # Segment X: lon=0, lat 0 to 1
  c(1, 0), c(1, 1)     # Segment Y: lon=1, lat 0 to 1
)

cat("segmode:", result$segmode, "- segments do not intersect\n")
cat("Closest point of intersection would be at:", result$lat, result$lon, "\n")

Next Intersection

From a known intersection point, find the next closest intersection. Geodesics on an ellipsoid can intersect multiple times:

# Two geodesics crossing at the origin
# Find where they cross again
next_int <- geodesic_intersect_next(c(0, 0), 45, 90)

data.frame(
  intersection = c("origin", "next"),
  lat = c(0, next_int$lat),
  lon = c(0, next_int$lon),
  distance_km = c(0, round(abs(next_int$x) / 1000))
)

All Intersections Within a Distance

Find all intersections within a specified distance from the starting points:

# Find all intersections within 5,000 km
all_ints <- geodesic_intersect_all(
  c(0, 0), 30,      # Geodesic X: from origin, azimuth 30°
  c(10, 0), 330,    # Geodesic Y: from (10, 0), azimuth 330°
  maxdist = 5e6     # Search radius: 5,000 km
)

cat("Found", nrow(all_ints), "intersections within 5,000 km\n")
all_ints

Practical Example: Great Circle Route Crossings

Determine where two flight routes cross:

# Route 1: Sydney to London (via typical great circle)
sydney <- c(151.21, -33.87)
london <- c(-0.13, 51.51)
route1 <- geodesic_inverse(sydney, london)

# Route 2: Tokyo to São Paulo
tokyo <- c(139.69, 35.69)
sao_paulo <- c(-46.63, -23.55)
route2 <- geodesic_inverse(tokyo, sao_paulo)

# Find intersection of extended routes
crossing <- geodesic_intersect(
  sydney, route1$azi1,
  tokyo, route2$azi1
)

cat("Routes cross at:", round(crossing$lat, 2), "lat,", 
    round(crossing$lon, 2), "lon\n")

# Distance from each origin to crossing
cat("From Sydney:", round(crossing$x / 1000), "km\n")
cat("From Tokyo:", round(crossing$y / 1000), "km\n")

Vectorized Operations

All intersection functions are vectorized for efficiency:

# Multiple intersection problems at once
results <- geodesic_intersect(
  cbind(c(0, 0, 0), c(0, 10, 20)),  # Three starting points for X
  c(45, 60, 75),                     # Three azimuths for X
  cbind(c(5, 5, 5), c(0, 10, 20)),  # Three starting points for Y
  c(315, 300, 285)                   # Three azimuths for Y
)

data.frame(
  problem = 1:3,
  lat = round(results$lat, 4),
  lon = round(results$lon, 4),
  dist_x_km = round(results$x / 1000, 1),
  dist_y_km = round(results$y / 1000, 1)
)

Nearest Neighbor Search

Find the closest points in a dataset using true geodesic distances on the WGS84 ellipsoid. This is useful for spatial queries, clustering, and matching points across datasets.

Basic k-Nearest Neighbors

Find the k closest points in a dataset for each query point:

# Australian cities dataset
oz_cities <- cbind(
  lon = c(151.21, 144.96, 153.03, 115.86, 138.60, 149.13),
  lat = c(-33.87, -37.81, -27.47, -31.95, -34.93, -35.28)
)
rownames(oz_cities) <- c("Sydney", "Melbourne", "Brisbane", 
                          "Perth", "Adelaide", "Canberra")

# Find 3 nearest cities to Canberra
query <- c(149.13, -35.28)  # Canberra
result <- geodesic_nn(oz_cities, query, k = 3)

# Show nearest cities
data.frame(
  rank = 1:3,
  city = rownames(oz_cities)[result$index[, 1]],
  distance_km = round(result$distance[, 1] / 1000, 1)
)

Multiple Query Points

Search for multiple query points at once:

# New cities to find neighbors for
queries <- cbind(
  lon = c(147.32, 130.84),
  lat = c(-42.88, -12.46)
)
rownames(queries) <- c("Hobart", "Darwin")

result <- geodesic_nn(oz_cities, queries, k = 2)

# Results for each query
for (i in 1:nrow(queries)) {
  cat(rownames(queries)[i], "- nearest cities:\n")
  for (j in 1:2) {
    cat("  ", j, ". ", rownames(oz_cities)[result$index[j, i]], 
        " (", round(result$distance[j, i] / 1000), " km)\n", sep = "")
  }
}

Radius Search

Find all points within a specified distance:

# Find all cities within 500 km of Adelaide
adelaide <- c(138.60, -34.93)
within_500km <- geodesic_nn_radius(oz_cities, adelaide, radius = 500000)

if (nrow(within_500km[[1]]) > 0) {
  data.frame(
    city = rownames(oz_cities)[within_500km[[1]]$index],
    distance_km = round(within_500km[[1]]$distance / 1000, 1)
  )
}

Self-Search for Clustering

Search a dataset against itself to find each point's nearest neighbors (useful for clustering or outlier detection):

# Find each city's nearest neighbor (excluding itself)
result <- geodesic_nn(oz_cities, oz_cities, k = 2)

# The second neighbor (k=2) is the nearest *other* city
data.frame(
  city = rownames(oz_cities),
  nearest = rownames(oz_cities)[result$index[2, ]],
  distance_km = round(result$distance[2, ] / 1000, 1)
)

Distance Matrix

Create a complete distance matrix between two sets of points:

# Distance from each Australian city to key world cities
world_cities <- cities[c("Sydney", "London", "New York", "Tokyo"), ]

# Get distances (k = all world cities)
result <- geodesic_nn(world_cities, oz_cities, k = nrow(world_cities))

# Format as distance matrix (km)
dist_mat <- matrix(
  round(result$distance / 1000),
  nrow = nrow(world_cities),
  dimnames = list(rownames(world_cities), rownames(oz_cities))
)
dist_mat

Performance

The nearest neighbor search uses a vantage-point tree, which provides efficient O(log n) queries after O(n log n) construction. This is much faster than computing all pairwise distances for large datasets.

# Example with larger dataset (not run)
set.seed(42)
large_dataset <- cbind(
  lon = runif(10000, 110, 155),
  lat = runif(10000, -45, -10)
)
queries <- cbind(
  lon = runif(100, 110, 155),
  lat = runif(100, -45, -10)
)

system.time(geodesic_nn(large_dataset, queries, k = 10))
#>  user  system elapsed 
#>  0.15    0.00    0.15

geographiclib Access to 'GeographicLib'

Geodesic Calculations In geographiclib: Access to 'GeographicLib'

Contents