Calculation of direction between two geographical points

1.3k Views Asked by At

I am trying to capture the relative position (bearing) from one geographical point to another.

I found an online tool which does exactly what I am looking for but for some reason I am not able to translate the same logic in simple python function.

The formula being used is as follows (apologies, I am new here and not sure how to covert this into proper math format so I am inserting the screenshot): formula

I have tried to translate this into python code as per my understanding but my output is nowhere close to the actual output.

Python code:

import math

def calc_bearing(pointA, pointB):
    latA = pointA[0]
    latB = pointB[0]
    lonA = pointA[1]
    lonB = pointB[1]

    delta_ratio = math.log(math.tan(latB / 2 + math.pi / 4) / math.tan(latA / 2 + math.pi / 4))
    delta_lon = abs(lonA - lonB)

    if delta_lon > 180:
        delta_lon = delta_lon % 180
    bearing = math.atan2(delta_lon, delta_ratio)
    return bearing

Two geo points for testing, in decimal coordinates:

g_point1 = (51.51154, -0.0029163)

g_point2 = (35.6709056, 139.7577372)

Expected output is:

98.97 ° degrees

The natural log is also failing since the value is in negative. I am new to trigonometry and geo data plotting.

Any help would be appreciated.

1

There are 1 best solutions below

10
On BEST ANSWER

It seems you're using degrees to measure angles, but the trig functions are defined for radian measure, so you just need a small adjustment as I show below. I think you can also omit the if statement and just use delta_lon %= math.pi, thanks @enzotip for correcting my mistake regarding the conversion.

import math

def calc_bearing(pointA, pointB):
    deg2rad = math.pi / 180
    latA = pointA[0] * deg2rad 
    latB = pointB[0] * deg2rad 
    lonA = pointA[1] * deg2rad 
    lonB = pointB[1] * deg2rad 

    delta_ratio = math.log(math.tan(latB/ 2 + math.pi / 4) / math.tan(latA/ 2 + math.pi / 4))
    delta_lon = abs(lonA - lonB)

    delta_lon %= math.pi
    bearing = math.atan2(delta_lon, delta_ratio)/deg2rad
    return bearing


g_point1 = (51.51154, -0.0029163)
g_point2 = (35.6709056, 139.7577372)

print(calc_bearing(g_point1, g_point2))

Try it online!