Geopy: calculating GPS heading / bearing

Question

First time poster here.

I am doing some data analyses on collected GPS data for a bridge inspection ROV octorotor. We have the octorotor running on ROS using a 3D scanning LIDAR, stereo vision, INS, and some other neat tech. I'm currently using a ublox LEA-6T in a similar setup as Doug Weibel's setup to collect raw GPS data like carrier phase, doppler shift, and satellite ephemeris. Then I use an opensource project RTKLIB to do some DGPS post processing with local NOAA CORS stations to obtain cm accuracy for better pose estimation when reconstructing the 3D point cloud of the bridge.

Anyhow, I'm using most of scipy to statistically verify my test results.
Specifically for this portion though, I'm just using:

• python-3.3
• numpy
• geopy

I've been studding my positional covariance with respect to offset from my measured ground truth using geopy's handy distance function. With little massaging the arguments, I can find the distance respect to each direction depicted by each standard deviation element in the matrix; North, East, Up and the three directions between.

However, these distances are absolute and do not describe direction.
Say: positive, negative would correlate to northward or southward respectively.

I could simply use the latitude and longitude to detect polarity of direction,
But I'd like to be able to find the precise point to point bearing of the distance described instead,
As I believe a value of global heading could be useful for further applications other than my current one.

I've found someone else pose a similar question
But its seem to be assuming a great circle approximation
Where I would prefer using at least the WGS-84 ellipsoidal model, or any of the same models that can be used in geopy:
Jump to Calculating distances

Any suggestion appreciated,
-ruffsl

Sources if interested:

• python-3.3: http:// www.python.org/download/releases/3.3.0/
• numpy: http:// www.numpy.org/
• geopy: https:// code.google.com/p/geopy/
• scipy: http:// www.scipy.org/
• ublox LEA-6T: http:// www.u-blox.com/en/gps-modules/u-blox-6-timing-module/lea-6t.html
• Doug Weibel's: http:// diydrones.com/profiles/blogs/proof-of-concept-test-extremely-accurate-3d-velocity-measurement
• RTKLIB: http:// www.rtklib.com/
• NOAA CORS: http:// geodesy.noaa.gov/CORS/
• ROS: http:// www.ros.org/wiki/

Solution 1

Use the geographiclib package for python. This computes distances and bearings on the ellipsoid and much more. (You can interpolate paths, measure areas, etc.) For example, after

pip install geographiclib

you can do

>>> from geographiclib.geodesic import Geodesic
>>> Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
{'lat1': -41.32, 'a12': 179.6197069334283, 's12': 19959679.26735382, 'lat2': 40.96, 'azi2': 18.825195123248392, 'azi1': 161.06766998615882, 'lon1': 174.81, 'lon2': -5.5}

This computes the geodesic from Wellington, New Zealand (41.32S 174.81E) to Salamanca, Spain (40.96N 5.50W). The distance is given by s12 (19959679 meters) and the initial azimuth (bearing) is given by azi1 (161.067... degrees clockwise from north).

Solution 2

Bearing between two lat/long coordinates: (lat1, lon1), (lat2, lon2)

In the code below lat1,lon1,lat2,lon2 are asumend to be in radians.
Convert before from degrees to radians.

dLon = lon2 - lon1;
y = Math.sin(dLon) * Math.cos(lat2);
x = Math.cos(lat1)*Math.sin(lat2) -
        Math.sin(lat1)*Math.cos(lat2)*Math.cos(dLon);
brng = Math.atan2(y, x).toDeg();

Bearing is now in range -180/180.

to normalize to compass bearing (0-360)

if brng < 0: brng+= 360

Solution 3

@AlexWien's answer in Python

import math, numpy as np

def get_bearing(lat1,lon1,lat2,lon2):
    dLon = lon2 - lon1;
    y = math.sin(dLon) * math.cos(lat2);
    x = math.cos(lat1)*math.sin(lat2) - math.sin(lat1)*math.cos(lat2)*math.cos(dLon);
    brng = np.rad2deg(math.atan2(y, x));
    if brng < 0: brng+= 360
    return brng