If you're using geopy, then the great_circle and vincenty distances are
equally convenient to obtain. In this case, you should almost always
use the one that gives you the more accurate result, i.e., vincenty.
The two considerations (as you point out) are speed and accuracy.
Vincenty is two times slower. But probably in a real application the
increased running time is negligible. Even if your application called
for a million distance calculations, we are only talking about a
difference in times of a couple of seconds.
For the points you use, the error in vincenty is 6 μm and the error in
the great circle distance is 0.75 m. I would then say that vincenty is
120000 times more accurate (rather than 0.17% more accurate). For general points, the error in the great circle distance can be as much as 0.5%. So can
you live with a 0.5% error in distances? For casual use (what's the
distance from Cape Town to Cairo?), probably you can. However, many GIS
applications have much stricter accuracy requirements. (0.5% is 5m over
1km. That really does make a difference.)
Nearly all serious mapping work is carried out on the reference
ellipsoid and it therefore makes sense that distances should be measured
on the ellipsoid too. Maybe you can get away with great-circle
distances today. But for each new application, you will have to check
whether this is still acceptable. Better is just to use the ellipsoidal
distance from the start. You'll sleep better at night.
ADDENDUM (May 2017)
In reply to the answer given by @craig-hicks. The vincenty() method in
geopy does have a potentially fatal flaw: it throws an error for nearly
antipodal points. The documentation in the code suggests increasing the
number of iterations. But this is not a general solution because the
iterative method used by vincenty() is unstable for such points
(each iteration takes you further from the correct solution).
Why do I characterize the problem as "potentially fatal"? Because any
use of the distance function within another software library needs to be
able to handle the exception. Handling it by returning a NaN or the
great-circle distance may not be satisfactory, because the resulting
distance function will not obey the triangle inequality which precludes
its use, e.g., in vantage-point trees.
The situation isn't completely bleak. My python package
geographiclib computes the geodesic distance accurately without any
failures. The geopy pull request #144 changes the geopy's distance
function to use geographiclib package if it's available. Unfortunately
this pull request has been in limbo since Augest 2016.
ADDENDUM (May 2018)
geopy 1.13.0 now uses the geographiclib package for computing distances.
Here's a sample call (based on the example in the original question):
>>> from geopy.distance import great_circle
>>> from geopy.distance import geodesic
>>> p1 = (31.8300167,35.0662833) # (lat, lon) - https://goo.gl/maps/TQwDd
>>> p2 = (31.8300000,35.0708167) # (lat, lon) - https://goo.gl/maps/lHrrg
>>> geodesic(p1, p2).meters
>>> great_circle(p1, p2).meters