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