Add reference Sphere class

[santisoler]

Add a new Sphere class to represent reference ellipsoids with zero flattening.
The Sphere class is a subclass of Ellipsoid, inheriting its methods and properties.
Some methods were overridden to avoid division by zero errors due to singularities caused by zero flattening.
Override normal_gravity method for computing the normal gravity generated by the sphere as the self gravitational field plus the centrifugal term.
Raise warning when initializing an Ellipsoid with zero (or almost zero) flattening.
Update Ellipsoid class example: replace sphere with an oblate ellipsoid.
Add test functions for the overridden and inherited methods.

Fixes #18

Reminders:

• Run make format and make check to make sure the code follows the style guide.
• Add tests for new features or tests that would have caught the bug that you're fixing.
• Add new public functions/methods/classes to doc/api/index.rst and the base __init__.py file for the package.
• Write detailed docstrings for all functions/classes/methods. It often helps to design better code if you write the docstrings first.
• If adding new functionality, add an example to the docstring, gallery, and/or tutorials.
• Add your full name, affiliation, and ORCID (optional) to the AUTHORS.md file (if you haven't already) in case you'd like to be listed as an author on the Zenodo archive of the next release.

[santisoler]

@leouieda I was thinking that maybe we would like to add some coordinate conversions capabilities to the new Sphere class.
Even though the latitudes will be the same, the geodetic height is equal to the height above the sphere surface, while the radius (in spherical coordinates) is the distance to the center of the sphere. The computations are very simple, but I think it may be nice to have it. For example, if you are working on the Moon (flattening = 0) and want to define tesseroids, you should do it in spherical coordinates. One possible way is to manually add or subtract the Moon radius, but having the methods with the same syntax that we use for Ellipsoids would make things easier. And in fact this open the possibility to reuse some workflow designed for the Earth, but on the Moon. What do you think?

[leouieda]

@santisoler do you mean adding a spherical_to_geodetic method? That seems a bit odd but we could have some coordinate conversions here for sure. Might be best for a follow up PR since it's not immediately required. What do you think?

[leouieda]

Sphere looks good @santisoler! Some comments there, particularly about checking the flattening on the ellipsoid.

[santisoler]

@leouieda After the talk we had today, I've redefined Sphere as a subclass of Ellipsoid.
I had to use attrs's post-init hooks for redefining flattening and semimajor_axis attributes and configured them to not be asked by the constructor (through init=False). I think this is the best solution, but if you know about something better, feel free to change it.

Luckily, only two methods from the original Ellipsoid were raising divisions by zero due to zero flattening.
I have also added some test functions for the inherited methods and properties.
BTW, I solved a huge issue: normal gravity was being returned in m/s^2! 😨

[leouieda]

@santisoler made a suggestion for defining the sphere attributes but after that this looks good to go. It will be good to test the gravity values against published ones when we have the Moon ellipsoid later.

[leouieda]

BTW, I solved a huge issue: normal gravity was being returned in m/s^2!

We really need to get #23 done 🙁

[leouieda]

@MarkWieczorek does this look OK to you?

[MarkWieczorek]

I just browsed through the code and have a couple of comments:

First, I think that the check in Ellipsoid to see if the flattening is close to zero is a little odd. What if you do the following? When you initialize an ellipsoid with exactly zero flattening, you return a sphere class instead. In my opinion, people will probably use Ellipsoid to initialize all cases: biaxial (which corresonds to Ellipsoid, triaxial (to implement! This might be a little complicated), and sphere.

Second, the normal gravity returns a result in units of mGals. Maybe this is ok, but I would probably prefer to have the default be SI, and then give an option to return mGals. I know that if you plot results you will always use mGals, but when you do calculations, it is usually easier just to assume that all quantities are SI and then do the conversion at the end.

Third, for the "units" suggestion: I don't have much experience with this, but there is also the astropy project that defines constants, and quantities and units. I would have to look in more detail to see what the difference is with pint (its not immediately clear if their quantity is the same as the astropy quantity).

[santisoler]

Thanks @MarkWieczorek for the feedback!

First, I think that the check in Ellipsoid to see if the flattening is close to zero is a little odd. What if you do the following? When you initialize an ellipsoid with exactly zero flattening, you return a sphere class instead. In my opinion, people will probably use Ellipsoid to initialize all cases: biaxial (which corresonds to Ellipsoid, triaxial (to implement! This might be a little complicated), and sphere.

This is tricky. We defined a new Sphere class because some of the Ellipsoid methods and properties will found divisions by zero or divisions by very small numbers, which obviously will raise problems. So, we might want to warn users if the passed flattening is very little, not only when it's equal to zero. This raises the problem of "what is a little flattening?". Because it's dimensionless I kind of feel comfortable of setting a threshold of 1e-7 for raising the warning. But setting a threshold for the constructor to decide if it should be an Ellipsoid or a Sphere doesn't sound a good idea to me.
This being said, I agree that having a single class constructor would be great. But on the other hand, only developers and users that want to define their own ellipsoids will face this issue of choosing between the Ellipsoid and theSphere class, so it doesn't seem to much of a problem. Any other user that needs to define a spherical reference ellipsoid (for the Moon, for example) would just use boule.MOON as they would use boule.WGS84, without even noticing that they are instances of different classes.

Second, the normal gravity returns a result in units of mGals. Maybe this is ok, but I would probably prefer to have the default be SI, and then give an option to return mGals. I know that if you plot results you will always use mGals, but when you do calculations, it is usually easier just to assume that all quantities are SI and then do the conversion at the end.

I totally agree to keep everything in SI. For some reason we've been breaking this rule for gravity accelerations, using mGal. This dates from Fatiando 0.5 and has been inherited by Harmonica, which returns gravity accelerations in mGal. I would keep them in mGal just for my convenience, but I cannot find a better argument, so I'm open to change it if we agree to do so. @leouieda do you have a very strong reason why we should stick with mGal?

Third, for the "units" suggestion: I don't have much experience with this, but there is also the astropy project that defines constants, and quantities and units. I would have to look in more detail to see what the difference is with pint (its not immediately clear if their quantity is the same as the astropy quantity).

Great @MarkWieczorek! We've been delaying it for quite a long time, so would be nice to implement it. BTW, Harmonica could also make a good use of it, but as Boule is a thinner package, it seems a better place for experimentation.

[leouieda]

Hi @MarkWieczorek thanks for input!

When you initialize an ellipsoid with exactly zero flattening, you return a sphere class instead.

I'm not too keen on this since we should be able to trust the users to know what they're doing. We don't actually expect end users to instantiate Ellipsoid and Sphere all that often unless they need to define a custom one. And this type of behaviour (implicitly assuming something) is usually where bugs tend to lurk in my experience.

The check is added as a warning so that it's more of a suggestion that the gravity calculations might not work. But a lot of other things would if you really wanted to make an Ellipsoid with tiny flattening. The tricky part with assuming that 0 flattening means a Sphere is floating point round off. How small do we consider to be close enough to 0?

@santisoler maybe we should remove that from this PR and leave it for when we tackle #43 instead? Then we can discuss the sanity checks all in one place.

Maybe this is ok, but I would probably prefer to have the default be SI, and then give an option to return mGals.

Yeah, I agree that maybe we should stick with SI for everything. I wouldn't want to include an option for conversion since then we go down the rabbit hole of which units do we support. It would be best to say that for now everything is SI and when we have a units package, then you can do whatever you want.

I don't have much experience with this, but there is also the astropy project that defines constants, and quantities and units.

There are 3 main projects out there for unit support:

• astropy.units which is tied with astropy. It's nice but it would mean adding astropy as a dependency just for the units, which seems a bit overkill.
• pint is a dedicated units package. I know that MetPy uses it and have had some good (and bad) experiences. But it's independent of any domain specific package and projects like xarray are looking into supporting it fully (the MetPy devs do a lot of work on that).
• unyt was born out of the astrophysics package yt. It's an independent project but not widely used outside of yt. The devs are good and it looks well designed but not that widely supported.

[MarkWieczorek]

This is tricky. We defined a new Sphere class because some of the Ellipsoid methods and properties will found divisions by zero or divisions by very small numbers, which obviously will raise problems. So, we might want to warn users if the passed flattening is very little, not only when it's equal to zero. This raises the problem of "what is a little flattening?". Because it's dimensionless I kind of feel comfortable of setting a threshold of 1e-7 for raising the warning. But setting a threshold for the constructor to decide if it should be an Ellipsoid or a Sphere doesn't sound a good idea to me.

What I would probably do is something like this when Ellispoid is called:

if f == 0:
    return Sphere( ... input parameters...)
if f < small and not quiet:
    raise warning(....)

By definition f=0.0 is a Sphere. The problem is, if f is small, should we return a sphere instead of an ellipsoid? I would say no. I would just give a warning when f is "small", and let the user proceed with caution. They could set quiet=True if they verified that this is not an issue for their use case. If the user wants to make the choice to approximate an ellipsoid by a sphere, that should be their choice. [In any case, I don't let my opinions hold this PR up!]

We don't actually expect end users to instantiate Ellipsoid and Sphere all that often unless they need to define a custom one.

For the Earth, perhaps that is true, but there are a lot of planets and moons!

As for mGals: what you could do is keep your current convention (so that nothing breaks) and then add an option like mgals=True/False.

[santisoler]

@santisoler maybe we should remove that from this PR and leave it for when we tackle #43 instead? Then we can discuss the sanity checks all in one place.

That could work, but it would depend on the design decisions we make on the Sphere class. So it might be better to keep it here for now...

[leouieda]

[In any case, I don't let my opinions hold this PR up!]

Not at all! It's easier in the long run to have these discussions before we commit to a particular API. This is a PR at its best.

By definition f=0.0 is a Sphere.

That is true. But it would be strange to instantiate a Ellipsoid class and get a Sphere instead. So I would propose raising an exception if flattening is 0 from Ellipsoid pointing people to Sphere but only issue a warning if it's small.

For the Earth, perhaps that is true, but there are a lot of planets and moons!

We'll have ALL of those 🙂

As for mGals: what you could do is keep your current convention (so that nothing breaks) and then add an option like mgals=True/False.

That is what I don't really want to do since we already plan on adding unit support. But if this is something that would be useful to have, we could add si_units=False to the normal gravity methods and then deprecate/remove that when we have units since it might take a while still.

but it would depend on the design decisions we make on the Sphere class.

@santisoler what do you mean? Those checks are independent of the Sphere class, aren't they?

[MarkWieczorek]

That is true. But it would be strange to instantiate a Ellipsoid class and get a Sphere instead. So I would propose raising an exception if flattening is 0 from Ellipsoid pointing people to Sphere but only issue a warning if it's small.

True. But if Ellipsoid was a superclass that only was used to instantiate subclasses, then it would make more sense. Eventually, I think it would be useful to have a triaxial ellipsoid class because all of the moons of our solar system are triaxial. So it might be interesting to convert the current Ellipsoid to BiAxialEllipsoid. Then you would have a structure like:

Ellipsoid

• Sphere
• BiAxialEllipsoid
• TriAxialEllipsoid

Ellipsoid would then determine which one to use based on the value of f and/or whether a, b and c are given.

[santisoler]

So I would propose raising an exception if flattening is 0 from Ellipsoid pointing people to Sphere but only issue a warning if it's small.

Ok, I think my first draft of this PR used to raise an error under zero flattening. So if we really make users to use Sphere on those cases, I agree.

@santisoler what do you mean? Those checks are independent of the Sphere class, aren't they?

Yes and no, the checks are independent of the Sphere class, but adding exceptions for zero flattening would remove the possibility of creating zero flattening ellipsoids until the Sphere class is merged. Besides, the exception would break some tests functions, which are already modified on this PR. So I think it just would be easier to keep this little check here.

BTW I opened #45 for checking parameter values, independent of the choice we make on how to handle zero flattening.

[leouieda]

Ah OK, I see your point now. This type of thing is better suited to "factory functions" which are basically what you described but in a function instead of class. The main advantage is that a class version of this requires some Python black magic while the function is pretty straight-forward.

But this can be for a later PR I guess when we have the triaxial case as well. @MarkWieczorek would you mind opening an issue for the triaxial ellipsoid? I suspect this will be a lot more challenging for gravity and coordinate conversions so it might take some digging to find the proper equations.

So I think it just would be easier to keep this little check here.

I'm fine with leaving it like that for now and then revisit in #45.

I created #46 to track the discussion about SI units.

[leouieda]

@santisoler unless you have anything else to add, I think this is good to go. We have some TODOs open in new issues to moving forward.

[santisoler]

I'm fine with leaving it like that for now and then revisit in #45.

Ok, I agree.

@santisoler unless you have anything else to add, I think this is good to go. We have some TODOs open in new issues to moving forward.

Sure, I can merge it right away.