diff --git a/doc/references.rst b/doc/references.rst index 7ebd8eb31..013ca3182 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 `__ .. [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 `__ +.. [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 b138586bb..b56c87675 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 diff --git a/harmonica/tests/test_forward_utils.py b/harmonica/tests/test_forward_utils.py index d8986cb6a..db1d0ab21 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, + ), + )