Geopy: GPS / bearing heading calculation

Question

The first poster is here.

I am doing some analysis of GPS data collected to check the ROV octorotor. We have an octorotor working on ROS using 3D LIDAR scanning, stereo vision, INS and some other neat technology. I am currently using ublox LEA-6T in the same setup as Doug Weibel to collect raw GPS data such as carrier phase, Doppler shift and satellite ephemeris. Then I use the open source project RTKLIB to do some post-processing DGPS with local NOAA CORS , to see the accuracy see for better posture estimation when reconstructing a 3D point cloud of the bridge.

Anyway, I use most of scipy to statistically check my test results.
In particular, for this part, I just use:

• python-3.3
• numpy
• geopy

I studied my positional covariance regarding offset from my measured truth using a convenient distance function. By massaging the arguments a little, I can find the distance corresponding to each direction indicated by each element of the standard deviation in the matrix; North, East, Up and three directions between them.

However, these distances are absolute and do not describe the direction. Say: positive, negative will correlate with north or south, respectively.

I could just use latitude and longitude to detect the polarity of direction,
But I would like to find the exact point bearing of the described distance,
Since I believe that the value of the global header may be useful for other applications than my current one.

I found someone else asked a similar question
But it seems that this involves a large circle approaching
Where I would prefer to use at least the ellipsoidal model WGS-84 or any of the same models that can be used in geophysics:
Go to distance calculation

Any suggestion appreciated,
-ruffsl

Sources if interested:

• python-3.3: http://www.python.org/download/releases/3.0.0/
• numpy: http://www.numpy.org/
• geography: 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/

Answer 1

Use geographiclib package for python. It calculates distances and bearings on an 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} 

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

Answer 2

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

In the code below, lat1, lon1, lat2, lon2 are expected to be in radians.
Convert to degrees in 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(); 

The bearing is now in the range -180/180.

for normalizing to a compass bearing (0-360)

  if (brng < 0) { brng+= 360; }