Overlap and collaboration with Harmonica

[leouieda]

Hi @scivision, I've been developing a package for gravity and magnetic geophysics called Harmonica with @santisoler. We've been having the need for some coordinate conversions not found in proj4 and ended up implementing our own functions and an ellipsoid class. I recently thought of pymap3d and was poking around the code base again and there are some similarities with what we've been doing. I think there is opportunity for collaboration here and maybe porting some of our code into pymap3d. I always felt like these things were a bit out of place in Harmonica. Here are a few links:

• Our ReferenceEllipsoid class and handling functions: https://github.com/fatiando/harmonica/blob/master/harmonica/ellipsoid.py
• Conversion from geodetic to geocentric spherical: https://github.com/fatiando/harmonica/blob/master/harmonica/coordinates.py
• Using ReferenceEllipsoid to compute Normal gravity: https://github.com/fatiando/harmonica/blob/master/harmonica/gravity_corrections.py#L10

I really like our way of handling ellipsoids across the library. We have a global default ellipsoid and ways of setting different (or custom) ellipsoids using context managers. For example, you can change ellipsoids like this:

import harmonica as hm

# Set the ellipsoid globally
hm.set_ellipsoid("GRS80")

# Calculate with GRS80
normal_grav = hm.normal_gravity(latitude, height)

# Set the ellipsoid locally
with hm.set_ellipsoid("WGS84"):
    # Calculate with WGS84
    normal_grav = hm.normal_gravity(latitude, height)

# Set a custom ellipsoid
ell =  ReferenceEllipsoid(name="sphere", semimajor_axis=1, inverse_flattening=1,
                          geocentric_grav_const=10, angular_velocity=1)
with hm.set_ellipsoid(ell):
    normal_grav = hm.normal_gravity(latitude, height)

Functions that need an ellipsoid can get the currently set one by calling get_ellipsoid() and everything works. See these examples.

I really like this mechanic and I feel like it could work well with pymap3d. I don't what's the best way of moving forward but I wanted to make sure we're all aware of the two projects 🙂 Notice that we need more than a, b, f for our ellipsoids because of the gravity calculations (which could be made optional).

How committed are you to supporting the Matlab toolkit interface? If your primary goal is compatibility then it might not work well but I think there is a lot of room for improvement on the API front (many of these Matlab toolkits are a bit poorly design, IMO).

What do you think?

[scivision]

Just wanted to make an initial reply:
This sounds nice. I would want to keep the dependency optional, falling back to the simple constants I have now for Ellipsoids.
This is a good time as well to restore an original feature of pymap3d--that it have zero prereqs as a fallback. Right now it depends on Numpy, but with increasing amounts of embedded systems running Python in streaming applications, the high speed and low resource usage of math vs. numpy for streaming scalars is becoming more important.
So I'd be addressed two issues with these changes--a flexible scaleup to Harmonica and a scaledown to scalars with math only.

[scivision]

Hi @leouieda I took a quick look. From the functions you linked above, I don't initially see a conflict with the style of API I use in PyMap3D. That is, where each input/output coordinate is a separate variable, that can be scalar or N-dimensional.

The point that needs harmonization is the ability to call the PyMap3D functions as methods of the Harmonica base class. In any case, I think this object-oriented approach, where the ellipsoid is the object is nice. However, I have other users who may want to make the sensing platform the object.

What is the key thing stopping you from using PyMap3D in Harrmonica? Is it how the ellipsoid is handled? Almost all PyMap3D functions don't directly use the ellipsoid, they just pass None or the requested string value along to the core numerical functions that actually use it. So maybe if I implemented your get_ellipsoid instead of explicitly passing the ellipsoid as an argument, that would be better in general as good software engineering.

[leouieda]

Hi @scivision sorry I let the thread die here. 2019 was an eventful year for me.

Some recent progress on this front:

I've been chatting with the developers of SHTOOLS and pygeoid, which has some overlapping functionality with regards to the ellipsoids. In particular the calculation of normal gravity, not so much the coordinate conversions. See ioshchepkov/pygeoid#7

To try to pool our efforts, I create a prototype package that would house the ellipsoid definition and normal gravity calculations: Boule. The idea is that the 3 projects could share this common functionality and then implement the more specific parts in separately (shtools does spherical harmonics, pygeoid does geoid calculations, and harmonica is more applied geophysics). In the end, I ditched the get_ellipsoid mechanism since it was overly complex and would not work across several different projects.

Looking back at pymap3D, I think the boule.Ellipsoid could be passed to your functions and "just work" since you only really need the geometric parameters and they have the same name between the 2 packages. For gravity, we also need GM and angular velocity so we needed some extra functionality.

We're starting to need some simple coordinate conversions (fatiando/boule#29) and could leverage pymap3d for those:

1. Add methods to the boule.Ellipsoid that call pymap3d functions under the hood.
2. Make sure that our ellipsoid works with pymap3d (by including tests for that) and just instruct users to use these functions instead.

I'm more prone to 1 because I can never remember what things like enu2ecef mean without looking up the documentation. From what I gather, these names are inherited from the matlab package, is that right?

I understand that you want to keep the dependencies minimal (non-existent) so depending on Boule would not be an option. At the same time, I'm more interested in integrating with things like xarray and unit packages (like pint and unyt). But I'd like to leverage what you have here and hopefully contribute back to pymap3d 🙂

[scivision]

It's worth a Wiki entry, but what I'm aiming for is to have pymap3d functions with only Python stdlib i.e.math for single point conversions. This is useful for embedded systems e.g. robots that are looping and computing the next point.

Transparently, if Numpy is available then Numpy is used for N-dimensional arrays. Finally, if the input arrays are xarray, we try to pass that through too.

Yes the names are made to match Matlab Mapping toolbox, to help those looking to transition from Matlab to Python. I don't think there's any overhead to map multiple names to a single function if it's more intuitive to have a name for your specialty.

[leouieda]

Closing this since with fatiando/boule#121 we should be able to pass boule.Ellipsoid to pymap3d functions transparently. Thank you making such a good package!