Wrong earth radius in haversine_destination?

[Boscop]

WGS84 says earth radius is 6378137.0 meters (like the comment says) but the code uses 6371000.0 meters: https://georust.github.io/rust-geo/src/geo/algorithm/haversine_destination.rs.html#29

[urschrei]

IIRC, it was because

the derivation of this ellipsoidal quadratic mean radius is wrong (the averaging over azimuth is biased). When applying these examples in real applications, it is better to use the mean earth radius, 6371 km. This value is recommended by the International Union of Geodesy and Geophysics and it minimizes the RMS relative error between the great circle and geodesic distance. (emphasis mine)

See also: https://en.wikipedia.org/wiki/Earth_radius#Mean_radius, which references ftp://athena.fsv.cvut.cz/ZFG/grs80-Moritz.pdf section 4a.

[Boscop]

Ah ok. And is it correct that bearing should be in degrees, clockwise with 0° at north for haversine_destination()?

I'm trying to convert lat/lng coords into flattened cartesian coords for the purpose of vehicle routing within a confined area (less than 50km in diameter) using the first point as origin of the cartesian coord system, but somehow I get a huge error when converting back. This is my code:

pub type F = f64;
pub type Vec2 = Vector2<F>; // km

fn to_cartesian(origin: &Location, loc: Location) -> Vec2 {
    use geo::Point;
    use geo::algorithm::haversine_distance::HaversineDistance;

    let orig = Point::new(origin.longitude, origin.latitude);
    let delta_x_m = orig.haversine_distance(&Point::new(loc.longitude, origin.latitude));
    let delta_y_m = orig.haversine_distance(&Point::new(origin.longitude, loc.latitude));
    Vec2::new(delta_x_m as F * 0.001, delta_y_m as F * 0.001)
}

fn to_lat_lng(origin: &Location, loc: &Vec2) -> Location {
    use geo::Point;
    use geo::algorithm::haversine_destination::HaversineDestination;

    let origin = Point::new(origin.longitude, origin.latitude);
    let ang_deg = loc.x.atan2(loc.y).to_degrees() as f64;
    let dist_km = loc.norm() as f64;
    let p = origin.haversine_destination(ang_deg, dist_km * 1000.);
    Location { latitude: p.lat(), longitude: p.lng() }
}

(This code uses the bearing I mentioned above.) But when I convert some points to cartesian and back (with first point as origin for all)..

info!("{:?}", to_lat_lng(&origin, &to_cartesian(&origin, Location {
    latitude: 49.37502237377056, longitude: 7.049561676562462
})));
info!("{:?}", to_lat_lng(&origin, &to_cartesian(&origin, Location {
    latitude: 49.27655120272033, longitude: 7.1099476812499915
})));
info!("{:?}", to_lat_lng(&origin, &to_cartesian(&origin, Location {
    latitude: 49.33386356715103, longitude : 6.930042858203137
})));

I get this output:

Location { latitude: 49.37502237377057, longitude: 7.049561676562462 }
Location { latitude: 49.473477779951104, longitude: 7.1100689471748435 }
Location { latitude: 49.41611951259804, longitude: 7.169180569636215 }

Which is more different than it should be. Why? I'm using f64 throughout.

I'd appreciate any help on this :)

[urschrei]

Hmm I'm not sure, maybe @mkulke can shed some light here. In general, I would be wary of using Haversine for anything involving precise calculations. I'd project into a local coordinate system (I suppose State Plane and the appropriate zone for the US), and decide whether the loss of accuracy introduced by the proj4 Helmert transform (between 5m and 15m, generally) is acceptable, using a local correction file if not – it's probably fine for transport, but definitely not OK for surveying – and then use Euclidean distances, reprojecting back to WGS84 on the way out. This means an additional dependency, but rust-proj should work fine.

[Boscop]

Ah, I forgot to multiply the x and y deltas by the signum of the original distances. It seems to work with this:

fn to_cartesian(origin: &Location, loc: Location) -> Vec2 {
    use geo::Point;
    use geo::algorithm::haversine_distance::HaversineDistance;

    let orig = Point::new(origin.longitude, origin.latitude);
    let delta_x_m = orig.haversine_distance(&Point::new(loc.longitude, origin.latitude)) * (loc.longitude - origin.longitude).signum();
    let delta_y_m = orig.haversine_distance(&Point::new(origin.longitude, loc.latitude)) * (loc.latitude  - origin.latitude).signum();
    Vec2::new(delta_x_m as F * 0.001, delta_y_m as F * 0.001)
}

fn to_lat_lng(origin: &Location, loc: &Vec2) -> Location {
    use geo::Point;
    use geo::algorithm::haversine_destination::HaversineDestination;

    let origin = Point::new(origin.longitude, origin.latitude);
    let r_lng = origin.haversine_destination(90., loc.x * 1000.);
    let r_lat = origin.haversine_destination( 0., loc.y * 1000.);
    Location { latitude: r_lat.lat(), longitude: r_lng.lng() }
}

In to_lat_lng() it seems to be more accurate with two haversine destinations in x and y dir and then only taking that component from each, instead of just one haversine destination with an angle.

original:
Location { latitude: 49.37502237377056, longitude: 7.049561676562462 }
Location { latitude: 49.27655120272033, longitude: 7.1099476812499915 }
Location { latitude: 49.33386356715103, longitude : 6.930042858203137 }

separate axes to_lat_lng:
Location { latitude: 49.37502237377057, longitude: 7.049561676562462 }
Location { latitude: 49.27655120272034, longitude: 7.10994766676014 }
Location { latitude: 49.33386356715102, longitude: 6.930042970550615 }

angular to_lat_lng:
Location { latitude: 49.37502237377057, longitude: 7.049561676562462 }
Location { latitude: 49.27653551631723, longitude: 7.109826989994984 }
Location { latitude: 49.33380202783896, longitude: 6.930142948880955 }

But thanks for the idea, I'm looking into doing it the proper way like you suggested. Which coordinate system would you recommend for Europe? And does rust-proj work well on Windows? Or would this work better? https://crates.io/crates/proj5

[urschrei]

Europe is a big place, so it very much depends. For the DACH countries and their immediate neighbours, ETRS89 LAEA (epsg:3035) should work. I don't think proj5 supports conversions between ellipsoids yet, so you wouldn't be able to use it (ETRS89 uses the GRS80 datum). As long as you've installed proj.4 on Windows, rust-proj should work fine.

[mkulke]

@Boscop: i'd agree with @urschrei is right, it depends on your usecase whether haversine works or you need sth more specific/sophisticated. that said, there should not be a big delta between indirect (0deg, then 90deg) and direct (45deg) destinations with haversine_destination. I added a test-case to verify this in #177, PTAL.

[Boscop]

My usecase is vehicle routing, and I basically just need to have x,y coords so that I can calculate euclidean distances in kilometers..

[frewsxcv]

Regarding the earth radius in Haversine, hopefully https://github.com/georust/rust-geo/pull/247 clarifies things a little. Haversine treats the earth as a sphere, so it'll use the mean earth radius. The 6378137.0 number is for ellipsoidal distances, like Vincenty.

If there are any other outstanding questions here, let me know and I'll reopen, otherwise feel free to open a new issue.

[mirabilos]

I’ve researched this for the bc(1) implementation. I found that the earth is defined in WGS84 as a larger radius and a flattening, and the IUGG says to use the (2x+y)/3 mean, so the bc(1) code does:

/* WGS84 reference ellipsoid: große Halbachse (metres), Abplattung */
i = 6378137.000
x = 1/298.257223563
/* other axis */
j = i * (1 - x)
/* mean radius resulting */
r = (2 * i + j) / 3

With sufficient precision, r = 6371008.7714150598325213221998778850522660571 m.

---

Incidentally, latest data indicates that i = 6378136.600 so r = 6371008.3718621012544876529625239555347354645644595172700 metres.