Unify the Ellipsoids into a single class

[leouieda]

Description of the desired feature

47 talks about implementing a factory function to instantiate ellipsoids. Recently, #73 got me thinking that we do some strange things to make sure that the Sphere and Ellipsoid classes are compatible, like asking for arguments to methods that are not used.

So how about we merge everything into a single Ellipsoid class? It would take 3 axis arguments (semimajor, semimedium, semiminor) with only the semimajor being mandatory. It can then dispatch the computations appropriately depending on the number of axis given (1 goes to sphere calculations, 2 goes to ellipsoid, 3 to triaxial). The calculations themselves can be moved to functions or private methods of this class. The public facing methods would all be the same method in reality so we wouldn't have to artificially make them compatible (or forget to do so).

This would make the process of creating a new ellipsoid simpler: just use Ellipsoid and give it the appropriate amount of axis. It would also remove a lot of duplicate code and documentation.

Are you willing to help implement and maintain this feature? Yes

[leouieda]

Here is what I would imagine this class would look like.

Defining attributes:

• name
• geocentric_grav_const
• angular_velocity
• semimajor_axis
• semimedium_axis=None
• semiminor_axis=None
• flattening=None (can be given instead of semiminor_axis)
• long_name=None
• reference=None

Derived attributes:

• kind: one of ["sphere", "triaxial", "oblate" ] depending on the defined axis
• mean_radius
• volume
• gravity_pole (can be calculated independently or using normal_gravity)
• gravity_equator (can be calculated independently or using normal_gravity)
• first_eccentricity (only defined for oblate ellipsoids, others result in an exception)
• second_eccentricity (only defined for oblate ellipsoids, others result in an exception)
• linear_eccentricity (only defined for oblate ellipsoids, others result in an exception)

Public methods:

• geocentric_radius(longitude, latitude, geodetic=True) or just radius(...) or surface_radius(...) (I like the last one better): distance from center to surface. For a sphere this is the semimajor_axis. Latitude is geocentric or geodetic (only makes a difference for oblate ellipsoids).
• normal_gravity(longitude, latitude, height): the norm of the gradient of the gravity potential (gravitational + centrifugal). Latitude is geodetic
• normal_gravitation(longitude, latitude, height): the norm of the gradient of the gravitational potential
• geocentric_to_geodetic(longitude, latitude, radius): convert geocentric coordinates to geodetic coordinates. For sphere and triaxial, this changes radius to height. For oblate ellipsoids, also changes latitude from geocentric to geodetic.
• geodetic_to_geocentric(longitude, latitude, height): the inverse of the above method.

The tricky parts:

• Geodetic latitudes only really exist for oblate ellipsoids. For sphere and triaxial, we can consider this to be the same as geocentric latitudes.
• prime_vertical_radius is mostly used internally and only defined for oblate ellipsoids so it could be a private method.

The methods can detect which type of ellipsoid we have using the kind attribute and dispatch computations appropriately to private methods or functions.

To me, this has many benefits:

• The API is straight-forward and consistent across all our methods
• Having a single class makes it trivial to maintain compatibility since we have 1 of each method instead of 3
• Creating ellipsoids is easier since you only need a single class and specify the number of axis you have
• Documentation is easier since all methods take the same arguments in the same order (with tab completion). The only exception is surface_radius but it's obvious why it only takes horizontal coordinates.

The issue of users thinking that sphere gravity varies with latitude and longitude to me is not relevant. It can be resolved with a single plot. Also, we should assume that our users are trained geoscientists who know what they are doing instead of ignorant robots 🤖.

About passing unused coordinates, that is not a huge problem since in practice when do people have latitude and not longitude and height/radius? I would think that it would happen in niche cases (modelling, testing, or internally within another application) so it would be easy enough to pass None or 0. If you have all coordinate arrays, it costs nothing to pass them along.

[leouieda]

@santisoler @dabiged @MarkWieczorek @birocoles PRs #73 and #72 raised a lot of issues with the current API and I think that this would resolve all of them. Could we move some of the discussions in those PRs here instead to avoid confusion?

What do you think? I'm sure I've overlooked some edge cases or potential issues but the more I think about it the more this seems like the best way forward.

[dabiged]

I agree this is a good way forward.

This will also lower the size of the codebase considerably (dropping redundant docstrings) and should make the testing system simpler.

[santisoler]

This looks good to me! I agree with @dabiged that this would reduce the amount of docstrings and probably testing functions will be much simpler.

I cannot see any caveats or potential issues, but having a standard set of public methods and attributes, I think we can work any future issue under the hood without the need to change the public API.

I do like that users won't need to think which class to use when defining a new Ellipsoid, and having a common way to use any type of Ellipsoid will prevent a lot of headaches.

[dabiged]

I started working in this last night. At present all I have done is modify about 1/3 of the initial docstring of the ellipsoid class to encompass spheres and triaxial ellipsoids. I am not working on the actual code until I have a clear idea what we want with examples which hopefully the docstring/doctests will provide.

One thing I noticed is that it is much harder to write clear documentation with a unified class particularly around initialisation. We now need 3 examples of initialisation rather than 1 and loops over each different type of ellipsoid to show the public methods work for each type.

It is slow going but will submit a PR when I get a chance.

[dabiged]

I have just submitted #77 .

This updates the ellipsoid class to encompass spheres, triaxial ellipsoids and oblate spheroids.

I have reached an impasse with the code that is resulting in a recursionerror that I cannot seem to get past. I am asking for help resolving this. Error Summary:

• For kind = 'oblate' spheroids, the semiminor_axis method takes the semimajor_axis and flattening value to calculate the semiminor axis.
• The kind method, looks at the value of the semi-minor axis to see what type of ellipsoid it is.
• Both methods are relying on the other, and neither can work resulting in a recursion error.
• The only way I can see around this is to force the user to specify the semi-minor axis, which breaks the API. Any suggestions are welcome.

This PR is definitely the largest I have submitted and is probably against the rules. If there is a way I am supposed to break this into more manageable pieces I am open to that, I just don't know how.

[leouieda]

After much thought and discussion, we ended up deciding against doing this. The class would have to delegate to specialized functions any way so the "general" class wouldn't be very general and would just add boiler plate code for ellipsoid type checking.