https://software.intel.com/en-us/blogs/2012/11/30/calculating-a-bearing-between-points-in-location-aware-apps

Submitted by John Mechalas (... on Fri, 11/30/2012 - 08:37

Earlier this week I wrote about how to calculate the distance between two points in a location-aware app. Today, I am going to discuss a related topic: how to calculate the bearing between two points.

Like the shortest-distance problem, the bearing between two points on the globe is calculated using the great circle arc that connects them. With the exception of lines of latitude and longitude, great circle arcs do not follow a constant direction relative to true north and this means that as you travel along the arc your heading will vary.

This is made clear in the figure below, which is a gnomonic projection of the earth, showing our route from Portland to London (the gnomonic projection has a very special property: straight lines on the map correspond to great circle arcs). As you can see, the direction of travel changes along the path. The initial bearing, or forward azimuth, is about 33.6° but the final bearing as we approach London is about 141.5°.

As you travel along a great circle route your bearing to your destination changes. The dotted lines represent the direction of true north relative to the starting and ending points.

To calculate the initial bearing bearing we use the following formula. Note the use of the two-argument form of the arctangent, atan2(y,x), which ensures that the resulting angle is in the correct quadrant:

Θ = atan2( sin(Δλ) * cos(Φ₂), cos(Φ₁) * sin (Φ₂) * cos(Δλ) )

This function will return the angle in radians from -π to π but what we want is an angle in degrees from 0 to 360. To accomplish this, we convert to degrees, add 360, and take the modulo 360:

Θ_d = ( Θ * 180 / π + 360 ) % 360

To get the final bearing, you reverse the latitudes and longitudes, and then take the angle that is in the opposite direction (180 degrees around).

Unlike our great circle distance calculation, the bearing calculation makes use of atan and it contains a singularity: when the two points converge, the angle becomes undefined. This makes perfect sense in the physical world, as if the source and the destination are exactly the same then there is no bearing between them. In practice, rounding errors would probably prevent a perfect equality from occurring, but it would still be good form to assume the points are coincident if their distance is below a threshold distance of a meter or two.

Code

Below are some code snippets that can be used to calculate the bearing between two points. You pass the latitude and longitude (in decimal degrees) for the first point as lat1 and long1, and for the second point in lat2 and long2.

For Windows developers, here is an implementation in C#:

class GreatCircleBearing
        {
            static Double degToRad = Math.PI / 180.0;
     
            static public Double initial (Double lat1, Double long1, Double lat2, Double long2)
            {
                return (_bearing(lat1, long1, lat2, long2) + 360.0) % 360;
            }
     
            static public Double final(Double lat1, Double long1, Double lat2, Double long2)
            {
                return (_bearing(lat2, long2, lat1, long1) + 180.0) % 360;
            }
     
            static private Double _bearing(Double lat1, Double long1, Double lat2, Double long2)
            {
                Double phi1 = lat1 * degToRad;
                Double phi2 = lat2 * degToRad;
                Double lam1 = long1 * degToRad;
                Double lam2 = long2 * degToRad;
     
                return Math.Atan2(Math.Sin(lam2-lam1)*Math.Cos(phi2),
                    Math.Cos(phi1)*Math.Sin(phi2) - Math.Sin(phi1)*Math.Cos(phi2)*Math.Cos(lam2-lam1)
                ) * 180/Math.PI;
            }
        }

And in Javascript:

function bearingInitial (lat1, long1, lat2, long2)
{
    return (bearingDegrees(lat1, long1, lat2, long2) + 360) % 360;
}

function bearingFinal(lat1, long1, lat2, long2) {
    return (bearingDegrees(lat2, long2, lat1, long1) + 180) % 360;
}

function bearingDegrees (lat1, long1, lat2, long2)
{
    var degToRad= Math.PI/180.0;

    var phi1= lat1 * degToRad;
    var phi2= lat2 * degToRad;
    var lam1= long1 * degToRad;
    var lam2= long2 * degToRad;

    return Math.atan2(Math.sin(lam2-lam1) * Math.cos(phi2),
        Math.cos(phi1)*Math.sin(phi2) - Math.sin(phi1)*Math.cos(phi2)*Math.cos(lam2-lam1)
    ) * 180/Math.PI;
}

And for Android developers, an implementation in Java:

    class GreatCircleBearing
    {
        static public double initial (double lat1, double long1, double lat2, double long2)
        {
            return (_bearing(lat1, long1, lat2, long2) + 360.0) % 360;
        }

        static public double final(double lat1, double long1, double lat2, double long2)
        {
            return (_bearing(lat2, long2, lat1, long1) + 180.0) % 360;
        }

        static private double _bearing(double lat1, double long1, double lat2, double long2)
        {
            static double degToRad = Math.PI / 180.0;
            double phi1 = lat1 * degToRad;
            double phi2 = lat2 * degToRad;
            double lam1 = long1 * degToRad;
            double lam2 = long2 * degToRad;

            return Math.atan2(Math.sin(lam2-lam1)*Math.cos(phi2),
                Math.cos(phi1)*Math.sin(phi2) - Math.sin(phi1)*Math.cos(phi2)*Math.cos(lam2-lam1)
            ) * 180/Math.PI;
        }
    }

As with our distance calculations, the assumption behind these formulas is a spherical earth. This is sufficiently accurate for casual use but scientific applications will need a more sophisticated model.