Geodesic Calculations

The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

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))
#>     lon1   lat1  lon2  lat2      s12      azi1     azi2     m12        M12
#> 1 151.21 -33.87 -0.13 51.51 16989429 -40.70198 -119.645 2907749 -0.8927693
#>          M21          S12
#> 1 -0.8859132 -5.58329e+13

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)
)
#>   destination azimuth distance_km
#> 1   Due North     0.0        1205
#> 2    Due East    95.3        1736
#> 3   Due South   180.0        1792
#> 4    Due West   -96.0        1959

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)
)
#>     from       to distance_km bearing
#> 1 Sydney   Hobart        1056  -162.4
#> 2 Sydney   London       16989   -40.7
#> 3 Sydney New York       15988    65.7
#> 4 Sydney    Tokyo        7793   -10.0
#> 5 Sydney      Rio       13540   164.5
#> 6 Sydney  Ushuaia        9481   158.4
#> 7 Sydney  McMurdo        4954   175.4

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)
)
#>    station distance_from_pole_km bearing
#> 1  McMurdo                  1357   166.7
#> 2    Davis                  2501    77.8
#> 3   Palmer                  2504   -68.1
#> 4 Progress                  2299    39.6

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)
#>     lon1   lat1 azi1   s12   lon2      lat2 azi2      m12       M12       M21
#> 1 151.21 -33.87    0 1e+06 151.21 -24.84822    0 995893.8 0.9876953 0.9876879
#>   S12
#> 1   0

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)
)
#>    bearing direction    lon    lat
#> 1        0         N 151.21 -24.85
#> 2       30       NNE 156.19 -25.96
#> 3       60       ENE 160.10 -29.04
#> 4       90         E 161.98 -33.40
#> 5      120       ESE 161.08 -37.99
#> 6      150       SSE 157.19 -41.53
#> 7      180         S 151.21 -42.88
#> 8      210       SSW 145.23 -41.53
#> 9      240       WSW 141.34 -37.99
#> 10     270         W 140.44 -33.40
#> 11     300       WNW 142.32 -29.04
#> 12     330       NNW 146.23 -25.96

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)
)
#>   distance_km    lon    lat
#> 1           0 166.67 -77.85
#> 2         100 165.49 -78.71
#> 3         200 164.11 -79.57
#> 4         300 162.49 -80.42
#> 5         400 160.56 -81.26
#> 6         500 158.23 -82.09

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
#>          lon         lat        azi        s
#> 1  151.21000 -33.8700000  -40.70198        0
#> 2  141.44062 -22.9133364  -36.02650  1544494
#> 3  133.16885 -11.4325913  -33.56792  3088987
#> 4  125.54992   0.2719042  -32.82205  4633481
#> 5  117.91613  11.9713350  -33.64160  6177974
#> 6  109.59580  23.4353367  -36.18856  7722468
#> 7   99.73102  34.3579933  -40.98645  9266961
#> 8   87.07053  44.2399329  -49.05480 10811455
#> 9   69.92381  52.1865274  -61.91318 12355948
#> 10  47.28697  56.7435266  -80.42969 13900442
#> 11  21.99232  56.4856288 -101.65444 15444935
#> 12  -0.13000  51.5100000 -119.64500 16989429

# 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
#>    lon       lat azi         s
#> 1    0 -70.00000 180       0.0
#> 2    0 -74.44687 180  496218.9
#> 3    0 -78.89194 180  992437.8
#> 4    0 -83.33570 180 1488656.7
#> 5    0 -87.77866 180 1984875.6
#> 6  180 -87.77866   0 2481094.4
#> 7  180 -83.33570   0 2977313.3
#> 8  180 -78.89194   0 3473532.2
#> 9  180 -74.44687   0 3969751.1
#> 10 180 -70.00000   0 4465970.0

# 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
#>           lon       lat       azi         s
#> 1   166.67000 -77.85000 142.72445       0.0
#> 2   179.21381 -80.36430 130.40195  384322.2
#> 3  -161.23686 -82.14906 111.07320  768644.4
#> 4  -135.26397 -82.65254  85.32054 1152966.7
#> 5  -110.96764 -81.64004  61.24451 1537288.9
#> 6   -94.11385 -79.54283  44.61151 1921611.1
#> 7   -83.41828 -76.87040  34.13761 2305933.3
#> 8   -76.43775 -73.90629  27.38006 2690255.5
#> 9   -71.62526 -70.78443  22.79228 3074577.8
#> 10  -68.13000 -67.57000  19.52390 3458900.0

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)
#>       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]
#> [1,]     0  1056 16989 15988  7793 13540  9481  4954
#> [2,]  1056     0 17397 16612  8735 12641  8667  3992
#> [3,] 16989 17397     0  5585  9582  9253 13361 17018
#> [4,] 15988 16612  5585     0 10873  7733 10597 15069
#> [5,]  7793  8735  9582 10873     0 18561 16994 12727
#> [6,] 13540 12641  9253  7733 18561     0  4110  8657
#> [7,]  9481  8667 13361 10597 16994  4110     0  4818
#> [8,]  4954  3992 17018 15069 12727  8657  4818     0

Antarctic Station Distances

# Distances between Antarctic stations
dist_matrix <- geodesic_distance_matrix(antarctic) / 1000
round(dist_matrix)
#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    0 2804 3459 3292 1357
#> [2,] 2804    0 4775 1552 2501
#> [3,] 3459 4775    0 3848 2504
#> [4,] 3292 1552 3848    0 2299
#> [5,] 1357 2501 2504 2299    0

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)
#>          McMurdo Davis Palmer Progress South Pole
#> Sydney      4954  5872   8285     7304       6253
#> Hobart      3992  4834   7399     6250       5253
#> London     17018 14696  14360    13796      15711
#> New York   15069 16561  12017    15085      14510
#> Tokyo      12727 12557  16046    14058      13953
#> Rio         8657  8875   5256     7357       7467
#> Ushuaia     4818  6142   1423     5061       3927
#> McMurdo        0  2804   3459     3292       1357

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))
)
#>     method distance_km starting_bearing extra_distance_km
#> 1 Geodesic       16989            -40.7                 0
#> 2    Rhumb       17681            -57.7               692

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)
)
#>     method distance_km bearing
#> 1 Geodesic    1736.113  95.269
#> 2    Rhumb    1738.550  90.000
# 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)
)
#>        type point       lon    lat
#> 1  Geodesic     1 151.21000 -33.87
#> 2  Geodesic     2 132.94716 -48.29
#> 3  Geodesic     3 103.69426 -57.87
#> 4  Geodesic     4  66.00286 -57.89
#> 5  Geodesic     5  36.70887 -48.33
#> 6  Geodesic     6  18.42000 -33.92
#> 7     Rhumb     1 151.21000 -33.87
#> 8     Rhumb     2 124.65820 -33.88
#> 9     Rhumb     3  98.10330 -33.89
#> 10    Rhumb     4  71.54530 -33.90
#> 11    Rhumb     5  44.98420 -33.91
#> 12    Rhumb     6  18.42000 -33.92
# 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
#> [1] 147332224442
#> 
#> $perimeter
#> [1] 2532820
#> 
#> $n
#> [1] 9

# Area in km²
abs(result$area) / 1e6
#> [1] 147332.2

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)
)
#>            region area_km2
#> 1        Tasmania   147189
#> 2 South Island NZ   427261

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
)
#>             winding  area_km2
#> 1 Counter-clockwise  12308.78
#> 2         Clockwise -12308.78
# 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)
)
#>   version distance_m difference_mm
#> 1   Exact   16989429   0.00000e+00
#> 2    Fast   16989429  -3.72529e-06
# 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
#>          x        y coincidence       lat lon
#> 1 78712.76 78712.76           0 0.5033503 0.5

# 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)
)
#>   geodesic  start azimuth displacement_km
#> 1        X (0, 0)      45           78.71
#> 2        Y (1, 0)     315           78.71

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
#>          x        y segmode coincidence lat lon
#> 1 55287.15 55659.75       0           0   0   0

# 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")
}
#> Segments intersect at: 0 0

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")
#> segmode: 4 - segments do not intersect
cat("Closest point of intersection would be at:", result$lat, result$lon, "\n")
#> Closest point of intersection would be at: 90 180

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))
)
#>   intersection           lat       lon distance_km
#> 1       origin  0.000000e+00    0.0000           0
#> 2         next -4.978254e-15 -179.5734       19987

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")
#> Found 1 intersections within 5,000 km
all_ints
#>         x       y coincidence      lat lon
#> 1 1104789 1104789           0 8.641145   5

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")
#> Routes cross at: 1.37 lat, 124.84 lon

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

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)
)
#>   problem     lat lon dist_x_km dist_y_km
#> 1       1  2.5144 2.5     393.3     393.3
#> 2       2 11.4144 2.5     315.1     315.1
#> 3       3 20.6130 2.5     269.8     269.8

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)
)
#>   rank      city distance_km
#> 1    1  Canberra         0.0
#> 2    2    Sydney       246.8
#> 3    3 Melbourne       467.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 = "")
  }
}
#> Hobart - nearest cities:
#>   1. Melbourne (598 km)
#>   2. Canberra (858 km)
#> Darwin - nearest cities:
#>   1. Adelaide (2609 km)
#>   2. Perth (2647 km)

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)
  )
}
#>       city distance_km
#> 1 Adelaide           0

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)
)
#>        city   nearest distance_km
#> 1    Sydney  Canberra       246.8
#> 2 Melbourne  Canberra       467.1
#> 3  Brisbane    Sydney       730.6
#> 4     Perth  Adelaide      2135.4
#> 5  Adelaide Melbourne       653.9
#> 6  Canberra    Sydney       246.8

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
#>          Sydney Melbourne Brisbane Perth Adelaide Canberra
#> Sydney        0       714      731  3298     1165      247
#> London     7793      8156     7130  7890     7819     7918
#> New York  15988     16672    15500 14470    16259    16225
#> Tokyo     16989     16898    16522 18700    17099    16977

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

Geodesic Calculations

Contents

Example Locations

The Inverse Problem: Distance and Bearing

Basic Distance Calculation

Understanding Azimuth

Distances Between Cities

Antarctic Distances

The Direct Problem: Finding Destinations

Basic Direct Calculation

Creating a “Circle” of Destinations

Antarctic Traverse Planning

Geodesic Paths

Sydney to London Great Circle

Antarctic Circle Path

Trans-Antarctic Path

Distance Matrices

City Distance Matrix

Antarctic Station Distances

Cross-Matrix: Cities to Antarctic Stations

Rhumb Lines (Loxodromes)

Geodesic vs Rhumb

East-West Travel

Rhumb Path

Polygon Area

Antarctic Ice Shelf Area

Multiple Polygons

Winding Direction Matters

Exact vs Fast Geodesic

Performance Comparison

Geodesic Intersections

Closest Intersection

Segment Intersection

Non-Intersecting Segments

Next Intersection

All Intersections Within a Distance

Practical Example: Great Circle Route Crossings

Vectorized Operations

Nearest Neighbor Search

Basic k-Nearest Neighbors

Multiple Query Points

Radius Search

Self-Search for Clustering

Distance Matrix

Performance

See Also