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:

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/

You are talking much, but asking little. What exactly do you need? A function that calculates bearing between two lat/long WG84 coordinates? — AlexWien, Jul 15 '13 at 13:50
@AlexWien, Yes, that in short is what I am aiming for. I find it more helpful when looking for answers from past posts when people enplane themselves, its make for more relevant and searchable keywords for others. — ruffsl, Jul 15 '13 at 18:26

cffk · Accepted Answer · 2018-01-26T20:16:25.807

13

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).

edited Jan 26 '18 at 20:16

answered Aug 20 '15 at 20:11

cffk

1,839
1
12
17

2

It would be helpful if you give a bit more details on the result set of `Geodesic.WGS84.Inverse` and to actually calculate the GPS heading/bearing asked by the question. – kovac Jan 25 '18 at 03:52
This link to the docs explains the dictionary returned. https://geographiclib.sourceforge.io/1.52/python/interface.html – Levi Baguley Jul 06 '21 at 18:58

score 4 · Answer 2 · edited Oct 15 '18 at 12:51

4

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

edited Oct 15 '18 at 12:51

Mr_and_Mrs_D

32,208
39
178
361

answered Jul 15 '13 at 19:34

AlexWien

28,470
6
53
83

3

I do like the elegant math, but I believe this is using spherical coordinates, thus still the great circle approximation. For navigation and general direction application, including my own purpose for checking polarity, I would guess this is sufficient compared to the WG84. But for surveying and site inspection this might too rough, or not applicable for local episode estimations. Awesome though. – ruffsl Jul 15 '13 at 23:22
Yes these formulas asume the earth as a spheroid. More accurate you have to use ellipsoid based formulas / implementations. – AlexWien Jul 16 '13 at 12:09
9

i'm confused... why doesn't geopy have a function for computing this in line with the vincenty function which computes distance? – Michael Oct 07 '14 at 16:58

score 4 · Answer 3 · answered Nov 22 '18 at 15:27

@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

Geopy: calculating GPS heading / bearing

3 Answers3

Linked