Distance between points given in geodetic coordinates

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.

harmonica/forward/utils.py

@@ -1,37 +1,44 @@
 """
-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
 
  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 +50,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 +187,130 @@ 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):
+ """
+ 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_geodetic(point_p, point_q, ellipsoid)
+
+ """
+ # 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
+ distance = 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 distance

harmonica/tests/test_forward_utils.py

@@ -3,6 +3,7 @@
 """
 import pytest
 import numpy.testing as npt
+import boule as bl
 
 from ..forward.utils import distance, check_coordinate_system
 
@@ -18,6 +19,13 @@ 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,
+ )
 
 
 def test_distance_invalid_coordinate_system():
@@ -32,3 +40,27 @@ 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)