Обсуждение: input out of error with haversine formula

Поиск
Список
Период
Сортировка

input out of error with haversine formula

От
Vince Carney
Дата:
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--

Re: input out of error with haversine formula

От
Dean Rasheed
Дата:
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

Re: input out of error with haversine formula

От
Dean Rasheed
Дата:
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