inst/doc/geodesics.R

In geographiclib: Access to 'GeographicLib'

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

## ----setup--------------------------------------------------------------------
library(geographiclib)

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

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

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

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

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

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

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

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

## ----path-basic---------------------------------------------------------------
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

## ----path-antarctic-----------------------------------------------------------
# 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

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

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

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

## ----matrix-cross-------------------------------------------------------------
# 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-compare------------------------------------------------------------
# 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))
)

## ----rhumb-eastwest-----------------------------------------------------------
# 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-------------------------------------------------------------
# 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

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

## ----polygon-winding----------------------------------------------------------
# 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------------------------------------------------------------
# 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, eval=FALSE--------------------------------------------------
# # 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

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

## ----intersect-segment--------------------------------------------------------
# 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")
}

## ----intersect-segment-parallel-----------------------------------------------
# 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")

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

## ----intersect-all------------------------------------------------------------
# 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

## ----intersect-routes---------------------------------------------------------
# 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")

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

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

## ----nn-multiple--------------------------------------------------------------
# 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 = "")
  }
}

## ----nn-radius----------------------------------------------------------------
# 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)
  )
}

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

## ----nn-distmat---------------------------------------------------------------
# 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

## ----nn-performance, eval=FALSE-----------------------------------------------
# # 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