Thread: input out of error with haversine formula
The following will return an input out of error as the acos() function cannot be -1 <= x <= 1.
SELECT * FROM
(SELECT *, (3959 * acos(cos(radians(37.7438640)) * cos(radians(37.7438640)) * cos(radians(-97.4631299) -
radians(-97.4631299)) + sin(radians(37.7438640)) * sin(radians(37.7438640))))
AS distance
FROM foo) AS distances
WHERE distance < 10
ORDER BY distance
If I break this down the following returns 1:
SELECT (cos(radians(37.7438640)) * cos(radians(37.7438640)) * cos(radians(-97.4631299) - radians(-97.4631299)) + sin(radians(37.7438640)) * sin(radians(37.743864000)));
acos(1) would give me 0.
Thoughts or workaround suggestions?
Thanks.
--Vince--
On 15 October 2010 06:03, Vince Carney <vincecarney@gmail.com> wrote: > The following will return an input out of error as the acos() function > cannot be -1 <= x <= 1. > SELECT * FROM > (SELECT *, (3959 * acos(cos(radians(37.7438640)) * > cos(radians(37.7438640)) * cos(radians(-97.4631299) - > radians(-97.4631299)) + sin(radians(37.7438640)) * > sin(radians(37.7438640)))) > AS distance > FROM foo) AS distances > WHERE distance < 10 > ORDER BY distance > If I break this down the following returns 1: > SELECT (cos(radians(37.7438640)) * cos(radians(37.7438640)) * > cos(radians(-97.4631299) - radians(-97.4631299)) + sin(radians(37.7438640)) > * sin(radians(37.743864000))); > acos(1) would give me 0. > Thoughts or workaround suggestions? > Thanks. > --Vince-- This form of the Haversine formula is known to suffer from large rounding errors when the distance between the points is small. It is much better to use the arcsin(..) form of this formula, which has much greater accuracy, particularly for the common case of small distances. See http://en.wikipedia.org/wiki/Great-circle_distance Regards, Dean
On 16 October 2010 21:13, Vince Carney <vincecarney@gmail.com> wrote: > Is it possible to do an if/else statement inside acos() to account for a > 1 > or < -1. I can't seem to get if/else working in postgres? > Try to stay on-list, so that others can benefit from the discussion. Yes you could use CASE..WHEN, but that would add extra complexity to your code without solving the underlying problem with acos(). It won't give you accurate answers when the distance between the points is small, which is usually the most common case. <begin rant> It's a shame that so many people quote this formula when telling people how to compute great circle distances, because it really isn't a very good approach. If you're prepared to consider a little maths, then think about this: Suppose y = cos(x) Then for small values of x (absolute value much less than 1), this can be approximated by a power series expansion: y = 1 - x^2 / 2 + O(x^4) So y is approximately 1, and small changes to x have a little effect on y. Thus computing x = acos(y) when y is approximately 1 is numerically very inaccurate. For small distances, computing Haversine using acos() involves an intermediate value for y, which is very close to 1. Rounding errors might push it slightly over 1, but even if they didn't, computing acos(y) is going to amplify those rounding errors, giving an inaccurate answer. For this reason, it is much better to use the Haversine formula that uses asin(). The power series for sin() looks like this: y = sin(x) = x - x^3 / 6 + O(x^5) So if x is very small, so is y, and small changes to x are matched by small changes to y, and the full accuracy of the computation is preserved. There's also another variant of the formula that uses atan2(), which is good for both small and large distances, but is a little more complicated. Aside from being slightly simpler, the acos() formula really has nothing to recommend it. <end rant> Regards, Dean http://en.wikipedia.org/wiki/Great-circle_distance