WIP: Add option to use geocentric latitude with Ellipsoid.normal_gravity

[MarkWieczorek]

This PR adds the option to use geocentric latitude with Ellipsoid.normal_gravity().

The following changes were made:

1. Added the optional argument geodetic to the method normal_gravity.
2. When False, the geocentric latitude is converted to geodetic latitude using the exact same code as in the method spherical_to_geodetic. The only difference is that it is first necessary to compute the radius, and we don't care about the longitude and height parameters.
3. The documention for the routine has been modified, mostly in a subtle manner. A note was added stating that if you already have the geodetic latitude, that you should use it (otherwise there is a useless conversion).

I also made a very minor change in the doc string for geocentric_radius: "geocentric geodetic" was changed to just "geodetic" (consistent with other routines). I also note that the normal gravity doc string referred to "geometric" height, but this should be either ellipsoidal height or geodetic height: I mentioned explicitly that this is normal to the ellipsoid, as its meaning could be misinterpreted when using geocentric latitude.

I have tried my best to make the documentation as clear as possible, but there is probably room for improvement.

And to prove that it works:

In [1]: import boule

In [2]: geodetic_lat = np.linspace(-90, 90, 19)

In [3]: normal_grav = boule.WGS84.normal_gravity(geodetic_lat, height=0)

In [4]: longitude, geocentric_lat, radius = boule.WGS84.geodetic_to_spherical(0., geodetic_lat, 0.)

In [5]: normal_grav_2 = boule.WGS84.normal_gravity(geocentric_lat, height=0, geodetic=False)

In [6]: normal_grav_2 - normal_grav
Out[6]: 
array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0., 0.])

In [7]: geodetic_lat
Out[7]: 
array([-90., -80., -70., -60., -50., -40., -30., -20., -10.,   0.,  10.,
        20.,  30.,  40.,  50.,  60.,  70.,  80.,  90.])

In [8]: geocentric_lat
Out[8]: 
array([-90.        , -79.93397881, -69.87599344, -59.83307615,
       -49.81038953, -39.81061055, -29.83363581, -19.87662986,
        -9.93439421,   0.        ,   9.93439421,  19.87662986,
        29.83363581,  39.81061055,  49.81038953,  59.83307615,
        69.87599344,  79.93397881,  90.        ])

Relevant issues/PRs:
resolves #162

[MarkWieczorek]

On further reflection, I'm not sure if this is the best approach. At this point, one can enter either

• geodetic latitude and geodetic height, or
• geocentric latitude and geodetic height.
  The second might seem a bit odd.

A different way to handle this would be to allow to enter the following:

• geodetic latitude and geodetic height, or
• geocentric latitude and geocentric height (i.e., height above the ellipsoid along the vector directed to the ellipsoid center), or
• ellipsoidal-harmonic latitude and u (the semiminor axis of the ellipsoid that passes through the point).

We could then tell the normal gravity routine which coordinates we are using by specifying an optional parameter like coords or coordinates that could take values of "geodetic" (default), "geocentric", or "ellipsoidal" / "ellipsoidal-harmonic".

[leouieda]

I'm in favor of adding a coordinate_system argument instead of geodetic. The docstring can mention that the calculations are done in the geodetic system so can be more efficient to do the conversions outside the method. The options could be "geodetic", "spherical" and "ellipsoidal-harmonic" (since ellipsoidal coordinates also exist and are slightly different). I'd prefer to use spherical instead of geocentric since by definition the geodetic system is also geocentric. But I'm open to discussion if you disagree 🙂

[leouieda]

@MarkWieczorek coming back to Boule after a while and I wanted to see if we can wrap this up. How about this:

We change all gravity calculation functions take a coordinates and a coordinate_system="geodetic" argument. The coordinates is a tuple of 3 arrays. What they are depends on coordinate_system. If it's geodetic, then the 3 arrays are longitude, latitude, height. If geocentric then longitude, spherical_latitude, radius. If ellipsoidal harmonic, then longitude, beta, u. I'm proposing using the argument coordinates so we avoid a situation with latitude=None, height=None, radius=None, ....

Then we handle conversions between coordinate systems internally. For spherical, we'd have to go from spherical to geodetic to ellipsoidal harmonic (all using methods we already have). We can then add to the docstrings saying that equations are naturally in ellipsoidal harmonic so that will be the most efficient and to do conversions independently if running gravity calculations in a loop.

This would be backward incompatible since we'd change the signature of the gravity functions. To be honest, I was already a bit bothered that they didn't take longitude in the first place for consistency of the API. So I'm OK with this change breaking things.

What do you think?

cc @santisoler for another opinion.

[MarkWieczorek]

This sounds good to me. If someone wants to implement this, please feel free to do so. I don't have much free time these days...

If you wanted to retain backwards compatibility, you could allow to input either coordinates, or keywords for geocentric_lat, geodetic_lat, height, etc. This might be annoying to parse and check the validity of the input, but many people might prefer to just enter geodetic_lat, and not specify the height (assumed zero) and longitude (not used).