# What tools in Python are available for doing great circle distance + line creation?

I need to use Python to create a great circle distance -- both a number, and preferably some kind of 'curve' that I can use to draw in a client-side map. I don't care about the format of the curve -- be it WKT, or a set of pairs of coordinates -- but just want to get the data out.

What tools are there out there? What should I use?

pyproj has the Geod.npts function that will return an array of points along the path. Note that it doesn't include the terminal points in the array, so you need to take them into account:

``````import pyproj
# calculate distance between points
g = pyproj.Geod(ellps='WGS84')
(az12, az21, dist) = g.inv(startlong, startlat, endlong, endlat)

# calculate line string along path with segments <= 1 km
lonlats = g.npts(startlong, startlat, endlong, endlat,
1 + int(dist / 1000))

# npts doesn't include start/end points, so prepend/append them
lonlats.insert(0, (startlong, startlat))
lonlats.append((endlong, endlat))
``````
• Thanks! Solution provided by well-known & massively used library here :) Commented Nov 6, 2015 at 2:14
• Great answer. I would just add that by setting the initial_idx and terminus_idx arguments within g.npts to 0, the output list will also include the initial and terminal points, which eliminates the need to manually insert them into your list. Commented May 3, 2022 at 18:18
• @KBurchfiel - my answer is very old. Pyproj may not have had that capability when I wrote this. But thanks for the update! Commented May 4, 2022 at 21:44

GeographicLib has a python interface:

This can computer geodesics on an ellipsoid (set flattening to zero to get great circles) and can generate intermediate points on a geodesic (see the "Line" commands in the sample).

Here's how to print out points on the geodesic line from JFK to Changi Airport (Singapore):

``````from geographiclib.geodesic import Geodesic
geod = Geodesic.WGS84

g = geod.Inverse(40.6, -73.8, 1.4, 104)
l = geod.Line(g['lat1'], g['lon1'], g['azi1'])
num = 15  # 15 intermediate steps

for i in range(num+1):
pos = l.Position(i * g['s12'] / num)
print(pos['lat2'], pos['lon2'])

->
(40.60, -73.8)
(49.78, -72.99)
(58.95, -71.81)
(68.09, -69.76)
(77.15, -65.01)
(85.76, -40.31)
(83.77, 80.76)
(74.92, 94.85)
...
``````
• The python port of GeographicLib is now available at pypi.python.org/pypi/geographiclib
– cffk
Commented Oct 8, 2011 at 12:28
• See also this paper: C. F. F. Karney, Algorithms for Geodesics, J. Geod, DOI: dx.doi.org/10.1007/s00190-012-0578-z
– cffk
Commented Sep 1, 2012 at 19:12

The answers provided by others are a little more elegant, but here's an ultrasimple, somewhat unpythonic, bit of Python that provides the basics. The function takes two coordinate pairs and a user-specified number of segments. It yields a set of intermediate points along a great circle path. Output: text ready to write as KML. Caveats: The code does not consider antipodes, and assumes a spherical earth.

Code by Alan Glennon http://enj.com July 2010 (author places this code in the public domain. Use at your own risk).

--

def tweensegs(longitude1,latitude1,longitude2,latitude2,num_of_segments):

``````import math

ptlon1 = longitude1
ptlat1 = latitude1
ptlon2 = longitude2
ptlat2 = latitude2

numberofsegments = num_of_segments
onelessthansegments = numberofsegments - 1
fractionalincrement = (1.0/onelessthansegments)

# 6371.009 represents the mean radius of the earth
# shortest path distance

mylats = []
mylons = []

# write the starting coordinates
mylats.append([])
mylons.append([])
mylats[0] = ptlat1
mylons[0] = ptlon1

f = fractionalincrement
icounter = 1
while (icounter <  onelessthansegments):
icountmin1 = icounter - 1
mylats.append([])
mylons.append([])
# f is expressed as a fraction along the route from point 1 to point 2
newlat=math.atan2(z,math.sqrt(math.pow(x,2)+math.pow(y,2)))
newlon=math.atan2(y,x)
newlat_degrees = math.degrees(newlat)
newlon_degrees = math.degrees(newlon)
mylats[icounter] = newlat_degrees
mylons[icounter] = newlon_degrees
icounter += 1
f = f + fractionalincrement

# write the ending coordinates
mylats.append([])
mylons.append([])
mylats[onelessthansegments] = ptlat2
mylons[onelessthansegments] = ptlon2

# Now, the array mylats[] and mylons[] have the coordinate pairs for intermediate points along the geodesic
# My mylat[0],mylat[0] and mylat[num_of_segments-1],mylat[num_of_segments-1] are the geodesic end points

# write a kml of the results
zipcounter = 0
kmlheader = "<?xml version=\"1.0\" encoding=\"UTF-8\"?><kml xmlns=\"http://www.opengis.net/kml/2.2\"><Document><name>LineString.kml</name><open>1</open><Placemark><name>unextruded</name><LineString><extrude>1</extrude><tessellate>1</tessellate><coordinates>"
while (zipcounter < numberofsegments):
outputstuff = repr(mylons[zipcounter]) + "," + repr(mylats[zipcounter]) + ",0 "
print outputstuff
zipcounter += 1
kmlfooter = "</coordinates></LineString></Placemark></Document></kml>"
print kmlfooter
``````
• note that the koders.com link is dead, redirecting to ohloh. Commented Jul 26, 2013 at 12:04

geopy A Geocoding Toolbox for Python

• +1 Didn't see functionality to do intermediate points along the geodesic, but very good toolset to know about (handles ellipsoid calcs). Thx. Commented Aug 6, 2010 at 17:09

I haven't used this package, but it seems interesting, and a possible solution: http://trac.gispython.org/lab/wiki/Shapely

• AFAIK, shapely doesn't do spheroidal calculations: `Shapely is a Python package for set-theoretic analysis and manipulation of **planar** features` Commented Jul 25, 2010 at 20:03