Calculating bearing between two CLLocationCoordinate2Ds

Question

Very "simple" problem: given two CLLocationCoordinate2Ds, how can I get the bearing (as radians) from the first to the second? I've done a lot of research and studying on this, both the general problem and Objective-C/Cocoa Touch/iOS specifically.

Here's my implementation:

- (float) getHeadingForDirectionFromCoordinate:(CLLocationCoordinate2D)fromLoc toCoordinate:(CLLocationCoordinate2D)toLoc
{
    float fLat = fromLoc.latitude;
    float fLng = fromLoc.longitude;
    float tLat = toLoc.latitude;
    float tLng = toLoc.longitude;

    return atan2(sin(fLng-tLng)*cos(tLat), cos(fLat)*sin(tLat)-sin(fLat)*cos(tLat)*cos(fLng-tLng));         
}

However, this method isn't returning consistant results for me. If the bearing is close to due north or due south, it seems to be fine, however, any other direction seems to return inconsistant data, for example:

From 50.405018, 8.437500

To 51.339802, 12.403340

My method returns: 5.918441 radians

Should be 1.18660576 radians

(see http://www.movable-type.co.uk/scripts/latlong.html and http://www.movable-type.co.uk/scripts/latlong-map.html?lat1=50.405018&long1=8.437500&lat2=51.339802&long2=12.403340)

I've double and triple checked the formula is correct. I've also spot checked a bunch of values like the example above, some correct, some wrong. I've played around with various modulos or bounding of the return value, also no luck.

Any ideas? Is there an issue with my code? Maybe I've misunderstood something about how math functions work?

Answer 1

Here the code modified with the changes suggested by Oren Trutner and from myself:

#define degreesToRadians(x) (M_PI * x / 180.0)
#define radiansToDegrees(x) (x * 180.0 / M_PI)

- (float)getHeadingForDirectionFromCoordinate:(CLLocationCoordinate2D)fromLoc toCoordinate:(CLLocationCoordinate2D)toLoc
{
    float fLat = degreesToRadians(fromLoc.latitude);
    float fLng = degreesToRadians(fromLoc.longitude);
    float tLat = degreesToRadians(toLoc.latitude);
    float tLng = degreesToRadians(toLoc.longitude);

    float degree = radiansToDegrees(atan2(sin(tLng-fLng)*cos(tLat), cos(fLat)*sin(tLat)-sin(fLat)*cos(tLat)*cos(tLng-fLng)));

    if (degree >= 0) {
        return degree;
    } else {
        return 360+degree;
    }
}

Answer 2

Your math is correct, with the following exceptions:

1. Make sure to convert fLat, fLon, tLat, and tLon to radians before applying any sin() or cos() to them. Divide by 180.0 and multiply by PI.

2. Enter the delta between tLng and fLng as tLng-fLng, and not the other way around. Note that this difference appears twice in the expression.

With those changes, I am getting 1.18660677830947 radians with double precision math and the values in the question.

Answer 3

Swift 3:

extension CLLocationCoordinate2D {
    func bearing(to point: CLLocationCoordinate2D) -> Double {
        func degreesToRadians(_ degrees: Double) -> Double { return degrees * Double.pi / 180.0 }
        func radiansToDegrees(_ radians: Double) -> Double { return radians * 180.0 / Double.pi }

        let lat1 = degreesToRadians(latitude)
        let lon1 = degreesToRadians(longitude)

        let lat2 = degreesToRadians(point.latitude);
        let lon2 = degreesToRadians(point.longitude);

        let dLon = lon2 - lon1;

        let y = sin(dLon) * cos(lat2);
        let x = cos(lat1) * sin(lat2) - sin(lat1) * cos(lat2) * cos(dLon);
        let radiansBearing = atan2(y, x);

        return radiansToDegrees(radiansBearing)
    }
}

Answer 4

you can use my code.. it's work on my project with microcontroller that use GPS for data.

#define d2r ((22/7.0)/180.0)
#define r2d (180.0/(22/7.0))

double get_heading1(double lat1, double long1, double lat2, double long2)  
{
    double diff_lat, diff_long;
    double degree;

    diff_long =(double) (((long2*1000000)-(long1*1000000))/1000000) * d2r;
    diff_lat = (double) (((lat2*1000000)-(lat1*1000000))/1000000) * d2r;     

    degree = r2d     (atan2(sin(diff_long)*cos(d2r*lat2),cos(d2r*lat1)*sin(d2r*lat2)-sin(d2r*lat1)*cos(d2r*lat2)    *cos(diff_long)));

    if (degree >= 0) {
        return degree;
    } else {
        return 360+degree;
    }                                                                 
}