How to Find a Destination Point from Distance and Bearing in Python

On this page

Problem statement

Start at a point, travel a given distance in a given direction: where do you end up? The question behind range rings around an airport, search areas around a last known position, sensor footprints and "everything within 500 km" circles has a simple answer on a sphere and an exact one on the ellipsoid. Getting it wrong in degrees produces circles that are ellipses and rings that are hundreds of kilometres off at high latitudes.

This page computes destination points with the spherical formula and with pyproj.Geod.fwd, compares them, draws true range rings around real airports from OurAirports, shows what a "circle" drawn in degrees really covers, and finds which large airports are within an aircraft's range of a hub.

Quick answer

Use pyproj.Geod.fwd for the exact point on WGS84, or the spherical formula for a fast approximation:

import numpy as np

def destination(lat, lon, bearing_deg, dist_km, r=6371.0088):
    p1, l1, b, d = np.radians(lat), np.radians(lon), np.radians(bearing_deg), dist_km / r
    p2 = np.arcsin(np.sin(p1) * np.cos(d) + np.cos(p1) * np.sin(d) * np.cos(b))
    l2 = l1 + np.arctan2(np.sin(b) * np.sin(d) * np.cos(p1), np.cos(d) - np.sin(p1) * np.sin(p2))
    return np.degrees(p2), (np.degrees(l2) + 540) % 360 - 180

lhr = AP.loc["LHR"]
lat2, lon2 = destination(lhr.latitude_deg, lhr.longitude_deg, 90, 1000)
lon_e, lat_e, _ = geod.fwd(lhr.longitude_deg, lhr.latitude_deg, 90, 1000e3)
print(f"1,000 km due east of LHR: sphere {lat2:.4f}, {lon2:.4f}; ellipsoid {lat_e:.4f}, {lon_e:.4f}")
1,000 km due east of LHR: sphere 50.5946, 13.7955; ellipsoid 50.5977, 13.7515

Setting off due east from Heathrow and travelling 1,000 km along a great circle ends at about 50.6°N, 13.8°E, near the German-Czech border. Note the latitude: a great circle that starts heading east immediately curves south of the starting parallel, here by almost 0.9°. The spherical and ellipsoidal answers differ by 0.04° of longitude.

Map of western Europe in longitude and latitude with true range rings of 500 and 2,000 kilometres around Heathrow, both stretched east-west into egg shapes, and a dashed 5 degree buffer that is round on the map but reaches 557 kilometres north and only 347 east. Map of western Europe in longitude and latitude with true range rings of 500 and 2,000 kilometres around Heathrow, both stretched east-west into egg shapes, and a dashed 5 degree buffer that is round on the map but reaches 557 kilometres north and only 347 east.
Range rings and a degree buffer.

Spherical against exact answers, range rings, degree circles

1. How far apart are the two answers?

Compare the spherical destination with the exact one for many bearings and distances, measuring the gap with the exact distance between them:

rows = []
for dist in (10, 100, 1000, 5000):
    b = np.arange(0, 360, 15)
    la, lo = destination(lhr.latitude_deg, lhr.longitude_deg, b, dist)
    lo_e, la_e, _ = geod.fwd(np.full(len(b), lhr.longitude_deg), np.full(len(b), lhr.latitude_deg), b, np.full(len(b), dist * 1e3))
    gap = geod.inv(lo, la, lo_e, la_e)[2] / 1000
    rows.append((dist, gap.max(), gap.max() / dist * 100))
for dist, g, pct in rows:
    print(f"{dist:>5} km from LHR: largest gap {g:7.3f} km ({pct:.2f}% of the distance)")
   10 km from LHR: largest gap   0.032 km (0.32% of the distance)
  100 km from LHR: largest gap   0.318 km (0.32% of the distance)
 1000 km from LHR: largest gap   3.176 km (0.32% of the distance)
 5000 km from LHR: largest gap  16.787 km (0.34% of the distance)

The spherical formula drifts from the exact answer by about 0.32% of the distance in the worst direction, whatever the distance: 32 m after 10 km, 3.2 km after 1,000 km, 16.8 km after 5,000 km.

2. Range rings

A ring of destination points at every bearing is a true circle on the Earth. Build one as a polygon and check its area against a flat circle of the same radius:

from shapely.geometry import Polygon
def ring(lat, lon, km, n=360):
    b = np.linspace(0, 360, n, endpoint=False)
    lo, la, _ = geod.fwd(np.full(n, lon), np.full(n, lat), b, np.full(n, km * 1e3))
    return Polygon(zip(lo, la))
for km in (500, 2000):
    poly = ring(lhr.latitude_deg, lhr.longitude_deg, km)
    area = abs(geod.geometry_area_perimeter(poly)[0]) / 1e6
    print(f"{km:>4} km ring around LHR: {area:,.0f} km² on the ellipsoid; flat circle {np.pi * km ** 2:,.0f} km²; "
          f"spans {poly.bounds[0]:.1f}° to {poly.bounds[2]:.1f}° E, {poly.bounds[1]:.1f}° to {poly.bounds[3]:.1f}° N")
 500 km ring around LHR: 784,957 km² on the ellipsoid; flat circle 785,398 km²; spans -7.7° to 6.7° E, 47.0° to 56.0° N
2000 km ring around LHR: 12,463,294 km² on the ellipsoid; flat circle 12,566,371 km²; spans -30.1° to 29.2° E, 33.5° to 69.4° N

The 500 km ring around Heathrow encloses 784,957 km², almost exactly a flat circle's area. The 2,000 km ring encloses 12.46 million km², 0.8% less than a flat circle, because a circle drawn on a curved surface holds less area than one on a plane. In degrees the 2,000 km ring spans 59° of longitude but 36° of latitude: on a latitude-longitude map it looks stretched, though every point on it is exactly 2,000 km away.

Largest gap between spherical and ellipsoidal destination points from Heathrow, by distance travelled. sphere against ellipsoid: largest gap, km10 km0.032 km100 km0.318 km1,000 km3.176 km5,000 km16.787 km06 km12 km18 kmAbout 0.32% of the distance in the worst direction spatialworkflow.io
How far the sphere drifts.

3. A "circle" in degrees is not a circle

Buffering a point by a number of degrees draws a circle in latitude-longitude space. Measure how far its edge really is in each direction:

deg = 5
for name, b, (dla, dlo) in (("north", 0, (deg, 0)), ("east", 90, (0, deg)), ("south", 180, (-deg, 0))):
    d = geod.inv(lhr.longitude_deg, lhr.latitude_deg, lhr.longitude_deg + dlo, lhr.latitude_deg + dla)[2] / 1000
    print(f"5° buffer edge, {name:>5} of LHR: {d:,.0f} km")
5° buffer edge, north of LHR: 557 km
5° buffer edge,  east of LHR: 347 km
5° buffer edge, south of LHR: 556 km

A 5° buffer around Heathrow reaches 557 km to the north and 556 km to the south but only 347 km to the east: at 51.5°N a degree of longitude is 62% of a degree of latitude, so the "circle" is really an ellipse on the ground.

Which airports are within range of a hub

The data used on this page

Coordinates of the large airports with an IATA code in OurAirports, a public-domain database. Exact destination points, distances and areas 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: which airports are within range of a hub?

An aircraft with a range of 9,000 km based in Dubai: which large airports can it reach non-stop, by exact distance?

dxb = AP.loc["DXB"]
d = geod.inv(np.full(len(AP), dxb.longitude_deg), np.full(len(AP), dxb.latitude_deg), AP.longitude_deg.values, AP.latitude_deg.values)[2] / 1000
inside = pd.Series(d < 9000, index=AP.index)
print(f"large airports within 9,000 km of DXB: {inside.sum() - 1} of {len(AP) - 1}")
print("by continent:", AP[inside].drop("DXB").continent.value_counts().to_dict())
print(f"farthest within range: {pd.Series(d, index=AP.index)[inside].idxmax()} at {pd.Series(d, index=AP.index)[inside].max():,.0f} km")
large airports within 9,000 km of DXB: 856 of 1171
by continent: {'AS': 421, 'EU': 260, 'AF': 170, 'OC': 4, 'NA': 1}
farthest within range: YAP at 8,842 km

From Dubai, 856 of the other 1,171 large airports are within 9,000 km: nearly all of Asia, Europe and Africa, four in Oceania and one in North America. The farthest within range is Yap in Micronesia, 8,842 km away.

World map of a 9,000 kilometre range ring around Dubai drawn from destination points at every bearing, with 856 large airports inside in red, covering Europe, Africa and Asia, and 315 outside in grey, mostly in the Americas. World map of a 9,000 kilometre range ring around Dubai drawn from destination points at every bearing, with 856 large airports inside in red, covering Europe, Africa and Asia, and 315 outside in grey, mostly in the Americas.
Within range of Dubai.

The haversine problem in reverse

Computing a destination point from a start, a bearing and a distance. from start, bearing and distance to a pointstartLHR 51.47N 0.46Wbearing,distance90°, 1,000 kmangled / Rlatitude,longitudeasin, atan2destination50.60N 13.75EExact: Geod.fwd(lon, lat, az, metres); the sphere drifts 0.32% spatialworkflow.io
From start, bearing and distance to a point.

The destination formula is the haversine problem run in reverse: given a start, an initial bearing and an angular distance (the distance divided by the radius), spherical trigonometry gives the end point's latitude directly and its longitude through atan2. On the ellipsoid there is no closed form; Geod.fwd solves the direct geodesic problem with Karney's algorithm to nanometre accuracy and also returns the azimuth at the destination.

The practical lesson is about circles. A ring of points at the same distance in every direction is the true circle, and on a latitude-longitude map it looks like an egg, stretched east-west at high latitudes. A buffer of so many degrees, by contrast, is a true circle only in degree space; on the ground its east-west radius is shorter than its north-south radius by the cosine of the latitude. Build rings with destination points, or buffer in a projected CRS centred on the point.

Rings across the antimeridian and sphere error

Destination points and range rings: problems and fixes. range ring trapsTrapInsteadbuffer of N degreesring of destination pointskilometres into Geod.fwdit takes metresring crosses 180°split before making a polygonflat circle areaGeod area: 0.8% less at 2,000 kmA 5° buffer at LHR reaches 557 km north, 347 km east spatialworkflow.io
Range ring traps.
  • Rings that cross the antimeridian or enclose a pole need splitting before they become valid polygons in longitude-latitude.
  • Sphere error: the spherical destination drifts by a few tenths of a percent of the distance; fine for display, not for surveying.
  • Many points at once: both destination above and Geod.fwd accept arrays.
  • Units: Geod.fwd takes metres, not kilometres.

FAQ

How do I find a point at a given distance and bearing in Python?

Use pyproj.Geod(ellps="WGS84").fwd(lon, lat, azimuth, distance_m), which returns the destination longitude, latitude and back azimuth.

How do I draw a circle of a given radius around a coordinate?

Compute destination points at bearings from 0° to 360° and join them into a polygon; on a latitude-longitude map it will look stretched east-west.

Why does my circle look like an ellipse on the map?

Because a true circle on the Earth spans more degrees of longitude than of latitude away from the equator; it is the map, not the circle, that is distorted.

How accurate is the spherical formula?

Within a few tenths of a percent of the distance travelled, as step 1 measures.