A comprehensive Python program that calculates the shortest path between any two points using two methods:
- Great Circle Distance - Theoretical shortest distance on Earth's surface (as the crow flies)
- Road Routing - Actual driving/cycling/walking routes with turn-by-turn directions
Calculates the theoretical shortest distance between two geographic coordinates using spherical geometry and the Haversine formula. This represents the path an airplane or bird would take.
Provides real-world driving, cycling, or walking directions using actual road networks. Includes turn-by-turn navigation, distance, and estimated travel time.
- Haversine Formula: Accurate spherical distance computation
- Multiple Units: Kilometers, miles, and nautical miles
- Initial Bearing: Compass direction to travel
- Waypoint Generation: Intermediate points along the optimal great circle path
- Address Geocoding: Convert addresses to coordinates automatically
- Multiple Transport Modes: Driving, cycling, or walking
- Turn-by-Turn Directions: Detailed navigation instructions
- Distance & Duration: Actual road distance and estimated travel time
- Real Road Networks: Uses OpenStreetMap data via OSRM
- Comparison Mode: Compare great circle vs actual road distance
- Flexible Input: Accept coordinates or street addresses
- Command-line Interface: Easy-to-use CLI with helpful examples
No external dependencies required! Uses only Python standard library and free public APIs.
# Make the script executable
chmod +x shortest_path.pyBasic distance calculation:
python3 shortest_path.py 40.7128 -74.0060 51.5074 -0.1278With waypoints and miles:
python3 shortest_path.py 40.7128 -74.0060 51.5074 -0.1278 --unit mi --waypoints 3Named points:
python3 shortest_path.py 40.7128 -74.0060 51.5074 -0.1278 --name1 "New York" --name2 "London"Driving directions between addresses:
python3 shortest_path.py --mode road --address "Times Square, New York, NY" "Statue of Liberty, New York, NY"With turn-by-turn directions:
python3 shortest_path.py --mode road --address "Central Park, NYC" "Empire State Building, NYC" --directionsCycling route:
python3 shortest_path.py --mode road --transport cycling --address "Golden Gate Bridge, SF" "Fisherman's Wharf, SF" --directionsWalking route between coordinates:
python3 shortest_path.py --mode road --transport walking 40.7128 -74.0060 40.7589 -73.9851 --directionsCompare straight-line vs road distance:
python3 shortest_path.py --mode both --address "San Francisco, CA" "Los Angeles, CA"Compare for a cycling trip:
python3 shortest_path.py --mode both --transport cycling --address "Seattle, WA" "Portland, OR"Sydney to Los Angeles:
python3 shortest_path.py -33.8688 151.2093 34.0522 -118.2437 --unit nm --waypoints 5--mode {great-circle,road,both}: Choose routing method (default: great-circle)great-circle: Calculate theoretical shortest distance on globeroad: Calculate actual road route with directionsboth: Show both and compare
--address: Interpret inputs as addresses instead of coordinates- Coordinates:
lat1 lon1 lat2 lon2(four numeric values) - Addresses:
--address "Start Address" "End Address"(two quoted strings)
--transport {driving,cycling,walking}: Transport mode for road routing (default: driving)
--unit {km,mi,nm}: Unit of measurement (default: km)--waypoints N: Generate N intermediate waypoints along great circle--name1 NAME: Name for the starting point--name2 NAME: Name for the ending point
--directions: Show detailed turn-by-turn directions
When traveling between two points on Earth, the shortest path is not a straight line on a flat map. Because Earth is approximately spherical, the shortest path follows a great circle — a circle on the sphere's surface whose center coincides with the center of the sphere.
Think of it this way: if you stretch a string tightly between two points on a globe, the string follows a great circle path.
The Haversine formula calculates the great circle distance between two points specified by latitude and longitude coordinates.
Given two points with coordinates (φ₁, λ₁) and (φ₂, λ₂):
a = sin²(Δφ/2) + cos(φ₁) · cos(φ₂) · sin²(Δλ/2)
c = 2 · atan2(√a, √(1−a))
d = R · c
Where:
- φ = latitude (in radians)
- λ = longitude (in radians)
- Δφ = φ₂ - φ₁ (difference in latitude)
- Δλ = λ₂ - λ₁ (difference in longitude)
- R = Earth's radius (6,371 km, 3,959 mi, or 3,440 nm)
- d = distance between the two points
The Haversine formula is preferred over simpler formulas because:
- Numerical stability: Avoids errors when calculating distances between close points
- Accuracy: Works well for both short and long distances
- No singularities: Unlike some formulas, it doesn't break down at the poles or equator
def haversine_distance(point1, point2, unit="km"):
# Convert degrees to radians
lat1_rad = math.radians(point1.lat)
lon1_rad = math.radians(point1.lon)
lat2_rad = math.radians(point2.lat)
lon2_rad = math.radians(point2.lon)
# Calculate differences
dlat = lat2_rad - lat1_rad
dlon = lon2_rad - lon1_rad
# Haversine formula
a = math.sin(dlat / 2) ** 2 + \
math.cos(lat1_rad) * math.cos(lat2_rad) * math.sin(dlon / 2) ** 2
c = 2 * math.asin(math.sqrt(a))
# Multiply by Earth's radius
return radius * cThe formula uses the haversine function: hav(θ) = sin²(θ/2), which is particularly stable for small angles.
The initial bearing (also called forward azimuth) tells you which compass direction to start traveling from point A to reach point B along the great circle.
θ = atan2(sin(Δλ) · cos(φ₂), cos(φ₁) · sin(φ₂) − sin(φ₁) · cos(φ₂) · cos(Δλ))
Where θ is the initial bearing in radians, converted to degrees (0-360°).
The bearing along a great circle path is not constant (except when traveling due north/south or along the equator). This is why航海 navigation requires continuous course corrections or waypoint navigation.
To generate waypoints along the great circle path, we use spherical linear interpolation (slerp).
For a point at fraction f (0 ≤ f ≤ 1) along the path:
A = sin((1−f) · d) / sin(d)
B = sin(f · d) / sin(d)
x = A · cos(φ₁) · cos(λ₁) + B · cos(φ₂) · cos(λ₂)
y = A · cos(φ₁) · sin(λ₁) + B · cos(φ₂) · sin(λ₂)
z = A · sin(φ₁) + B · sin(φ₂)
φᵢ = atan2(z, √(x² + y²))
λᵢ = atan2(y, x)
Where:
- d = angular distance between the two points (in radians)
- f = fraction of the distance (0 = start, 1 = end, 0.5 = midpoint)
- (φᵢ, λᵢ) = coordinates of the intermediate point
- Convert to 3D Cartesian: Transform latitude/longitude into 3D unit vectors on a sphere
- Interpolate: Use weighted combination based on the fraction
- Convert back: Transform the 3D vector back to latitude/longitude
This ensures the waypoints lie exactly on the great circle path, not on a rhumb line or straight line on a flat projection.
Simple linear interpolation of latitude and longitude coordinates would give you points that do not lie on the great circle path. On a 2D map they'd look reasonable, but on a sphere, you'd be traveling a longer route!
Example:
- Linear interpolation between (0°, 0°) and (0°, 180°) might pass through (0°, 90°)
- But the shortest path between these points goes over one of the poles!
- Latitude: -90° (South Pole) to +90° (North Pole)
- Longitude: -180° (West) to +180° (East)
- Positive: North and East
- Negative: South and West
The program treats Earth as a perfect sphere with radius:
- 6,371 km (mean radius)
- 3,959 miles
- 3,440 nautical miles
Note: In reality, Earth is an oblate spheroid (slightly flattened at the poles). For most applications, the spherical model provides sufficient accuracy (within 0.5%). For survey-grade accuracy, use the Vincenty formula with the WGS84 ellipsoid model.
Haversine accuracy:
- Excellent for distances < 10,000 km
- Good for distances up to ~20,000 km (half Earth's circumference)
- Small rounding errors possible for antipodal points (exactly opposite sides of Earth)
Error sources:
- Earth is not a perfect sphere (~0.3% error from ellipsoid shape)
- Terrain elevation not considered
- Floating-point arithmetic (~10⁻¹⁵ relative error)
For most practical purposes, these errors are negligible.
The road routing system consists of three main components:
- Geocoding Service (Nominatim)
- Routing Engine (OSRM)
- Route Processing (Local)
Nominatim is OpenStreetMap's geocoding service that converts addresses to geographic coordinates.
Address String → HTTP Request → Nominatim API → JSON Response → Coordinates
Example:
- Input:
"Empire State Building, New York" - Output:
(40.7484, -73.9857, "Empire State Building, 350 5th Avenue, NYC...")
- API Endpoint:
https://nominatim.openstreetmap.org/search - Rate Limiting: Maximum 1 request per second (enforced by the program)
- No API Key Required: Free to use with proper attribution
- Data Source: OpenStreetMap crowd-sourced data
GET https://nominatim.openstreetmap.org/search?
q=Empire+State+Building,+New+York
&format=json
&limit=1
&addressdetails=1
The program:
- URL-encodes the address
- Adds required
User-Agentheader - Enforces rate limiting with
time.sleep() - Parses JSON to extract latitude/longitude
- Returns display name for user-friendly output
OSRM (Open Source Routing Machine) is a high-performance routing engine using real road network data.
Start Coords + End Coords → OSRM API → Route Geometry + Instructions → Formatted Directions
OSRM supports three transport modes, each with different constraints:
| Profile | Road Types | Speed Assumptions | Turn Restrictions |
|---|---|---|---|
| Car | All drivable roads | Highway speeds | Yes (one-way, no left turn, etc.) |
| Bike | Bike lanes, roads, paths | ~15-20 km/h | Limited |
| Foot | Sidewalks, paths, roads | ~5 km/h | None (can go any direction) |
GET https://router.project-osrm.org/route/v1/{profile}/{lon1},{lat1};{lon2},{lat2}?
steps=true
&geometries=geojson
&overview=full
Important: OSRM uses longitude,latitude order (NOT latitude,longitude)!
{
"code": "Ok",
"routes": [{
"distance": 5234.6, // meters
"duration": 623.4, // seconds
"geometry": {
"coordinates": [[lon1,lat1], [lon2,lat2], ...]
},
"legs": [{
"steps": [{
"distance": 245.3,
"duration": 34.2,
"name": "5th Avenue",
"maneuver": {
"type": "turn",
"modifier": "left",
"location": [lon, lat]
}
}, ...]
}]
}]
}The program processes OSRM's response to create user-friendly output:
The program converts OSRM maneuver types into natural language:
| Maneuver Type | Modifier | Instruction |
|---|---|---|
depart |
north |
"Head north on Main Street" |
turn |
left |
"Turn left onto 5th Avenue" |
turn |
right |
"Turn right onto Broadway" |
merge |
slight left |
"Merge slight left" |
roundabout |
- | "Take roundabout onto Park Ave" |
arrive |
- | "Arrive at destination" |
Smart distance formatting based on magnitude:
- Metric: < 100m → meters; ≥ 100m → kilometers
- Imperial: < 528 feet (0.1 mi) → feet; ≥ 528 feet → miles
- < 60 seconds: "45s"
- 1-60 minutes: "23m 15s"
-
60 minutes: "2h 45m"
User Input:
python3 shortest_path.py --mode road --address "Boston, MA" "New York, NY" --directionsProcessing Steps:
-
Geocode Start
- Query: "Boston, MA"
- Result:
(42.3601, -71.0589, "Boston, Massachusetts, USA")
-
Geocode End
- Query: "New York, NY"
- Result:
(40.7128, -74.0060, "New York, New York, USA")
-
Route Request
- URL:
https://router.project-osrm.org/route/v1/car/-71.0589,42.3601;-74.0060,40.7128?steps=true... - Response: Route with ~50 steps, 346 km, 3.5 hours
- URL:
-
Process & Display
- Parse maneuvers into instructions
- Format distances and durations
- Display summary and turn-by-turn directions
- Free: No API keys, no quotas for reasonable use
- Open Data: Crowd-sourced, constantly updated
- Global Coverage: Addresses worldwide
- Limitation: 1 request/second rate limit
- Free: No authentication required
- Fast: Optimized C++ implementation
- Data: Based on OpenStreetMap road network
- Limitation: Public server for light use only
- Production Alternative: Self-host OSRM server
Road routing uses graph algorithms on a network where:
- Nodes = Intersections, endpoints
- Edges = Road segments with weights (distance, time, restrictions)
OSRM uses optimized variants of Dijkstra's algorithm with preprocessing:
1. Start at source node
2. Explore neighboring nodes, tracking cumulative distance
3. Always expand the closest unvisited node next
4. Continue until reaching destination
5. Backtrack to reconstruct the path
Optimizations used by OSRM:
- Contraction Hierarchies: Preprocess road network for faster queries
- Bidirectional Search: Search from both start and end simultaneously
- Heuristics: A* algorithm with geographic distance estimation
- Depends on address specificity
- City-level: ~center of city (could be kilometers off)
- Street address: Usually within 10-50 meters
- Specific POI: Usually very accurate
- Distance: Typically accurate within 1-2%
- Duration: Estimates based on speed limits and road types
- Does NOT account for: Traffic, weather, construction, driver behavior
- Real-world times can vary ±30% or more
- Live traffic conditions
- Road closures and construction
- Traffic signals and stop signs
- Weather conditions
- Driver skill and vehicle type
- Parking and final approach
- Aviation: Flight path planning and fuel calculation
- Maritime: Ship navigation and route optimization
- Telecommunications: Satellite coverage and signal path calculation
- Research: Studying global phenomena (migration patterns, climate, etc.)
- Logistics: International shipping distance estimates
- Trip Planning: Vacation route planning with distance and time estimates
- Delivery Services: Route optimization for last-mile delivery
- Emergency Response: Calculating actual travel time for ambulances, fire trucks
- Urban Planning: Analyzing accessibility and commute times
- Real Estate: Commute distance calculations from properties
- Cycling & Hiking: Planning routes for outdoor activities
- Logistics: Local and regional delivery route planning
- Ride Sharing: Estimating pickup times and routes
In differential geometry, a geodesic is the shortest path between two points on a curved surface. On a flat plane, geodesics are straight lines. On a sphere, geodesics are segments of great circles.
Consider any path on a sphere between two points. If the path deviates from the great circle, you can "straighten" it by moving it closer to the great circle, reducing its length. The only path that cannot be shortened further is the great circle itself.
A rhumb line (or loxodrome) is a path of constant bearing — easier to navigate but longer than the great circle. Ships and aircraft often use rhumb lines for short distances because they're simpler to follow, but switch to great circle routes (broken into waypoints) for long distances.
- Time Complexity: O(1) for distance calculation, O(n) for n waypoints
- Space Complexity: O(n) for storing n waypoints
- Typical Execution Time: < 1ms for single calculation
This is free and open-source software.
- Haversine Formula - Wikipedia
- Great Circle - Wikipedia
- Aviation Formulary - Ed Williams
- Movable Type Scripts - Geodesy Functions
- OSRM Documentation
- Nominatim API Documentation
- OpenStreetMap Wiki
- Dijkstra's Algorithm - Wikipedia
- Contraction Hierarchies - Wikipedia
Contributions are welcome! Some ideas for enhancements:
- Add support for the more accurate Vincenty formula
- Implement rhumb line calculations for comparison
- Add antipodal point detection and handling
- Include elevation data for actual surface distance
- Support for multiple waypoints (A → B → C → D routing)
- Alternative route suggestions
- Traffic data integration
- Avoid tolls/highways options
- Export routes to GPX/KML format
- Route optimization (traveling salesman problem)
- Batch geocoding from CSV files
- Add interactive map visualization with folium or leaflet
- Create a web interface with Flask/FastAPI
- Generate route elevation profiles
- Export routes to Google Maps / Apple Maps
- Mobile-friendly responsive design
- Caching for geocoding results
- Offline mode with local OSM data
- Support for custom OSRM server endpoints
- Parallel route calculation for comparison
- Integration with public transit APIs
For questions or issues, please open an issue on the GitHub repository.