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

Question

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?

score 15 · Answer 1 · answered Jul 26 '13 at 12:15

15

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

answered Jul 26 '13 at 12:15

scruss

1,989
17
23

Thanks! Solution provided by well-known & massively used library here :) – tdihp Nov 06 '15 at 02:14
1

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 03 '22 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 04 '22 at 21:44

score 11 · Answer 2 · edited Mar 06 '19 at 13:23

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

edited Mar 06 '19 at 13:23

bugmenot123

11,011
3
34
69

answered Aug 07 '11 at 19:59

cffk

3,271
18
24

The python port of GeographicLib is now available at http://pypi.python.org/pypi/geographiclib – cffk Oct 08 '11 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 01 '12 at 19:12

glennon · Answer 3 · 2010-08-01T09:24:00.760

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

score 3 · Accepted Answer · edited Aug 07 '11 at 20:59

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/

edited Aug 07 '11 at 20:59

djq

16,297
31
110
182

answered Jul 23 '10 at 00:29

B Johnson

74
1

note that the koders.com link is dead, redirecting to ohloh. – scruss Jul 26 '13 at 12:04

score 3 · Answer 5 · answered Aug 06 '10 at 14:42

3

geopy A Geocoding Toolbox for Python

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

answered Aug 06 '10 at 14:42

Googol

139
3

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 06 '10 at 17:09

score -2 · Answer 6 · answered Jul 22 '10 at 20:03

-2

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

answered Jul 22 '10 at 20:03

Hugo Estrada

296
3
4

8

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

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

6 Answers6

Linked