Calculating points along a line between two points using geopy

Question

I am using geopy to map the location of a simulated object as it travels between two geographical points. I would like to be able to calculate the lat/lon of this object at any point in its theoretical journey: so for an object travelling between Paris and New York, 0% would be Paris, 100% would be New York and 50% would be half way between the two when the object has travelled half the distance.

(I say a theoretical object because I am not interested in adjusting for eg the route that an airplane might take or tracking a real object, I just want to assume a straight line between the two points)

My code is as follows:

from geopy.distance import geodesic as GD

start = (48.8567, 2.3508) # Paris
end = (40.7128, 74.0060) # New York

distance_between = GD(start, end).km

print (f"{distance_between} km between Paris and New York")

pct_travelled = 50

new_lat = start[0] + (end[0] - start[0]) * (int(pct_travelled) / 100)
new_lon = start[1] + (end[1] - start[1]) * (int(pct_travelled) / 100)

print (f"Your current location is {new_lat}, {new_lon}")

distance_travelled_so_far = GD(start, (new_lat, new_lon)).km
distance_still_to_travel = GD((new_lat, new_lon), end).km

print (f"You have travelled {distance_travelled_so_far} km")
print (f"You have {distance_still_to_travel} km left to travel")

This code delivers an incorrect answer:

5529.689905151459 km between Paris and New York
Your current location is 44.78475, 38.178399999999996
You have travelled 2744.973813842307 km
You have 2943.5970959509054 km left to travel

(It is incorrect because in a correct answer the distance travelled and the distance left to travel at 50% completed would be equal)

I presume this is because my code assumes that the distance of a degree of lat and a degree of lon are the same, which is obviously not correct.

How can I do this correctly with geopy or other python libraries?

Answer 1

For the benefit of future generations (and myself when I Google this in future) here is the code I got working:

import geopy.distance
import math


def calculate_initial_compass_bearing(start, end):
    if (type(start) != tuple) or (type(end) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = math.radians(start[0])
    lat2 = math.radians(end[0])

    diffLong = math.radians(end[1] - start[1])

    x = math.sin(diffLong) * math.cos(lat2)
    y = math.cos(lat1) * math.sin(lat2) - (math.sin(lat1)
                                           * math.cos(lat2) * math.cos(diffLong))

    initial_bearing = math.atan2(x, y)

    # Now we have the initial bearing but math.atan2 return values
    # from -180° to + 180° which is not what we want for a compass bearing
    # The solution is to normalize the initial bearing as shown below
    initial_bearing = math.degrees(initial_bearing)
    compass_bearing = (initial_bearing + 360) % 360

    return compass_bearing


def distance_between(start, end):
    return geopy.distance.distance(start, end).km


def location_at_pct(
    start,
    end,
    pct_travelled
):
    bearing = calculate_initial_compass_bearing(start, end)

    distance_travelled = distance_between(start, end) * (pct_travelled / 100)

    current_location = geopy.distance.distance(
        kilometers=distance_travelled
    ).destination(
        start,
        bearing=bearing
    )

    return current_location

It's weird that geopy:

• does support calculating the distance between a start point and and end point
• does support calculating the location of an end point when given a start point, a distance and a bearing
• does not support calculating the bearing when given the start point and an end point