Search code examples
c++openframeworks

GeographicLib: moving a latLon coordinate and back creates an offset


Using the following code (emplyoing GeographicLib) which is moving a coordinate and back again creates an offset to the original starting point. The difference is growing with the distance of movement and depends on the Azimuth. The same is true for using both GeodesicExact and Geodesic. What i want to achieve in the end is creating a latLon shape by moving the starting coordinate.

Is there are an exact/better way of doing this or do I miss something fundamental?

inline double distanceInMeters(const GeoCoords &_c1, const GeoCoords &_c2) {
    GeodesicExact   geod = geodWGS84(); // GeodesicExact::WGS84()
    double          meters;
    geod.Inverse(_c1.Latitude(), _c1.Longitude(),
        _c2.Latitude(), _c2.Longitude(),
        meters);
    return meters;
}

// move coord _byMeters in direction _azimuth
// inexact with horiz moves !!!
inline GeoCoords move(const GeoCoords &_coords, const double &_azimuth, const double &_byMeters) {
    GeodesicExact   geod = geodWGS84();  // GeodesicExact::WGS84()
    double          latOut, lngOut;
    geod.Direct(_coords.Latitude(), _coords.Longitude(), _azimuth, _byMeters, latOut, lngOut);
    return GeoCoords(latOut, lngOut);
}

inline void testDistanceMove() {
    GeoCoords c(12.3456789, 12.3456789);
    GeoCoords cc = c;
    double dist = 123459998.6789; // meters
    bool bHorz = true; // <-- creates some offset???
    bool bVert = true; // almost exact
    if (bHorz) cc = move(cc, Azimuth::WEST, dist); // 270.
    if (bVert) cc = move(cc, Azimuth::SOUTH, dist); // 180
    if (bHorz) cc = move(cc, Azimuth::EAST, dist);  // 90.
    if (bVert) cc = move(cc, Azimuth::NORTH, dist); // 0.

    ofLogNotice((__func__)) << "c : " << toString(c);
    ofLogNotice((__func__)) << "cc: " << toString(cc);
    double diff = distanceInMeters(c, cc);
    ofLogNotice((__func__)) << "diff: " << ofToString(diff, 12) << " m";
}

Solution

  • Simple notions of planar notions of geometry don't apply to a sphere on an ellipsoid. For example the sum of the interior angles of a quadrilateral is more than 360°. You would get approximate closure if the distance were small (order of 1km) and you're not close to a pole; however your distance is more than 3 times the circumference of the earth so all bets are off.

    ADDENDUM

    To help picture the issues, consider starting 1 meter south of the north pole and draw the 4 edges (successively west, south, east, and north) of distance 1 meter. Because 1 meter is much less than the radius of the earth, this is a planar problem. The polyline then looks like (the dashed lines are meridians)

    enter image description here

    The picture looks even more strange if you start within 1 meter of the south pole.