24

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?

6 Answers 6

12

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))
3
  • Thanks! Solution provided by well-known & massively used library here :)
    – tdihp
    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.
    – KBurchfiel
    May 3 at 18:18
  • 1
    @KBurchfiel - my answer is very old. Pyproj may not have had that capability when I wrote this. But thanks for the update!
    – scruss
    May 4 at 21:44
11

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)
...
2
  • The python port of GeographicLib is now available at pypi.python.org/pypi/geographiclib
    – cffk
    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
    Sep 1, 2012 at 19:12
8

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)

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.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
        A=math.sin((1-f)*distance_radians)/math.sin(distance_radians)
        B=math.sin(f*distance_radians)/math.sin(distance_radians)
        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=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>"
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
3

Here are some links that might help:

http://www.koders.com/python/fid0A930D7924AE856342437CA1F5A9A3EC0CAEACE2.aspx http://code.activestate.com/recipes/393241-calculating-the-distance-between-zip-codes/

1
  • note that the koders.com link is dead, redirecting to ohloh.
    – scruss
    Jul 26, 2013 at 12:04
3

geopy A Geocoding Toolbox for Python

http://code.google.com/p/geopy/wiki/GettingStarted#Calculating_distances

1
  • 1
    +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, 2010 at 17:09
-2

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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