Compute distances in geodetic coordinates

[santisoler]

Description of the desired feature

Currently Harmonica has functions to compute distances between points defined in Cartesian coordinates or in spherical geocentric coordinates. Would be nice if we add a function to compute Euclidean distances between points whose coordinates are defined in geodetic (ellipsoidal) coordinates (don't confuse it with great circle distances).
Worth noting that the geodetic coordinates depend on the reference ellipsoid, so we would need to pass a boule.Ellipsoid to the new function along with the coordinates of the two points.
This may allow to compute gravitational and magnetic effects of point masses directly on geodetic coordinates without the need to convert them into spherical. Also it could enable us to build a separate equivalent layer class, a EQLHarmonicGeodetic for example, that interpolates harmonic data given in geodetic coordinates.

Some articles that may give insights of how to compute such distances are:

• Roussel et al. (2015). Complete gravity field of an ellipsoidal prism by Gauss–Legendre quadrature. doi: https://doi.org/10.1093/gji/ggv438
• Vajada et al. (2008). Global ellipsoid-referenced topographic, bathymetric and stripping corrections to gravity disturbance. doi: https://doi.org/10.1007/s11200-008-0003-5

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

[birocoles]

I can help

[leouieda]

That would be great @birocoles!

[santisoler]

I can help

Awesome! I couldn't find any closed-form expression for computing the Euclidean distance between two points defined in geocentric coordinates on the papers mentioned above. What I mean is that for computing the Euclidean distance both papers somewhat convert the coordinates to spherical and then compute the distance. @birocoles, are you aware of a closed-form expression for it?

[birocoles]

The closed form is very ugly. It can be found at Vajda et al (2004), eq 15.

Apparently, transforming coordinates to Cartesian and then compute the distance is more efficient. But it should be investigated.

[santisoler]

Thanks @birocoles for the reference.

It seems they are just converting to Cartesian (no projecting, Cartesian coordinates with the origin on the center of the Earth) and then replacing the X, Y and Z into the distance expression.

My intuition tells me this:

• Converting from geodetic to spherical and then computing distance should be less efficient that converting to Cartesian (geocentric Cartesian) and then computing the distance. I think the first approach involves the evaluation of more trigonometric functions.
• We should corroborate if the closed-form expression spare us from computing some trigonometric functions. If so, then it would definitely be more efficient than computing Cartesian coordinates and then the distance.

[santisoler]

I think the closed-form expression is a little bit more efficient.

Converting to geocentic Cartesian
If we want to convert geodetic coordinates to geocentric Cartesian and the compute the distance, we should use the following equations for both point_1 and point_2:

This means we need to evaluate 4 different trigonometric functions for each point, resulting in the evaluation of 8 trigonometric functions for computing the distance.

Closed-form formula
If we decide to use the closed-form formula we only need to evaluate 5 trigonometric functions:

So I suppose this will be more efficient than converting to Cartesian.
Do you agree?

[birocoles]

I think we need to implement both formulas and use them not only for comparing the computation times, but also the results. The results must be the same. What do you think?

[santisoler]

That's a very good idea.
So I think we should create a new distance_geodetic function that computes the Euclidean distance using the closed-form formula, and some test functions that uses the coordinates conversions to test its results. The code for performing coordinates conversion will be much easier to read, so we can safely use it to test the API function.

[santisoler]

@birocoles Would you like to tackle this (if you have time, obviously)? If not, don't worry, I can start writing the PR and we can continue the discussion over there. Again, no hurries.

[birocoles]

Ok @santisoler , I can tackle this. I will include more functions for converting coordinates in coordinates.py and then implement the distance in geodetic coordinates. I will create a first PR for this coordinate conversions.