Python for finding the GPS coordinates of a point x-distance between 2 other GPS coordinates

Question

I have 2 pairs of GPS coordinates A = (44.27364400,-121.17116400) and B = (44.27357900, -121.17006800) and a distance of 10m. I'm trying to find the GPS coordinates of the point that's 10m from A toward B. I know I can use haversine to find the distance between A and B coutesy of: https://stackoverflow.com/a/4913653/5369777 is there someway I can add or subtract the distance 10m to return the degrees for the new point?

Answer 1

Charles Karney wrote the powerfull GeographicLib program to solve all of this problems.

You can install the geographiclib Python package, it comes as a dependency with geopandas too.

from geographiclib.geodesic import Geodesic

A = (44.27364400, -121.17116400) #Point A (lat, lon)
B = (44.27357900, -121.17006800) #Point B (lat, lon)
s = 10 #Distance (m)

#Define the ellipsoid
geod = Geodesic.WGS84

#Solve the Inverse problem
inv = geod.Inverse(A[0],A[1],B[0],B[1])
azi1 = inv['azi1']
print('Initial Azimuth from A to B = ' + str(azi1))

#Solve the Direct problem
dir = geod.Direct(A[0],A[1],azi1,s)
C = (dir['lat2'],dir['lon2'])
print('C = ' + str(C))

Returns:

Initial Azimuth from A to B = 94.71831670212772
C = (44.27363659722271, -121.17103916871612)

Answer 2

I found that I needed compass bearing to get the correct direction for my distance. For bearing I found compassbearing.py @ https://gist.github.com/jeromer/2005586

The inputs of which are pointA, pointB). Output is bearing. From there I can use geopy Point and geopy.distance vincenty to figure out the distance and subsequent new coordinate.

import math from geopy import Point from geopy.distance import vincenty

def calculate_initial_compass_bearing(pointA, pointB):
    """
    Calculates the bearing between two points.
    The formulae used is the following:
    θ = atan2(sin(Δlong).cos(lat2),
              cos(lat1).sin(lat2) − 
    sin(lat1).cos(lat2).cos(Δlong))
    :Parameters:
    - `pointA: The tuple representing the 
    latitude/longitude for the
    first point. Latitude and longitude must be in 
    decimal degrees
    - `pointB: The tuple representing the latitude/longitude for the
    second point. Latitude and longitude must be in decimal degrees
    :Returns:
    The bearing in degrees
    :Returns Type:
    float
    """
    if (type(pointA) != tuple) or (type(pointB) != tuple):
        raise TypeError("Only tuples are supported as arguments")

    lat1 = math.radians(pointA[0])
    lat2 = math.radians(pointB[0])

    diffLong = math.radians(pointB[1] - pointA[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


bearing = calculate_initial_compass_bearing((44.27364400, -121.17116400),(44.27357900, -121.17006800))

distKm = .01
lat1 = 44.27364400
lon1 = -121.17116400

new_coord = vincenty(kilometers=distKm).destination(Point(lat1, lon1), bearing).format_decimal())