Repository: fatiando/harmonica

Issue:
Compute distances in geodetic coordinates
**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

<!--
Please be as detailed as you can in your description. If possible, include an
example of how you would like to use this feature (even better if it's a code
example).
-->

**Are you willing to help implement and maintain this feature?** Yes

<!--
Every feature we add is code that we will have to maintain and keep updated.
This takes a lot of effort. If you are willing to be involved in the project
and help maintain your feature, it will make it easier for us to accept it.
-->


Discussion:
santisoler: The function relies on the `boule.Ellipsoid._prime_vertical_radius()` method from https://github.com/fatiando/boule/pull/35.

Patch:
diff --git a/doc/references.rst b/doc/references.rst
index 7ebd8eb..013ca31 100644
--- a/doc/references.rst
+++ b/doc/references.rst
@@ -12,4 +12,5 @@ References
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 <https://doi.org/10.1007/s001900000116>`__
 .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2000). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy (2002) 76: 475. doi:`10.1007/s00190-002-0264-7 <https://doi.org/10.1007/s00190-002-0264-7>`__
+.. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314.
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.
diff --git a/harmonica/forward/utils.py b/harmonica/forward/utils.py
index b138586..b56c876 100644
--- a/harmonica/forward/utils.py
+++ b/harmonica/forward/utils.py
@@ -1,37 +1,47 @@
 """
-Utilities for forward modelling in Cartesian and spherical coordinates.
+Utilities for forward modelling
 """
 import numpy as np
 from numba import jit
 
 
-def distance(point_p, point_q, coordinate_system="cartesian"):
+def distance(point_p, point_q, coordinate_system="cartesian", ellipsoid=None):
     """
-    Distance between two points in Cartesian or spherical coordinates
+    Distance between two points in Cartesian, spherical or geodetic coordinates
+
+    Computes the Euclidean distance between two points given in Cartesian,
+    spherical or geodetic coordinates.
 
     Parameters
     ----------
     point_p : list or tuple or 1d-array
         List, tuple or array containing the coordinates of the first point in
         the following order: (``easting``, ``northing`` and ``upward``) if
-        given in Cartesian coordinates, or (``longitude``, ``latitude`` and
-        ``radius``) if given in a spherical geocentric coordiante system.
+        given in Cartesian coordinates, (``longitude``, ``latitude`` and
+        ``radius``) if given in a spherical geocentric coordiante system, or
+        (``longitude``, ``latitude`` and ``height``) if given in geodetic
+        coordinates.
         All ``easting``, ``northing`` and ``upward`` must be in meters.
         Both ``longitude`` and ``latitude`` must be in degrees, while
-        ``radius`` in meters.
+        ``radius``  and ``height`` in meters.
     point_q : list or tuple or 1d-array
         List, tuple or array containing the coordinates of the second point in
         the following order: (``easting``, ``northing`` and ``upward``) if
-        given in Cartesian coordinates, or (``longitude``, ``latitude`` and
-        ``radius``) if given in a spherical geocentric coordiante system.
+        given in Cartesian coordinates, (``longitude``, ``latitude`` and
+        ``radius``) if given in a spherical geocentric coordiante system, or
+        (``longitude``, ``latitude`` and ``height``) if given in geodetic
+        coordinates.
         All ``easting``, ``northing`` and ``upward`` must be in meters.
         Both ``longitude`` and ``latitude`` must be in degrees, while
-        ``radius`` in meters.
+        ``radius``  and ``height`` in meters.
     coordinate_system : str (optional)
         Coordinate system of the coordinates of the computation points and the
         point masses.
-        Available coordinates systems: ``cartesian``, ``spherical``.
-        Default ``cartesian``.
+        Available coordinates systems: ``cartesian``, ``spherical`` and
+       ``geodetic``. Default ``cartesian``.
+    ellipsoid : :class:`boule.Ellipsoid`
+        Reference ellipsoid for points coordinates. Ignored if
+        ``coordinate_system`` is not ``"geodetic"``. Default ``None``.
 
     Returns
     -------
@@ -43,11 +53,13 @@ def distance(point_p, point_q, coordinate_system="cartesian"):
         dist = distance_cartesian(point_p, point_q)
     if coordinate_system == "spherical":
         dist = distance_spherical(point_p, point_q)
+    if coordinate_system == "geodetic":
+        dist = distance_geodetic(point_p, point_q, ellipsoid)
     return dist
 
 
 def check_coordinate_system(
-    coordinate_system, valid_coord_systems=("cartesian", "spherical")
+    coordinate_system, valid_coord_systems=("cartesian", "spherical", "geodetic")
 ):
     """
     Check if the coordinate system is a valid one.
@@ -178,3 +190,132 @@ def distance_spherical_core(
     cospsi = sinphi_p * sinphi + cosphi_p * cosphi * coslambda
     dist = np.sqrt((radius - radius_p) ** 2 + 2 * radius * radius_p * (1 - cospsi))
     return dist, cospsi, coslambda
+
+
+def distance_geodetic(point_p, point_q, ellipsoid):  # pylint: disable=too-many-locals
+    """
+    Calculate the distance between two points in geodetic coordinates
+
+    Computes the Euclidean distance between two points given in geodetic
+    coordinates using the closed-form formula given by [Vajda2004]_.
+    All angles must be in degrees and height above the ellipsoid in meters.
+
+    Parameters
+    ----------
+    point_p : tuple or 1d-array
+        Tuple or array containing the coordinates of the first point in the
+        following order: (``longitude``, ``latitude`` and ``height``).
+        Both ``longitude`` and ``latitude`` must be in degrees, while
+        ``height`` in meters.
+    point_q : tuple or 1d-array
+        Tuple or array containing the coordinates of the second point in the
+        following order: (``longitude``, ``latitude`` and ``height``).
+        Both ``longitude`` and ``latitude`` must be in degrees, while
+        ``height`` in meters.
+    ellipsoid : :class:`boule.Ellipsoid`
+        Reference ellipsoid for the geodetic coordinates of points ``point_p``
+        and ``point_q``. Must be a instance of :class:`boule.Ellipsoid`.
+
+    Returns
+    -------
+    distance : float
+        Euclidean distance between ``point_p`` and ``point_q``.
+
+    Example
+    -------
+
+    >>> import boule as bl
+    >>> ellipsoid = bl.WGS84
+    >>> point_p = (-72.3, -33.3, 644)
+    >>> point_q = (-70.1, -31.6, 1024)
+    >>> distance = distance_geodetic(point_p, point_q, ellipsoid)
+    >>> print("{:.2f} m".format(distance))
+    279878.84 m
+
+    """
+    # Get coordinates of the two points
+    longitude, latitude, height = point_p[:]
+    longitude_p, latitude_p, height_p = point_q[:]
+    # Convert angles to radians
+    longitude, latitude = np.radians(longitude), np.radians(latitude)
+    longitude_p, latitude_p = np.radians(longitude_p), np.radians(latitude_p)
+    # Compute trigonometric quantities
+    cosphi = np.cos(latitude)
+    sinphi = np.sin(latitude)
+    cosphi_p = np.cos(latitude_p)
+    sinphi_p = np.sin(latitude_p)
+    coslambda = np.cos(longitude_p - longitude)
+    # Compute prime vertical radii for both points
+    prime_vertical_radius = ellipsoid.prime_vertical_radius(sinphi)
+    prime_vertical_radius_p = ellipsoid.prime_vertical_radius(sinphi_p)
+    # Compute the Euclidean distance using the close-form formula
+    return geodetic_distance_core(
+        cosphi,
+        sinphi,
+        height,
+        cosphi_p,
+        sinphi_p,
+        height_p,
+        coslambda,
+        prime_vertical_radius,
+        prime_vertical_radius_p,
+        ellipsoid.first_eccentricity ** 2,
+    )
+
+
+def geodetic_distance_core(
+    cosphi,
+    sinphi,
+    height,
+    cosphi_p,
+    sinphi_p,
+    height_p,
+    coslambda,
+    prime_vertical_radius,
+    prime_vertical_radius_p,
+    ecc_sq,
+):
+    """
+    Core computation of distance between two points in geodetic coordinates
+
+    Parameters
+    ----------
+    cosphi, sinphi : floats
+        Cosine and sine of the latitude angle for the first point
+    height : float
+        Height over ellipsoid of the first point (in meters).
+    cosphi_p, sinphi_p : floats
+        Cosine and sine of the latitude angle for the second point
+    height_p : float
+        Height over ellipsoid of the second point (in meters).
+    coslambda : float
+        Cosine of the difference between longitudes angles of both points.
+    prime_vertical_radius : float
+        Prime vertical radius for the latitude angle of the first point.
+    prime_vertical_radius_p : float
+        Prime vertical radius for the latitude angle of the second point.
+    ecc_sq : float
+        Square of ellipsoid first eccentricity.
+
+    Returns
+    -------
+    distance : float
+        Euclidean distance between both points.
+    """
+    upward_sum = prime_vertical_radius + height
+    upward_sum_p = prime_vertical_radius_p + height_p
+    dist = np.sqrt(
+        upward_sum_p ** 2 * cosphi_p ** 2
+        + upward_sum ** 2 * cosphi ** 2
+        - 2 * upward_sum * upward_sum_p * cosphi * cosphi_p * coslambda
+        + (upward_sum_p - ecc_sq * prime_vertical_radius_p) ** 2 * sinphi_p ** 2
+        + (upward_sum - ecc_sq * prime_vertical_radius) ** 2 * sinphi ** 2
+        - (
+            2
+            * (upward_sum_p - ecc_sq * prime_vertical_radius_p)
+            * (upward_sum - ecc_sq * prime_vertical_radius)
+            * sinphi
+            * sinphi_p
+        )
+    )
+    return dist


Test patch:
diff --git a/harmonica/tests/test_forward_utils.py b/harmonica/tests/test_forward_utils.py
index d8986cb..db1d0ab 100644
--- a/harmonica/tests/test_forward_utils.py
+++ b/harmonica/tests/test_forward_utils.py
@@ -2,7 +2,9 @@
 Test utils functions for forward modelling
 """
 import pytest
+import numpy as np
 import numpy.testing as npt
+import boule as bl
 
 from ..forward.utils import distance, check_coordinate_system
 
@@ -18,6 +20,14 @@ def test_distance():
     point_a = (32.3, 40.1, 1e4)
     point_b = (32.3, 40.1, 1e4 + 100)
     npt.assert_allclose(distance(point_a, point_b, coordinate_system="spherical"), 100)
+    # Geodetic coordinate system
+    point_a = (-71.3, 33.5, 1e4)
+    point_b = (-71.3, 33.5, 1e4 + 100)
+    npt.assert_allclose(
+        distance(point_a, point_b, coordinate_system="geodetic", ellipsoid=bl.WGS84),
+        100,
+        rtol=1e-6,
+    )
 
 
 def test_distance_invalid_coordinate_system():
@@ -32,3 +42,59 @@ def test_check_coordinate_system():
     "Check if invalid coordinate system is passed to _check_coordinate_system"
     with pytest.raises(ValueError):
         check_coordinate_system("this-is-not-a-valid-coordinate-system")
+
+
+def test_geodetic_distance_vs_spherical():
+    """
+    Compare geodetic distance vs spherical distance after conversion
+
+    Test if the closed-form formula for computing the Euclidean distance
+    between two points given in geodetic coordinates is equal to the same
+    distance computed by converting the points to spherical coordinates and
+    using the ``distance_spherical`` function.
+    """
+    # Initialize the WGS84 ellipsoid
+    ellipsoid = bl.WGS84
+    # Define two points in geodetic coordinates
+    point_a = (-69.3, -36.4, 405)
+    point_b = (-71.2, -33.3, 1025)
+    # Compute distance using closed-form formula
+    dist = distance(point_a, point_b, coordinate_system="geodetic", ellipsoid=ellipsoid)
+    # Convert points to spherical coordiantes
+    point_a_sph = ellipsoid.geodetic_to_spherical(*point_a)
+    point_b_sph = ellipsoid.geodetic_to_spherical(*point_b)
+    # Compute distance using these converted points
+    dist_sph = distance(point_a_sph, point_b_sph, coordinate_system="spherical")
+    npt.assert_allclose(dist, dist_sph)
+
+
+def test_geodetic_distance_equator_poles():
+    """
+    Check geodetic distance between points in equator and the poles
+
+    The Euclidean distance between a point on the equator and a point on the
+    pole, both located on the surface of the ellipsoid and with the same
+    longitude can computed through the semimajor and semiminor axis of the
+    ellipsoid.
+    """
+    # Initialize the WGS84 ellipsoid
+    ellipsoid = bl.WGS84
+    # Compute the expected distance between the two points
+    expected_distance = np.sqrt(
+        ellipsoid.semimajor_axis ** 2 + ellipsoid.semiminor_axis ** 2
+    )
+    # Compute distance for different longitudes and alternate between points on
+    # both poles
+    for pole_lat in (-90, 90):
+        for lon in np.linspace(0, 350, 36):
+            point_equator = (lon, 0, 0)
+            point_pole = (lon, pole_lat, 0)
+            npt.assert_allclose(
+                expected_distance,
+                distance(
+                    point_equator,
+                    point_pole,
+                    coordinate_system="geodetic",
+                    ellipsoid=ellipsoid,
+                ),
+            )