`inverse_haversine` precision for cardinal directions

[ngshya]

If I run
inverse_haversine((45.7597, 0.0), -300, 0.5*pi)
the output is
(45.69452167473055, -3.864098251954592)
while I expect that the first number remains the same as before since I am moving to the east direction (45.69452167473055 vs 45.7597).
Is this behavior correct?

[jdeniau]

I think it may be because pi seems to be an approximation:

from math import pi, cos

cos(0.5*pi)
// 6.123233995736766e-17 

Where the result should be 0

[jdeniau]

A possible solution is to hook for the directions directly on specific values (0, 0.5, 1 and 1.5) to do specific treatment, but it might be a little overkill, don't you think ?

[uri-rodberg]

I think if you go east or west, the latitude should be exactly the same. And if you go north or south, the longitude should be exactly the same. An error up to 1e-10 is acceptable, more than this means there is something wrong with the calculation.

[jdeniau]

I ended up on https://stackoverflow.com/questions/51425732/how-to-make-sinpi-and-cospi-2-zero that recommend to use sympy for this problem (that uses mpmath under the hood).

This seems like a huge dependency to add (both libs do 1000 more things that haversine itself !), licences might conflict, sympy is really slower, so it's not a easy decision to make, but if what's wanted in haversine is precision, then a switch might be a solution.

[jdeniau]

After a second thought, haversine will never be "precise" enough, as all calculations are made on a round sphere, whereas the earth is not a perfect sphere (+ it does not handle ).

But you are right : going on a cardinal point should not be weird, so I will push a fix for these specific cases

[jdeniau]

I tried to work on that and found out that it does not come from mathematical imprecision.
My mathematical skills are too rusty though.

I pushed my work on #59.

If I try to go west of a position for example at this line: https://github.com/mapado/haversine/pull/59/files#diff-aac46069c41f551819527065ee5717c7767fa1a3ec7fb2b4cde307d1a01e78c3R253

Then sin_lat * cos_d_r + cos_lat * sin_d_r * cos_brng is reduced to sin_lat * cos_d_r as cos_brng is 0, but the asin(sin_lat * cos_d_r) is not equal to lat at all.

I think it may be because we do not try to really go west, but instead to join the west point of the equator like in this picture (but I'm not really sure).

Maybe @CrapsJeroen that implement it can take a look ?

[CrapsJeroen]

I think if you go east or west, the latitude should be exactly the same. And if you go north or south, the longitude should be exactly the same. An error up to 1e-10 is acceptable, more than this means there is something wrong with the calculation.

From my understanding it's a limitation of the haversine function in itself, because it doesn't constantly adjust for the direction of WEST or EAST which is required at any latitude that isn't the equator.

So the more steps you take in your initial starting point, the further along the wrong direction you get and the more you move away from the latitude you started at. Imagine the extreme example of going 300km to the West from the exact North Pole, then you don't expect to still be on coordinate (90.0, 0.0). You would need to constantly adjust the direction of your movement.

What gives this away is that the closer you move to the equator the smaller the the aforementioned "error" becomes.

The reason this is not an issue for North or South is that the meridians of longitude are all of equal distance, namely the circumference of the earth.

So @jdeniau is right when he says "I think it may be because we do not try to really go west".

Sources:
Calculate distance, bearing and more between Latitude/Longitude points
Aviation Formulary V1.47 By Ed Williams

[jobh]

Fwiw: pyproj's geod.fwd() agrees with inverse_haversine, lat=45.69468 with back_azimuth=-92.8. So the calculation is consistent. Special-casing NSEW will not be consistent (consider the big jump from pi/2 to pi/2+epsilon).

[jobh]

Note: Great-circles (orthodromes), which Haversine computes, are not constant compass direction — which it would need to be to follow an non-Equator parallel. Such constant compass-direction courses are called rhumb-lines (loxodromes).

[blaylockbk]

Just ran into this issue myself and I wanted to share a visual example of the problem for anyone else who stumbles on this issue...

If you compute the inverse haversine moving EAST or WEST for successive points, your points will move closer to the equator and not along a constant latitude band as you might expect.

import cartopy.crs as ccrs
import cartopy.feature as cfeature
from haversine import inverse_haversine, Direction

ax = plt.subplot(projection=ccrs.PlateCarree())
ax.add_feature(cfeature.LAND, lw=0)
ax.add_feature(cfeature.OCEAN, lw=0)

start = (50,0)   # Starting point (lat, lon)
distance = 555   # Distance to move EAST (km)
n = 51           # Number of points

kwargs = dict(vmax=n, vmin=0, ec='k', lw=.5, zorder=1000)

ax.scatter(start[1], start[0], c=0, **kwargs)
ax.gridlines(ylocs=[start[0]], zorder=0, color='k', xlocs=None, lw=.5)

for i in range(n):
    lat, lon = inverse_haversine(start, distance, Direction.EAST, unit='km')
    ax.scatter(lon, lat, c=i+1, **kwargs)
    start = (lat, lon)
    
ax.set_global()

[jobh]

@blaylockbk If you plot just two points (n=1) 10000 km apart, and compare it with the great circle between those points, you will see that it does indeed fit the great circle that starts out due east.

By spacing the points closer, and correcting the compass course at each point, you will approximate the constant-bearing path (rhumb line), which is not a great circle. Only at the equator these two are the same.

Imagine this: You are standing close to the North Pole, and start walking in a straight line towards east. It will not be long until the Pole is behind you — your bearing is now almost due south! If you instead keep correcting towards east (constant latitude), you will walk in a small circle.

This is not a limitation of the Haversine formula, but the point of it: to measure and follow shortest-path geodesics ("straight lines" AKA great circles).

[jobh]

If what you want is the rhumb line (constant bearing), https://github.com/SpyrosMouselinos/distancly seems to give you that (I haven't tried it myself)

[blaylockbk]

Thanks, @jobh! I didn't know about the Rhumb line; this is very educational and helpful for my case.

I should have read your previous comment more carefully 😆

Note: Great-circles (orthodromes), which Haversine computes, are not constant compass direction — which it would need to be to follow an non-Equator parallel. Such constant compass-direction courses are called rhumb-lines (loxodromes).