Formulas to Calculate Geo Proximity

The Law of Cosines and the Haversine Formula will give identical results assuming a machine with infinite precision. The Haversine formula is more robust to floating point errors. However, today's machines have double precision of the order of 15 significant figures, and the law of cosines may work just fine for you. Both these formulas assume spherical earth, whereas Vicenty's iterative solution (most accurate) assumes ellipsoidal earth (in reality the earth is not even an ellipsoid - it is a geoid). Some references: http://www.movable-type.co.uk/scripts/gis-faq-5.1.html

It gets better: note the latitude to be used in the law of cosines as well as the Haversine is the geocentric latitude, which is different from geodetic latitude. For a sphere, these two are the same.

Which one is fastest to compute?

In order from fastest to slowest are: law of cosines (5 trig. calls) -> haversine (involves sqrt) -> Vicenty (have to solve this iteratively in a for loop)

Which one is most accurate?

Vicenty.

Which one is best when speed and accuracy are both considered?

If your problem domain is such that for the distances you are trying to calculate, the earth can be considered as flat, then you can work out (I am not going to give details) a formula of the form x = kx * difference in longitude, y = ky * difference in latitude. Then distance = sqrt(dxdx + dydy). If your problem domain is such that it can be solved with distance squared, then you won't have to take sqrt, and this formula will be as fast as you get possibly get. It has the added advantage that you can calculate the vector distance - x is distance in east direction, and y is distance in the north direction. Otherwise, experiment with the 3 and choose what works best in your situation.

So you want to:

• sort records by distance from p0
• select only records whose distance from p0 is less than r

The trick is that you don't exactly need to compute the great circle distance for that! You can do with any function from a pair of points to a real value that strictly grows with the great circle distance between the points. There are many such functions and some are much faster to compute than the various formulas for the exact great circle distance. One such function is the Euclidean distance in 3D. Converting latitude and longitude to a 3D point on the sphere doesn't involve inverse trigonometric functions.

Once you have x,Y,Z, you can realize that you don't actually need the distance from p0 to your point, because you can as well use the distance from the tangent plane at p0. That distance also strictly grows with the great circle distance, and is computed from X,Y,Z as a linear combination - not even a square root is needed. You just need to precompute the coefficients and the cutoff distance that corresponds to the desired great circle distance.