Best strategy for converting change in longitude and latitude to north, south, east, west in degrees

Question

This is not a class assignment but a conceptual question I'm struggling with for a personal project.

I need to convert a large paired group of long and lat coordinate values that show a two-dimensional directional change to degrees change (0-360).

What are the fastest ways to make this kind of transformation? I am hoping to avoid nesting a lot of conditional statements.

Answer 1

For accuracy better than the general tools use to gather the data is a Python module pyproj Feed your data into something like the odometer function extracted from waypoint.py

from pyproj import Geod  # sudo pip3 install pyproj
geoid = Geod(ellps='WGS84')  # See pyproj documentation for other fun options

def odometer(start, end):
    """Calculate bearings and distance between start and end.

    Arguments:
      start (tuple of float): latitude and longitude of starting position
      end (tuple of float): latitude and longitude of ending position

    Returns:
      float: True North bearing start -> end
      float: True North bearing end -> start
      float: distance start <-> end in meters
    """
    try:
        bearing_to, bearing_fro, distance = geoid.inv(*start + end)
    except Exception as error:  # more specific?
        print("odometer({}, {}) created error:{}".format(start, end, error))
        return None

        bearing_to %= 360
        bearing_fro %= 360
    return bearing_to, bearing_fro, distance

Answer 2

I used geographiclib to find the azymouth (0 - 360) between 2 points in the WGS84 ellipsoid using the Inverse method and filtering its results as follows:

from geographiclib.geodesic import Geodesic

def orientation(startpoint, endpoint): # each point is a tuple of (lon, lat)
    return (
        Geodesic.WGS84.Inverse(
            startpoint[1], 
            startpoint[0],
            endpoint[1],
            endpoint[0],
        )['azi1'] % 360
    )

The degrees can be converted to directions based on the compass points. Here is a more clean approach.

I have used the following code and verified some results with this online service:

from geographiclib.geodesic import Geodesic

def orientation(startpoint, endpoint): # each point is a tuple of (lon, lat)
    azi = (
        Geodesic.WGS84.Inverse(
            startpoint[1], 
            startpoint[0],
            endpoint[1],
            endpoint[0],
        )['azi1'] % 360
    )
    ranges = {
        (0, 11.25): 'N',
        (11.25, 33.75): 'NNE',
        (33.75, 56.25): 'NE',
        (56.25, 78.75): 'ENE',
        (78.75, 101.25): 'E',
        (101.25, 123.75): 'ESE',
        (123.75, 146.25): 'SE',
        (146.25, 168.75): 'SSE',
        (168.75, 191.25): 'S',
        (191.25, 213.75): 'SSW',
        (213.75, 236.25): 'SW',
        (236.25, 258.75): 'WSW',
        (258.75, 281.25): 'W',
        (281.25, 303.75): 'WNW',
        (303.75, 326.25): 'NW',
        (326.25, 348.75): 'NNW',
        (348.75, 360): 'N',
    }
    for i in ranges.keys():
        if i[0] < azi <= i[1]:
            return ranges[i]