Haversine Formula: Calculate Distance Between Points

The Haversine formula calculates the great-circle distance between two points on a sphere given their latitude and longitude coordinates, providing distance in kilometers when using appropriate Earth radius values. For WGS84 coordinates, the formula assumes a spherical Earth model, which introduces some error compared to the actual ellipsoidal shape of WGS84, but remains accurate enough for most practical applications with typical precision of about 0.5% or better.

Contents

• Mathematical Foundation
• Step-by-Step Implementation
• WGS84 Coordinate System Considerations
• Accuracy Comparison with Other Methods
• Practical Code Examples
• Common Pitfalls and Best Practices

Mathematical Foundation

The Haversine formula is based on spherical trigonometry and calculates the shortest distance between two points on the surface of a sphere. The fundamental formula is:

a = sin²(Δφ/2) + cos(φ₁) × cos(φ₂) × sin²(Δλ/2)
c = 2 × atan2(√a, √(1−a))
d = R × c

Where:

• φ₁, φ₂ are the latitudes of point 1 and point 2 (in radians)
• λ₁, λ₂ are the longitudes of point 1 and point 2 (in radians)
• Δφ = φ₂ - φ₁ (difference in latitudes)
• Δλ = λ₂ - λ₁ (difference in longitudes)
• R is Earth’s radius (typically 6371 km for kilometers)
• d is the distance between the two points

The haversine function itself is defined as hav(θ) = sin²(θ/2), which is why the formula bears this name.

Step-by-Step Implementation

1. Convert Degrees to Radians

Since most GPS coordinates are provided in degrees, you’ll need to convert them to radians first:

radians = degrees × (π/180)

2. Calculate Differences

Compute the differences in latitude and longitude:

Δφ = φ₂ - φ₁
Δλ = λ₂ - λ₁

3. Apply the Haversine Formula

Calculate the intermediate values:

a = sin²(Δφ/2) + cos(φ₁) × cos(φ₂) × sin²(Δλ/2)
c = 2 × atan2(√a, √(1−a))
distance = R × c

4. Choose Earth’s Radius

For kilometers, use R = 6371 km (mean Earth radius). For higher precision in specific regions, you might use values ranging from 6356.8 km at the poles to 6378.1 km at the equator.

WGS84 Coordinate System Considerations

Understanding the Limitation

The Haversine formula assumes a perfect spherical Earth, while WGS84 models Earth as an ellipsoid. This fundamental difference introduces some error in the calculations.

Radius Selection for WGS84

For minimum error with WGS84 coordinates, you should use:

• A spherical Earth radius that matches the radial distance on the ellipsoid near your GPS points of interest
• For GPS points randomly distributed on Earth, a radius of approximately 6371.0 km provides good balance
• If working with points in a specific region, consider using a radius appropriate for that latitude

Coordinate System Compatibility

The WGS84 coordinates can be directly used with the Haversine formula without conversion, as they represent the same angular measurements. However, be aware that:

• The coordinate system differences between datums (WGS84, NAD83, etc.) are negligible compared to the spherical vs. ellipsoidal error
• The formula expects coordinates in decimal degrees

Accuracy Comparison with Other Methods

Haversine Formula Accuracy

• Typical error: 0.5% or better for most applications
• Maximum error: Can be higher at extreme latitudes or for very long distances
• Strength: Computationally simple and fast
• Weakness: Assumes spherical Earth model

Vincenty’s Formula

• Accuracy: Considered the gold standard for geodetic calculations
• Method: Takes into account Earth’s ellipticity using WGS84 ellipsoid parameters
• Computational complexity: Significantly more complex than Haversine
• Use case: Surveying, high-precision applications where every meter counts

Other Methods Comparison

According to research findings, for most purposes, the Haversine formula is “accurate enough” and significantly faster than more complex alternatives.

Practical Code Examples

Python Implementation

import math

def haversine(lat1, lon1, lat2, lon2):
    """
    Calculate the great circle distance between two points 
    on the earth (specified in decimal degrees)
    """
    # Convert decimal degrees to radians
    lat1, lon1, lat2, lon2 = map(radians, [lat1, lon1, lat2, lon2])
    
    # Haversine formula
    dlat = lat2 - lat1
    dlon = lon2 - lon1
    a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
    c = 2 * asin(sqrt(a))
    
    # Radius of earth in kilometers. Use 3956 for miles
    r = 6371
    return c * r

# Example usage
lyon = (45.7597, 4.8422)
paris = (48.8567, 2.3508)
distance = haversine(*lyon, *paris)
print(f"Distance: {distance:.2f} km")  # Output: Distance: 392.22 km

JavaScript Implementation

function haversineDistance(lat1, lon1, lat2, lon2) {
    const R = 6371; // Earth's radius in kilometers
    const dLat = (lat2 - lat1) * Math.PI / 180;
    const dLon = (lon2 - lon1) * Math.PI / 180;
    const a = 
        Math.sin(dLat/2) * Math.sin(dLat/2) +
        Math.cos(lat1 * Math.PI / 180) * Math.cos(lat2 * Math.PI / 180) * 
        Math.sin(dLon/2) * Math.sin(dLon/2);
    const c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));
    return R * c;
}

// Example usage
const lyon = [45.7597, 4.8422];
const paris = [48.8567, 2.3508];
const distance = haversineDistance(lyon[0], lyon[1], paris[0], paris[1]);
console.log(`Distance: ${distance.toFixed(2)} km`);

SQL Implementation

-- PostgreSQL with PostGIS extension
SELECT ST_Distance_Sphere(
    ST_MakePoint(lon1, lat1)::geography,
    ST_MakePoint(lon2, lat2)::geography
) / 1000 AS distance_km;

-- Alternative without PostGIS
CREATE OR REPLACE FUNCTION haversine_km(lat1 float, lon1 float, lat2 float, lon2 float)
RETURNS float AS $$
DECLARE
    R float := 6371; -- Earth's radius in kilometers
    dlat float;
    dlon float;
    a float;
    c float;
BEGIN
    dlat := radians(lat2 - lat1);
    dlon := radians(lon2 - lon1);
    a := sin(dlat/2) * sin(dlat/2) +
         cos(radians(lat1)) * cos(radians(lat2)) *
         sin(dlon/2) * sin(dlon/2);
    c := 2 * atan2(sqrt(a), sqrt(1-a));
    RETURN R * c;
END;
$$ LANGUAGE plpgsql;

-- Usage
SELECT haversine_km(45.7597, 4.8422, 48.8567, 2.3508) AS distance_km;

Common Pitfalls and Best Practices

Common Errors to Avoid

1. Forgetting to convert degrees to radians - This is the most common mistake that leads to completely wrong results
2. Using inconsistent coordinate formats - Mixing DMS (degrees, minutes, seconds) with decimal degrees
3. Incorrect Earth radius - Using miles instead of kilometers or vice versa
4. Floating-point precision issues - For very small distances, precision errors can become significant

Performance Considerations

• Pre-compute trigonometric values if calculating many distances with the same starting point
• Batch processing - For multiple calculations, vectorized operations can provide significant speedups
• Approximation for short distances - For distances under 1km, you might use simpler methods

Library Recommendations

• Python: Use the haversine package from PyPI for pre-built, tested implementations
• JavaScript: Implement as shown above or use geographic libraries like Turf.js
• Database: Use built-in functions like ST_Distance_Sphere in PostgreSQL with PostGIS
• Java: Consider GeoTools or Apache Commons Math libraries

When to Use Alternatives

Consider Vincenty’s formula or other geodetic methods when:

• You need survey-level accuracy (better than 1 meter)
• Working with very long distances (>1000 km)
• Applications where small errors could have significant consequences
• Scientific research requiring high precision

Sources

1. Haversine formula - Wikipedia
2. Calculate distance between two latitude-longitude points using haversine formula - Movable Type
3. Haversine formula implementation - Rosetta Code
4. Is Haversine or Vincenty’s formula better for calculating distance? - Stack Overflow
5. Distance on a sphere: The Haversine Formula - Esri Community
6. Haversine formula in Firebird SQL - Synaptica
7. Which spatial reference system is used by the Haversine formula? - GIS Stack Exchange
8. Haversine package - PyPI

Conclusion

The Haversine formula provides an efficient and reasonably accurate method for calculating distances between latitude and longitude points in kilometers. While it assumes a spherical Earth model rather than the WGS84 ellipsoid, the typical 0.5% error is acceptable for most applications including navigation, web applications, and general geographic calculations.

For practical implementation, remember to:

• Convert degrees to radians before calculations
• Use Earth’s radius of 6371 km for kilometers
• Consider Vincenty’s formula only when maximum precision is required
• Leverage existing libraries when available for tested, optimized implementations

The formula’s balance of accuracy and computational efficiency makes it the go-to choice for most distance calculations between geographic coordinates in the WGS84 system.