Python code to calculate angle between three points (lat long coordinates)

Question

Can anybody suggest how to calculate angle between three points (lat long coordinates)

A : (12.92473, 77.6183)
B : (12.92512, 77.61923)
C : (12.92541, 77.61985)

Gerenuk

Possible duplicate of How to calculate an angle from three points? – Thereinafter

This will work for lat long too ? – Gerenuk

It might not be an exact solution, but your points seem close enough. In this case we can approximate the spherical surface of the earth to a plane. – Thereinafter

1. points are "close enough": use a map projection and compute the angle, using e.g. Cosines Law, using the Cartesian coordinates from the map. — 2: points are not close enough, you have to use great circles and spherical trigonometry – Drury

You cannot use those values as planar points and get the angle. For one, the earth is not flat, though it is approximately flat at a small scale. More important, the units for longitude and latitude are not the same. At the equator, a degree of latitude or longitude are the same. But move away from the equator, and a degree of longitude covers less distance than the degree of latitude. At (very close to) the north or south pole, a degree of longitude covers zero distance, while a degree of latitude is still about 69 miles. So this question is not a duplicate of that other question. – Charity

I second Rory's remark, that explicitly says what I said obliquely in my remark above (use map projection) – Drury

Actually, The problem statement is that if We already have a list of lat long. Now if I want to provide a navigation to a delivery boy. how to do the same ? – Gerenuk

@Gerenuk Does your "delivery boy" use a helicopter or is he bound to a road network? – Drury

@Drury : Road Network only – Gerenuk

@Drury it should be possible to get the angle between the normal lines at the Earth's surface, without needing to do a map projection. A projection might throw the angle off because projections are never perfect. – Alialia

OK, road network only. Now please ask yourself "What the … this angle is needed for?" ,--- Should you give yourself a positive answer, you can map the spherical coordinates if the points are close and use CL, or use ST. – Drury

@Mark, we are close to be put in chat... however what you proposed IS a mapping, a formula that puts points from a sphere onto a leaf of paper. However for 3 points in a, say, 10 km neighborhood the mapping doesn't make a difference especially when the problem is "to provide a navigation to a delivery boy". – Drury

What is the practical problem you are trying to solve here? You mention "delivery boy" and "road network", in which case an angle would be relatively meaningless since road networks do not cover all angles; in the extreme you might need U-turns to get to a destination. – Perea

For example we have 3 coordinates A, B, C. A is starting point and C is destination point. Now If a delivery boy starts to walk from point A to B straight lets say (100 meter) then turn left and walk straight for 200 meter and finally reach to destination C. Now the problem is here how to know after reaching point B we need to take turn left ? – Gerenuk

@Drury as the question is being refined it seems that an approximation is probably good enough, making your advice perfectly sound. But I think it's still an interesting theoretical question. Can we agree that there's a unique shortest arc that connects two points, and that the angle between two of those arcs is the angle between the normals at the joining point? – Alialia

But is this a real-life situation? If you need real navigation, you can have your own instance of OSRM set up in 15 minutes and it will give you turn-by-turn directions (walking or driving) in 1 or 2 milliseconds for huge journeys, without any limit and for free. Graphhopper is similar. If you're talking real-road, there is no hope in your current approach (my full-time job is in vehicle routing problems). Then you just need requests to communicate with your server. – Perea

Check this out Calculate distance, bearing and more between Latitude/Longitude points – Fungosity

Answer 1

I see two main ways to solve your problem, assuming you want angle ABC (B is the vertex of the angle). Since your three points are close to each other (less than 0.0007° latitude and 0.002° longitude apart), we can approximate the earth as a plane and use two-dimensional vector calculations. A degree of longitude and of latitude are not the same distance when we are away from the equator, but we can adjust for that. Another solution is to treat your points as in three-dimensional space and use three-dimensional vector calculations. Here we just need to convert the given spherical coordinates to 3D Cartesian coordinates.

Here is my code for your problem. I use the numpy module here for convenience, but this can be done pretty easily without it. This code is fairly wordy so you can see better just what is being done.

import numpy as np
import math

def latlong_to_3d(latr, lonr):
    """Convert a point given latitude and longitude in radians to
    3-dimensional space, assuming a sphere radius of one."""
    return np.array((
        math.cos(latr) * math.cos(lonr),
        math.cos(latr) * math.sin(lonr),
        math.sin(latr)
    ))

def angle_between_vectors_degrees(u, v):
    """Return the angle between two vectors in any dimension space,
    in degrees."""
    return np.degrees(
        math.acos(np.dot(u, v) / (np.linalg.norm(u) * np.linalg.norm(v))))

# The points in tuple latitude/longitude degrees space
A = (12.92473, 77.6183)
B = (12.92512, 77.61923)
C = (12.92541, 77.61985)

# Convert the points to numpy latitude/longitude radians space
a = np.radians(np.array(A))
b = np.radians(np.array(B))
c = np.radians(np.array(C))

# Vectors in latitude/longitude space
avec = a - b
cvec = c - b

# Adjust vectors for changed longitude scale at given latitude into 2D space
lat = b[0]
avec[1] *= math.cos(lat)
cvec[1] *= math.cos(lat)

# Find the angle between the vectors in 2D space
angle2deg = angle_between_vectors_degrees(avec, cvec)


# The points in 3D space
a3 = latlong_to_3d(*a)
b3 = latlong_to_3d(*b)
c3 = latlong_to_3d(*c)

# Vectors in 3D space
a3vec = a3 - b3
c3vec = c3 - b3

# Find the angle between the vectors in 2D space
angle3deg = angle_between_vectors_degrees(a3vec, c3vec)


# Print the results
print('\nThe angle ABC in 2D space in degrees:', angle2deg)
print('\nThe angle ABC in 3D space in degrees:', angle3deg)

This gives the results

The angle ABC in 2D space in degrees: 177.64369006

The angle ABC in 3D space in degrees: 177.643487338

Note that the results are very close (off by about one five-thousandth of a degree), as expected for three points so close together.

Charity

Answer 2

To get angle between two directions in lat/lon system, you can use difference of two bearings from this page:

Formula:    
θ = atan2( sin Δλ ⋅ cos φ2 , cos φ1 ⋅ sin φ2 − sin φ1 ⋅ cos φ2 ⋅ cos Δλ )
where   
φ1,λ1 is the start point, φ2,λ2 the end point (Δλ is the difference in longitude)

JavaScript:
(all angles in radians) 
var y = Math.sin(λ2-λ1) * Math.cos(φ2);
var x = Math.cos(φ1)*Math.sin(φ2) -
        Math.sin(φ1)*Math.cos(φ2)*Math.cos(λ2-λ1);
var brng = Math.atan2(y, x).toDegrees();

Obsequious