---
title: "Geodesic Calculations"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Geodesic Calculations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(geographiclib)
```

## Contents

- [Example Locations](#example-locations)
- [The Inverse Problem: Distance and Bearing](#the-inverse-problem-distance-and-bearing)
- [The Direct Problem: Finding Destinations](#the-direct-problem-finding-destinations)
- [Geodesic Paths](#geodesic-paths)
- [Distance Matrices](#distance-matrices)
- [Rhumb Lines (Loxodromes)](#rhumb-lines-loxodromes)
- [Polygon Area](#polygon-area)
- [Exact vs Fast Geodesic](#exact-vs-fast-geodesic)
- [Geodesic Intersections](#geodesic-intersections)
- [Nearest Neighbor Search](#nearest-neighbor-search)
- [See Also](#see-also)

## Example Locations

```{r 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

```{r inverse-basic}
# 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

```{r 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)
)
```

### Distances Between Cities

```{r 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)
)
```

### Antarctic Distances

```{r 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)
)
```

## The Direct Problem: Finding Destinations

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

### Basic Direct Calculation

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

### Creating a "Circle" of Destinations

```{r 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)
)
```

### Antarctic Traverse Planning

```{r 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)
)
```

## Geodesic Paths

Generate points along the shortest path between two locations.

### Sydney to London Great Circle

```{r 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
```

### Antarctic Circle Path

```{r 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
```

### Trans-Antarctic Path

```{r 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
```

## Distance Matrices

Compute all pairwise distances between sets of points.

### City Distance Matrix

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

### Antarctic Station Distances

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

### Cross-Matrix: Cities to Antarctic Stations

```{r 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 Lines (Loxodromes)

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

### Geodesic vs Rhumb

```{r 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))
)
```

### East-West Travel

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

```{r 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

```{r 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

```{r 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
```

### Multiple Polygons

```{r 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)
)
```

### Winding Direction Matters

```{r 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 Geodesic

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

```{r 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 Comparison

```{r 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
```

## 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:

```{r 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)
)
```

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:

```{r 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")
}
```

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

```{r 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")
```

### Next Intersection

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

```{r 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))
)
```

### All Intersections Within a Distance

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

```{r 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
```

### Practical Example: Great Circle Route Crossings

Determine where two flight routes cross:

```{r 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")
```

### Vectorized Operations

All intersection functions are vectorized for efficiency:

```{r 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)
)
```

## 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:

```{r 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)
)
```

### Multiple Query Points

Search for multiple query points at once:

```{r 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 = "")
  }
}
```

### Radius Search

Find all points within a specified distance:

```{r 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)
  )
}
```

### Self-Search for Clustering

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

```{r 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)
)
```

### Distance Matrix

Create a complete distance matrix between two sets of points:

```{r 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
```

### 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.

```{r 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
```

## See Also

- `vignette("grid-reference-systems")` for MGRS, Geohash, etc.
- `vignette("projections")` for map projections
- `vignette("local-coordinates")` for Local Cartesian and Geocentric
- [Geodesics on an ellipsoid (Wikipedia)](https://en.wikipedia.org/wiki/Geodesics_on_an_ellipsoid)
- [GeographicLib geodesic documentation](https://geographiclib.sourceforge.io/C++/doc/geodesic.html)