The more general problem, posed for an ellipsoid of revolution,
is considered in Section 8 of
http://arxiv.org/abs/1109.4448v2
This gives solutions of the interception problem (the problem at
hand) and the intersection problem using the ellipsoidal
gnomonic projection. The same technique will apply to a sphere,
of course.
According to Wikipedia, Vincenty's formula is slower but more accurate:
Vincenty's formulae are two related iterative methods used in geodesy
to calculate the distance between two points on the surface of a
spheroid, developed by Thaddeus Vincenty (1975a) They are based on the
assumption that the figure of the Earth is an oblate spheroid, and
hence are more accurate than methods such as great-circle distance
which assume a spherical Earth.
The accuracy difference is ~0.17% in a 428 meters distance in Israel. I've made a quick-and-dirty speed test:
<class 'geopy.distance.vincenty'> : Total 0:00:04.125913, (0:00:00.000041 per calculation)
<class 'geopy.distance.great_circle'> : Total 0:00:02.467479, (0:00:00.000024 per calculation)
Code:
import datetime
from geopy.distance import great_circle
from geopy.distance import vincenty
p1 = (31.8300167,35.0662833)
p2 = (31.83,35.0708167)
NUM_TESTS = 100000
for strategy in vincenty, great_circle:
before = datetime.datetime.now()
for i in range(NUM_TESTS):
d=strategy(p1, p2).meters
after = datetime.datetime.now()
duration = after-before
print "%-40s: Total %s, (%s per calculation)" % (strategy, duration, duration/NUM_TESTS)
To conclude: Vincenty's formula is doubles the calculation time compared to great-circle, and its accuracy gain at the point tested is ~0.17%.
Since the calculation time is negligible, Vincenty's formula is preferred for every practical need.
Update: Following the insightful comments by whuber and cffk's and cffk's answer, I agree that the accuracy gain should be compared with the error, not the measurement. Hence, Vincenty's formula is a few orders of magnitude more accurate, not ~0.17%.
Best Answer
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: