-jg vs project. A bit confusing

[joa-quim]

See these examples

C:\v>echo 0 0 | mapproject -G45/45 -jg
0       0       6663472.79365

C:\v>echo 0 0 | mapproject -G45/45 -jg --PROJ_MEAN_RADIUS=mean
0       0       6663474.45716

C:\v>echo 0 0 | mapproject -G45/45 -jg --PROJ_MEAN_RADIUS=authalic
0       0       6663472.79365

C:\v>project -C0/0 -E45/45 -G20000 -Q
-3.18055468146e-15      3.18055468146e-15       0
45      45      6671.70311851

C:\v>echo 0 0 | mapproject -G45/45 -jg --PROJ_ELLIPSOID=Sphere
0       0       6671704.78406

So by default -jg use the authalic radius of the selected ellipsoid, whilst project uses a Sphere. Reading again the docs, the -jg behavior is documented but not so much that of project. For consistency, shouldn't it allow, and default, to -jg?

[PaulWessel]

I am not sure it can. Walter wrote this decades ago and I think it is sphere only.

[joa-quim]

It's not what I meant. Sure making project work on the ellipsoid it's not the point. But the -jg options seems odd. Either default to Sphere, or even make mapproject default to -je.

[PaulWessel]

But we have supported this behavior forever: Ellipsoid is used for mapping and we approximate geodesics with great circles using the authalic radius. Would be backwards trouble if we changed this, no?

[joa-quim]

It's like mapproject needing a small -R to avoid spherical (autalic?). I think all of these root on 20 years CPU reasons where ellipsoidal was rather expensive. But now? -j is a GMT6 option so we would not be breaking anything. Only the default change to ellipsoidal.

[PaulWessel]

This is the only place -j enters (in _gmt_initdistaz where we initialize distances calculations for programs such as map project and others):

/* Determine if -j selected a particular mode */
if (gmt_M_is_geographic (GMT, GMT_IN) && GMT->common.j.active) {    /* User specified a -j setting */
    static char *kind[5] = {"Cartesian", "Flat Earth", "Great Circle", "Geodesic", "Loxodrome"};
    GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Spherical distance calculation mode: %s.\n", kind[GMT->common.j.active]);
    if (mode != GMT_GREATCIRCLE)    /* We override a selection due to deprecated leading -|+ signs before increment or radius */
        GMT_Report (GMT->parent, GMT_MSG_VERBOSE, "Your distance mode (%s) differs from your -j option (%s) which takes precedence.\n", kind[mode], kind[GMT->common.j.active]);
    mode = GMT->common.j.mode;  /* Override with what -j said */
}

So if -j is not given then we don't go there - hence hard to switch to a default of ellipsoidal since we let map project and other modules do this via the leading + on the distances unit (-G+u+k for instance). It is cleaner to stick with the long default and add -je to change to ellipsoidal. I would be very afraid of trying to find all the places where we assume great circle initializations and try to fix that versus sticking with the current model. I think it remains true that for most applications the great circle is perfectly good (plotting in particular) but when you know you need to top precision and accuracy then you also know you need -je.

[joa-quim]

OK, but why the great circle is not over Sphere instead of derived spheres of an ellipsoid? This comes up when comparing with others (or even ourselves, see project) and we get different answers not immediately knowing why.

[PaulWessel]

Well, because an approximate geodesic is closer to the truth than spherical great circle, but at no additional cost.

[joa-quim]

Yes, it has a cost. The one of being different (nor Sphere not Ellipsoid) and users (me included when I forget about this again) being confused if they compare with different softwares.

[PaulWessel]

Point is the cost has been paid all along. Sounds like the solution for now is to spell this out more clearly then. Where should it go? You want something in project and then where?

[joa-quim]

Yes, project clearly. Here too https://docs.generic-mapping-tools.org/dev/ternary.html#unit-attributes but even mapproject -G could say right away something like Calculate distances along track (default to ellipsoid spherical authalic radius)

[seisman]

Closed by #1550.