HaversineDistance raises DomainError about geostats.jl HOT 6 CLOSED

Thank you @basilkraft for reporting the issue. Please let me know if you are aware of another equivalent formula for the distance that is less prone to floating point errors. I will try to find the time to research more about it.

As a side note, I am in the process of migrating all the distances in GeoStats.jl to Distances.jl. This means that you will be able to use any distance defined therein.

from geostats.jl.

Thanks for the quick reaction @juliohm

Just found this:

Although this formula [haversine] is accurate for most distances on a sphere, it too suffers from rounding errors for the special (and somewhat unusual) case of antipodal points (on opposite ends of the sphere). A more complicated formula that is accurate for all distances is the following special case of the Vincenty formula for an ellipsoid with equal major and minor axes.
Wikipedia: Great-circle_distance

Check the link, you will find the Vincenty formula, which is, as written there, more complicated. You will also find the spherical law of cosines, which is more simple but has a drawback:

On computer systems with low floating-point precision, the spherical law of cosines formula can have large rounding errors if the distance is small (if the two points are a kilometer apart on the surface of the Earth, the cosine of the central angle comes out 0.99999999). For modern 64-bit floating-point numbers, the spherical law of cosines formula, given above, does not have serious rounding errors for distances larger than a few meters on the surface of the Earth.
Wikipedia: Great-circle_distance

I found this nice explanation (for R)
Great-circle distance calculations in R

It seems that the Vincenty formula is an overkill.

Maybe just add this ;)

    # haversine formula
    a = sin(Δφ/2)^2 + cos(φ₁)*cos(φ₂)*sin(Δλ/2)^2
    a > 1. && (a = 1.) # <==========================
    c = 2atan2(√a, √(1-a))

from geostats.jl.

Thank you for the links, very useful information. Could you maybe share your use case so that I can have a better sense of how many users would be impacted by such change? Are you doing global estimation (and hence antipodal points)?

I am willing to implement the change you suggested, but I am just trying to assess how often it occurs in practice. I'd rather avoid introducing a branch in the code, specially if there is an alternative solution.

from geostats.jl.

Everyone who works with global gridded data (e.g. satellite data) may be affected. A regular grid - besides some exotic case - will usually contain antipodal points.

from geostats.jl.

Thank you, I will investigate the issue a little further before patching the file. Perhaps there is a more accurate haversine implementation that we should use instead? Meanwhile, it is good to know that you can patch the code locally and proceed with your work.

from geostats.jl.

Thank you again for reporting the issue @bask0, I added a line for handling floating point errors. I will do the same with the implementation in Distances.jl

from geostats.jl.