Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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?

share|improve this question
up vote 3 down vote accepted

Here are some links that might help:

share|improve this answer
note that the link is dead, redirecting to ohloh. – scruss Jul 26 '13 at 12:04

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 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)

ptlon1_radians = math.radians(ptlon1)
ptlat1_radians = math.radians(ptlat1)
ptlon2_radians = math.radians(ptlon2)
ptlat2_radians = math.radians(ptlat2)

distance_radians=2*math.asin(math.sqrt(math.pow((math.sin((ptlat1_radians-ptlat2_radians)/2)),2) + math.cos(ptlat1_radians)*math.cos(ptlat2_radians)*math.pow((math.sin((ptlon1_radians-ptlon2_radians)/2)),2)))
# 6371.009 represents the mean radius of the earth
# shortest path distance
distance_km = 6371.009 * distance_radians

mylats = []
mylons = []

# write the starting coordinates
mylats[0] = ptlat1
mylons[0] = ptlon1 

f = fractionalincrement
icounter = 1
while (icounter <  onelessthansegments):
        icountmin1 = icounter - 1
        # f is expressed as a fraction along the route from point 1 to point 2
        x = A*math.cos(ptlat1_radians)*math.cos(ptlon1_radians) + B*math.cos(ptlat2_radians)*math.cos(ptlon2_radians)
        y = A*math.cos(ptlat1_radians)*math.sin(ptlon1_radians) +  B*math.cos(ptlat2_radians)*math.sin(ptlon2_radians)
        z = A*math.sin(ptlat1_radians) + B*math.sin(ptlat2_radians)
        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[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=\"\"><Document><name>LineString.kml</name><open>1</open><Placemark><name>unextruded</name><LineString><extrude>1</extrude><tessellate>1</tessellate><coordinates>"
print kmlheader
while (zipcounter < numberofsegments):
        outputstuff = repr(mylons[zipcounter]) + "," + repr(mylats[zipcounter]) + ",0 "
        print outputstuff
        zipcounter += 1
kmlfooter = "</coordinates></LineString></Placemark></Document></kml>"
print kmlfooter
share|improve this answer

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 intermediate points on the geodesic line from JFK to Changi Airport (Singapore):

import sys;
from GeographicLib import Geodesic, GeodesicLine
p=Geodesic.WGS84.Inverse(40.6, -73.8, 1.4, 104)
for i in range(num+1):

(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)
share|improve this answer
The python port of GeographicLib is now available at – cffk Oct 8 '11 at 12:28
See also this paper: C. F. F. Karney, Algorithms for Geodesics, J. Geod, DOI: – cffk Sep 1 '12 at 19:12

geopy A Geocoding Toolbox for Python

share|improve this answer
+1 Didn't see functionality to do intermediate points along the geodesic, but very good toolset to know about (handles ellipsoid calcs). Thx. – glennon Aug 6 '10 at 17:09

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))
share|improve this answer
Thanks! Solution provided by well-known & massively used library here :) – tdihp Nov 6 '15 at 2:14

I haven't used this package, but it seems interesting, and a possible solution:

share|improve this answer
AFAIK, shapely doesn't do spheroidal calculations: Shapely is a Python package for set-theoretic analysis and manipulation of **planar** features – fmark Jul 25 '10 at 20:03

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.