calculate a third point lat/lon from two points, distance and angle

83 Views Asked by At

I'm working on a geospatial problem where I have two points defined by their latitude and longitude on a map, and I need to calculate the position of a third point based on these. Here's the scenario:

Point A and Point B are given, each defined by latitude and longitude. I have a distance (in meters or kilometers) and an angle (in degrees). What I want to achieve is to calculate the position of Point C. Point C should be located at the given distance from Point B and at the specified angle relative to the line from Point A to Point B.

I tried to use this:

But when I enter to extend_perpendicular_line a 105 degree, I get 75 degrees and 90 is 64 etc..

from geopy import Point
from pyproj import Proj, transform

from geopy.distance import Geodesic, geodesic
import math




def extend_perpendicular_line(pointA, pointB, distance,direction = 1,angle=90):

    # Calculate the bearing from pointA to pointB
    bearing = calculate_bearing(pointA, pointB)

    # Calculate the bearing that is perpendicular to the initial bearing
    perp_bearing = (bearing + angle*direction) % 360

    # Create a geodesic line from pointB and extend it by the given distance at the perpendicular bearing
    extended_point = geodesic().destination(Point(pointB.latitude, pointB.longitude), perp_bearing, distance)

    anga = calculate_angle(pointA,pointB,extended_point)


    return extended_point


def calculate_bearing(pointA, pointB):

    lat1, lon1 = pointA.latitude, pointA.longitude
    lat2, lon2 = pointB.latitude, pointB.longitude

    diff_lon = math.radians(lon2 - lon1)
    lat1 = math.radians(lat1)
    lat2 = math.radians(lat2)

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

    initial_bearing = math.atan2(x, y)

    # Normalize the bearing
    initial_bearing = math.degrees(initial_bearing)
    bearing = (initial_bearing + 360) % 360

    return bearing


def calculate_angle(pointA, pointB, pointC):
    # Calculate the sides of the triangle
    a = calculate_distance(pointB, pointC)
    b = calculate_distance(pointA, pointC)
    c = calculate_distance(pointA, pointB)

    # Law of cosines: c^2 = a^2 + b^2 - 2ab * cos(C)
    # Solve for cos(C)
    cos_B = (a ** 2 + c ** 2 - b ** 2) / (2 * a * c)
    angle = math.degrees(math.acos(cos_B))



    return angle

def calculate_distance(pointA, pointB):
    # Haversine formula to calculate distance between two lat/lon     points
    coords_1 = (pointA.latitude, pointA.longitude)
    coords_2 = (pointB.latitude, pointB.longitude)


    return geodesic(coords_1, coords_2).km

pointA = Point(46.9540700,7.4474400)
pointB = Point(46.9560700,7.4494400)

extend_perpendicular_line(pointA,pointB,0.25,1,105)
1

There are 1 best solutions below

0
BitsAreNumbersToo On

As you question is stated, it's slightly ambiguous from what point the angle between line AB and C is calculated, either as AB to AC or AB to BC. Here are a couple diagrams that show this ambiguity.

AB to AC

AB to BC

Since the first is rarely going to have a unique solution, I have provided the solution to the second example here, using mostly functions built into pyproj.

import pyproj

def extend_angled_line(point_a, point_b, distance_from_b, angle, geo=None):
    geo = geo or pyproj.Geod(ellps='WGS84')
    bearing_a_to_b, _, d = geo.inv(point_a[1], point_a[0], point_b[1], point_b[0])  # deg, deg, m
    new_bearing = bearing_a_to_b + angle
    c_lon, c_lat, _ = geo.fwd(point_b[1], point_b[0], new_bearing, distance_from_b)
    return c_lat, c_lon

def _main():
    a = (45, 45)
    b = (45.1, 45.1)
    c = extend_angled_line(a, b, 1000, 180)
    print(c)

if __name__ == '__main__':
    _main()