Rhumb destinations do not wrap angles after the first pole intersection

[joe-saronic]

Given that a rhumb line has finite length to a given pole, and that RhumbDestination appears to wrap correctly in one direction, one would expect that latitude would be bounded to +/-90 and longitude to +/-180 as the path length increased. Instead, we see the following:

Here is a rust program that dumps the data:

use geo::{Point, RhumbDestination};

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

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

and some python code to load and plot:

import pandas as pd
from matplotlib import pyplot as plt

df = pd.read_csv('rhumb.csv', header=None)
fig, ax = plt.subplots(2, 1, constrained_layout=True)
ax[0].plot(df[0], df[1])
ax[1].plot(df[0], df[2])
ax[0].set_title('Lon vs Dist')
ax[1].set_title('Lat vs Dist')

The expected fix would be that latitude would be clamped, either by taking the modulo of the path length from pole to pole, or empirically.

[urschrei]

@apendleton this is your feature. I don't think this is a bug, but I agree that clamping lon and lat would be preferable. WDYT, and do you have capacity to add it?

[joe-saronic]

Returning a latitude <-90 is pretty close to being a bug :) I will take a look if I have time myself.

[apendleton]

For reference, here's a quick port of this test program to Javascript so I could test it with Turf (which is what the Rust implementation was ported from):

const fs = require('fs');
const turf = require('@turf/turf');

function main() {
    const out = [];
    const origin = [0.0, 0.0];
    for (let i = 0; i < 1000000; i++) {
        let dist = 100.0 * i;
        let dest = turf.rhumbDestination(origin, dist, 45.0, {units: 'meters'}).geometry.coordinates;
        out.push(`${dist},${dest[0]},${dest[1]}`);
    }
    fs.writeFileSync("./turf_rhumb.csv", out.join('\n'));
}

main();

And it produces broadly similar output with respect to latitude, but probably-saner output with respect to longitude:

With respect to longitude, we do already attempt to ensure that we're wrapping correctly, but it looks like there's maybe some corner-case weirdness that results in us producing some occasional out-of-bounds values as evidenced by this graph. I think wrapping (not clamping) is the preferred behavior for longitude, so probably we need to figure out what's going on that results in us occasionally producing longitudes <-180 in the Rust version.

For latitude, I'm actually not sure what we want the behavior to be for very long lines. I think clamping at +/- 90 is defensible, but it's not clear if we just allow the longitude number to continue to spin meaninglessly in that case or what. As I understand it, conceptually, for non-horizontal lines, we should be spiraling towards the pole for awhile and then eventually reach it at some fixed/calculable distance, after which further movement along the line is undefined.

I can pitch in however is most useful, but it might be good to see if any other implementations exist out in the wild to see if there's a standard/expected behavior in this circumstance? I don't have a clear frame of reference here except for Turf (which seems to also be arguably wrong).

[joe-saronic]

Returning an Option<...> which is None when you've gone too far seems like a reasonable approach. The maximum distance is proportional to the difference in latitudes and the secant of the bearing, so is not hard to calculate up-front from a given latitude. For bearing = +/- 90, the option would always be Some, with longitude determined by the local radius.

[apendleton]

An Option seems reasonable to me too, though that'd be a breaking API change, and would make rhumb destination inconsistent with haversine destination in terms of signature, so I could go either way.

As to how to calculate: yeah, makes sense. My use of geo tends to be in a pretty performance-sensitive context and trig is slow, so given that in the typical case I think we'd probably expect not to hit the pole, I might sort of be inclined to avoid doing that secant every time if it usually won't matter, vs. checking after the fact. I suppose that's empirically testable though.

[urschrei]

I would say API inconsistency isn't a big deal; it's accounting for a different detail after all. The breaking change is a concern. On the other hand, the API isn't in wide use, and the change is certainly justifiable. I can't speak to the perf concerns – we should look at the numbers first.

[apendleton]

Makes sense. I'll probably have time to dig in in more detail either Friday or this weekend, but happy to cede this to someone else if there's interest in it being done sooner.

[joe-saronic]

Turns out to be a lot simpler than I expected. Early return is actually the faster/simpler approach here since the conditional is already being checked anyway.