How to Calculate the Bearing Between Two Coordinates in Python

On this page

Problem statement

The bearing from one point to another is the compass direction you would set off in: 0° north, 90° east, 180° south, 270° west. It is what a navigation display shows, what orients an arrow on a map, and what groups trips or wind or migration by direction. On a sphere it has a twist that surprises people: the direction you start in is not the direction you arrive in, because a great circle crosses each meridian at a different angle, so a flight from New York to London leaves heading north-east and arrives heading south-east.

This page computes the initial bearing with the standard formula, checks it against the exact azimuth on the WGS84 ellipsoid, converts bearings to compass points, shows how much the initial and final bearings differ for real routes between OurAirports airports, and computes the bearings from one airport to every other at once.

Quick answer

Use atan2 on the formula for the initial great-circle bearing and normalise to 0–360°:

import numpy as np

def bearing(lat1, lon1, lat2, lon2):
    p1, p2, dl = np.radians(lat1), np.radians(lat2), np.radians(lon2 - lon1)
    x = np.sin(dl) * np.cos(p2)
    y = np.cos(p1) * np.sin(p2) - np.sin(p1) * np.cos(p2) * np.cos(dl)
    return (np.degrees(np.arctan2(x, y)) + 360) % 360

jfk, lhr = AP.loc["JFK"], AP.loc["LHR"]
b = bearing(jfk.latitude_deg, jfk.longitude_deg, lhr.latitude_deg, lhr.longitude_deg)
az = geod.inv(jfk.longitude_deg, jfk.latitude_deg, lhr.longitude_deg, lhr.latitude_deg)[0] % 360
print(f"JFK to LHR: initial bearing {b:.2f}° (sphere), {az:.2f}° (WGS84 ellipsoid)")

JFK to LHR: initial bearing 51.35° (sphere), 51.38° (WGS84 ellipsoid)

A flight from New York JFK to London Heathrow sets off at a bearing of 51.35° on the sphere and 51.38° on the ellipsoid: north-east, a few hundredths of a degree apart.

World map of great circles from New York JFK to five airports with their initial bearings: London 51 degrees north-east, Paris 53 degrees north-east, Tokyo 333 degrees north-north-west over the Arctic, Sao Paulo 153 degrees south-south-east, and Sydney 266 degrees, due west. — Initial bearings from JFK.

Compass points, changing bearings, sphere against ellipsoid

1. From degrees to compass points

Divide the circle into 16 sectors of 22.5°, centred on the named directions:

POINTS = "N NNE NE ENE E ESE SE SSE S SSW SW WSW W WNW NW NNW".split()
def compass(deg):
    return np.array(POINTS)[(np.round(np.asarray(deg) / 22.5) % 16).astype(int)]
for dest in ("LHR", "CDG", "NRT", "GRU", "SYD"):
    d = AP.loc[dest]
    b = bearing(jfk.latitude_deg, jfk.longitude_deg, d.latitude_deg, d.longitude_deg)
    print(f"JFK to {dest}: {b:6.1f}° {compass(b)}")

JFK to LHR:   51.4° NE
JFK to CDG:   53.5° NE
JFK to NRT:  332.6° NNW
JFK to GRU:  153.2° SSE
JFK to SYD:  266.0° W

From JFK, London and Paris are north-east, Tokyo is north-north-west, over the Arctic, São Paulo is south-south-east, and Sydney is due west, 266.0°: great circles rarely point where a flat map suggests.

2. The bearing changes along the way

The final bearing is the reverse of the bearing from the destination back to the start. Compare initial and final bearings for a few long routes:

for a, z in (("JFK", "LHR"), ("LHR", "NRT"), ("LAX", "SYD"), ("JFK", "SIN")):
    A, Z = AP.loc[a], AP.loc[z]
    ini = bearing(A.latitude_deg, A.longitude_deg, Z.latitude_deg, Z.longitude_deg)
    fin = (bearing(Z.latitude_deg, Z.longitude_deg, A.latitude_deg, A.longitude_deg) + 180) % 360
    print(f"{a} to {z}: leaves at {ini:6.1f}° ({compass(ini)}), arrives at {fin:6.1f}° ({compass(fin)}), turn {((fin - ini + 540) % 360) - 180:+.1f}°")

JFK to LHR: leaves at   51.4° (NE), arrives at  107.9° (ESE), turn +56.6°
LHR to NRT: leaves at   30.9° (NNE), arrives at  156.8° (SSE), turn +125.9°
LAX to SYD: leaves at  241.0° (WSW), arrives at  241.0° (WSW), turn +0.0°
JFK to SIN: leaves at    3.3° (N), arrives at  177.5° (S), turn +174.2°

The heading changes along the way. JFK to London leaves heading north-east and arrives heading east-south-east, a turn of 56.6°; JFK to Singapore leaves almost due north and arrives almost due south, a turn of 174.2°, because the route passes close to the pole. Los Angeles to Sydney is the exception: its endpoints sit at nearly opposite latitudes, so the great circle is symmetric about where it crosses the equator and leaves and arrives at the same bearing, 241.0°.

The bearing changes along the way.

3. Spherical against ellipsoidal bearings

How large is the sphere's error in direction? Compare with pyproj.Geod azimuths for 20,000 random pairs:

rng = np.random.default_rng(3)
i, j = rng.integers(0, len(AP), (2, 20000))
a, z = AP.iloc[i], AP.iloc[j]
sph = bearing(a.latitude_deg.values, a.longitude_deg.values, z.latitude_deg.values, z.longitude_deg.values)
ell = geod.inv(a.longitude_deg.values, a.latitude_deg.values, z.longitude_deg.values, z.latitude_deg.values)[0] % 360
diff = (sph - ell + 540) % 360 - 180
ok = i != j
print(f"pairs: {ok.sum():,}; difference: median {np.median(np.abs(diff[ok])):.3f}°, 99th percentile {np.percentile(np.abs(diff[ok]), 99):.3f}°, max {np.abs(diff[ok]).max():.2f}°")

pairs: 19,992; difference: median 0.079°, 99th percentile 0.488°, max 6.68°

For 19,992 random pairs, the spherical bearing differs from the ellipsoidal azimuth by a median of 0.079°, and by less than 0.49° for 99% of pairs. The largest difference, 6.68°, is for U-Tapao in Thailand to Lima, 19,814 km apart and close to antipodal, where the direction of the shortest path is very sensitive to the shape of the Earth.

Directions of every airport from one hub

The data used on this page

Coordinates of the large airports with an IATA code in OurAirports, a public-domain database. Exact azimuths come from pyproj.Geod on WGS84.

from pathlib import Path

import numpy as np
import pandas as pd
import requests
from pyproj import Geod

Path("data").mkdir(exist_ok=True)
p = Path("data/airports.csv")
if not p.exists():
    r = requests.get("https://davidmegginson.github.io/ourairports-data/airports.csv", timeout=300)
    r.raise_for_status()
    p.write_bytes(r.content)
AP = pd.read_csv(p, keep_default_na=False, na_values=[""])
AP = AP[(AP.type == "large_airport") & AP.iata_code.notna()].set_index("iata_code")
geod = Geod(ellps="WGS84")
print(f"{len(AP):,} large airports")

1,172 large airports

The data has 1,172 large airports with an IATA code.

Example 1: the directions of every airport from one hub

Bearings work on arrays. Count large airports by compass sector as seen from Frankfurt:

fra = AP.loc["FRA"]
b = bearing(fra.latitude_deg, fra.longitude_deg, AP.drop("FRA").latitude_deg.values, AP.drop("FRA").longitude_deg.values)
sectors = pd.Series(compass(b)).value_counts().reindex(POINTS, fill_value=0)
print(sectors.to_string())
print(f"most common direction from FRA: {sectors.idxmax()} ({sectors.max()} airports)")

N       23
NNE     33
NE     107
ENE    154
E      158
ESE     90
SE      88
SSE     62
S       52
SSW     35
SW      74
WSW     40
W       62
WNW    119
NW      57
NNW     17
most common direction from FRA: E (158 airports)

Seen from Frankfurt, most large airports lie east to east-north-east, towards Asia, and the next largest group west-north-west, towards North America; the fewest lie due north and north-north-west.

Rose diagram of large airports by initial bearing from Frankfurt in 16 compass sectors: the longest petals point east and east-north-east, 158 and 154 airports, and west-north-west, 119, with few airports due north. — Directions from Frankfurt.

The initial bearing from spherical trigonometry

Sphere or ellipsoid?

The initial bearing is the angle between the meridian through the starting point and the great circle to the destination, measured clockwise from north. The formula comes from spherical trigonometry; atan2 returns the angle with the correct quadrant, in −180° to 180°, which the modulo turns into a compass bearing. Because meridians converge at the poles, a great circle crosses each one at a different angle, so except along a meridian or the equator, the bearing changes continuously along the route. A ship or plane that held its initial bearing would follow a rhumb line instead, a longer spiral towards the pole.

On the ellipsoid the azimuth differs slightly from the spherical bearing, by hundredths of a degree for most pairs and more near the poles and for near-antipodal points. For drawing arrows, binning directions or reporting a heading to the nearest degree, the spherical formula is more than enough; for surveying, use the ellipsoidal azimuth from pyproj.Geod.inv, which also returns the back azimuth.

Magnetic north, coincident points and poles

Before using a bearing.

Magnetic north: these are true bearings; a compass points to magnetic north, which differs by the local declination, more than 10° in parts of North America.
Coincident points: the bearing between identical points is undefined; atan2(0, 0) returns 0.
Poles: at a pole every direction is south or north; bearings from a pole are not meaningful.
Averaging bearings: average directions as unit vectors, not as numbers; the mean of 350° and 10° is 0°, not 180°.

Internal links

How to calculate the distance between two coordinates in Python — the companion formula
How to find a destination point from distance and bearing in Python — the inverse problem
Great circle and rhumb line explained — why the bearing changes

FAQ

How do I calculate the bearing between two points in Python?

Convert to radians and compute θ = atan2(sin Δλ · cos φ2, cos φ1 · sin φ2 − sin φ1 · cos φ2 · cos Δλ), then convert to degrees and take (θ + 360) mod 360.

Why is the bearing from A to B not the reverse of B to A?

Because the great circle crosses meridians at different angles at each end; the final bearing at B is the bearing from B to A plus 180°.

How do I convert a bearing to N, NE, E?

Divide by 22.5° for 16 points (or 45° for 8), round, and index a list of names.

Does pyproj give the bearing?

Yes: Geod.inv returns the forward azimuth, back azimuth and distance on the ellipsoid.