# 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")Given two points, find the distance and azimuths between them.
# 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+13The 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
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 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
# 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.6Given a starting point, bearing, and distance, find the destination.
# 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
# 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
# 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.09Generate points along the shortest path between two locations.
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
# 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
# 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.0Compute all pairwise distances between sets of points.
# 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
# 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
# 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 1357Rhumb lines maintain constant bearing. They’re longer than geodesics but easier to navigate.
# 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 692For 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 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 latitudeCalculate the area of polygons on the ellipsoid.
# 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
# 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
# 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 signsgeographiclib 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
# 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.04Find where two geodesics cross on the ellipsoid. This is useful for navigation, surveying, and geometric calculations.
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.71The 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
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 0The 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
# 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 180From 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 19987Find 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 5Determine 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 kmAll 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.8Find 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.
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.1Search 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)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 0Search 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.8Create 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 16977The 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.15vignette("grid-reference-systems") for MGRS, Geohash,
etc.vignette("projections") for map projectionsvignette("local-coordinates") for Local Cartesian and
Geocentric