GeodesicDestination produces cyclcically incosistent results

[joe-saronic]

My expectation is that repeatedly traversing a great circle with geodesic_destination will yield identical lat/lon/bearings up to floating point roundoff or so. I have made a small script to show some inconsistencies:

use geo::{GeodesicBearing, GeodesicDestination, Point};

use std::fs::OpenOptions;
use std::io::{BufWriter, Error, Write};

fn main() -> Result<(), Error> {
    let write = OpenOptions::new().write(true).create(true).open("./geodesic.csv");
    let mut writer = BufWriter::new(write.unwrap());
    let origin = Point::new(0.0, 0.0);
    writer.write("dist,lon,lat,bearing\n".as_bytes())?;
    for i in 0..1000000 {
        let dist = 100.0 * (i as f64);
        let dest = origin.geodesic_destination(45.0, dist);
        let lon = dest.x();
        let lat = dest.y();
        let bearing = (180.0 + dest.geodesic_bearing(origin)) % 360.0;
        writer.write(format!("{dist:.16},{lon:.16},{lat:.16},{bearing:.16}\n").as_bytes())?;
    }

    Ok(())
}

Here is a python script to plot the results against each other:

from matplotlib import pyplot as plt
import pandas as pd

df = pd.read_csv('./geodesic.csv')
fig, ax = plt.subplots(3, 2, constrained_layout=True)
labels = [('dist', 'lon'), ('lat', 'bearing'), ('dist', 'lat'), ('lon', 'lat'), ('dist', 'bearing'), ('lon', 'bearing')]
for a, (x, y) in zip(ax.ravel(), labels):
    a.set_xlabel(x)
    a.set_ylabel(y)
    a.plot(x, y, data=df)

The bearing may be a bit of a red herring in this case. However, the lat-lon relationship should in theory produce identical results with multiple passes, yet it does not appear to:

Zooming shows that the circle is rotating approximately one degree in longitude with every full rotation. This does not seem like something that is ascribable simply to round-off error.

[urschrei]

@phayes Could you have a look?

[phayes]

Hi @joe-saronic, thanks for the great bug report!

@urschrei, can you assign this ticket to me please? I'm not sure that I'll be able to get too deep into it today, but I'll have time next week.

[urschrei]

Thank you! Done.

[joe-saronic]

On a more careful reading of the documentation, I am beginning to think that what I was actually looking for is the Haversine destination. It looks like geodesic destination actually respects the non-spherical model, which correctly accounts for the rotation of the geodesic as it goes around the Earth.

[phayes]

For anyone else following up on this ticket, travelling around a geodesic at 45 degrees should look something like this (note that this is exaggerated in comparison to earth, earth has f = 1/298.25) :

Note that the only closed geodesic routes are where we are travelling at either 90 degrees (north-south), or 0 degrees (east-west) around the ellipsoid. Travelling these directions should result in a closed cyclic route.