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 🙂