Undo deprecation of coordinate conversion methods

boule/_ellipsoid.py

@@ -423,14 +423,6 @@ def geodetic_to_spherical(self, longitude, latitude, height):
         """
         Convert from geodetic to geocentric spherical coordinates.
 
-        .. warning::
-
-            **This method is deprecated and will be removed in Boule v0.5.0.**
-            Please use the equivalent function in
-            `pymap3d <https://github.com/geospace-code/pymap3d/>`__ instead
-            (``pymap3d.geodetic2spherical``) which is available since version
-            2.9.0.
-
         The geodetic datum is defined by this ellipsoid. The coordinates are
         converted following [Vermeille2002]_.
 
@@ -454,39 +446,22 @@ def geodetic_to_spherical(self, longitude, latitude, height):
             system in degrees.
         radius : array
             Converted spherical radius coordinates in meters.
-
         """
-        warn(
-            "Ellipsoid.geodetic_to_spherical is deprecated and will be removed "
-            "in Boule v0.5.0. Use pymap3d.geodetic2spherical instead.",
-            FutureWarning,
-        )
-
-        latitude_rad = np.radians(latitude)
-        coslat, sinlat = np.cos(latitude_rad), np.sin(latitude_rad)
-        prime_vertical_radius = self.prime_vertical_radius(sinlat)
+        sinlat = np.sin(np.radians(latitude))
+        coslat = np.sqrt(1 - sinlat**2)
+        prime_radius = self.prime_vertical_radius(sinlat)
         # Instead of computing X and Y, we only compute the projection on the
         # XY plane: xy_projection = sqrt( X**2 + Y**2 )
-        xy_projection = (height + prime_vertical_radius) * coslat
-        z_cartesian = (
-            height + (1 - self.first_eccentricity**2) * prime_vertical_radius
-        ) * sinlat
-        radius = np.sqrt(xy_projection**2 + z_cartesian**2)
+        xy_projection = (height + prime_radius) * coslat
+        z_cartesian = (height + (1 - self.eccentricity**2) * prime_radius) * sinlat
+        radius = np.hypot(xy_projection, z_cartesian)
         spherical_latitude = np.degrees(np.arcsin(z_cartesian / radius))
         return longitude, spherical_latitude, radius
 
     def spherical_to_geodetic(self, longitude, spherical_latitude, radius):
         """
         Convert from geocentric spherical to geodetic coordinates.
 
-        .. warning::
-
-            **This method is deprecated and will be removed in Boule v0.5.0.**
-            Please use the equivalent function in
-            `pymap3d <https://github.com/geospace-code/pymap3d/>`__ instead
-            (``pymap3d.spherical2geodetic``) which is available since version
-            2.9.0.
-
         The geodetic datum is defined by this ellipsoid. The coordinates are
         converted following [Vermeille2002]_.
 
@@ -511,43 +486,24 @@ def spherical_to_geodetic(self, longitude, spherical_latitude, radius):
             degrees.
         height : array
             Converted ellipsoidal height coordinates in meters.
-
         """
-        warn(
-            "Ellipsoid.spherical_to_geodetic is deprecated and will be removed "
-            "in Boule v0.5.0. Use pymap3d.spherical2geodetic instead.",
-            FutureWarning,
-        )
-
-        spherical_latitude = np.radians(spherical_latitude)
-        k, big_z, big_d = self._spherical_to_geodetic_terms(spherical_latitude, radius)
-        latitude = np.degrees(
-            2 * np.arctan(big_z / (big_d + np.sqrt(big_d**2 + big_z**2)))
-        )
-        height = (
-            (k + self.first_eccentricity**2 - 1)
-            / k
-            * np.sqrt(big_d**2 + big_z**2)
-        )
-        return longitude, latitude, height
-
-    def _spherical_to_geodetic_terms(self, spherical_latitude, radius):
-        "Calculate intermediate terms needed for the conversion."
-        # Offload computation of these intermediate variables here to clean up
-        # the main function body
-        cos_latitude = np.cos(spherical_latitude)
-        big_z = radius * np.sin(spherical_latitude)
-        p_0 = radius**2 * cos_latitude**2 / self.semimajor_axis**2
-        q_0 = (1 - self.first_eccentricity**2) / self.semimajor_axis**2 * big_z**2
-        r_0 = (p_0 + q_0 - self.first_eccentricity**4) / 6
-        s_0 = self.first_eccentricity**4 * p_0 * q_0 / 4 / r_0**3
+        sinlat = np.sin(np.radians(spherical_latitude))
+        coslat = np.sqrt(1 - sinlat**2)
+        big_z = radius * sinlat
+        p_0 = radius**2 * coslat**2 / self.semimajor_axis**2
+        q_0 = (1 - self.eccentricity**2) / self.semimajor_axis**2 * big_z**2
+        r_0 = (p_0 + q_0 - self.eccentricity**4) / 6
+        s_0 = self.eccentricity**4 * p_0 * q_0 / 4 / r_0**3
         t_0 = np.cbrt(1 + s_0 + np.sqrt(2 * s_0 + s_0**2))
         u_0 = r_0 * (1 + t_0 + 1 / t_0)
-        v_0 = np.sqrt(u_0**2 + q_0 * self.first_eccentricity**4)
-        w_0 = self.first_eccentricity**2 * (u_0 + v_0 - q_0) / 2 / v_0
+        v_0 = np.sqrt(u_0**2 + q_0 * self.eccentricity**4)
+        w_0 = self.eccentricity**2 * (u_0 + v_0 - q_0) / 2 / v_0
         k = np.sqrt(u_0 + v_0 + w_0**2) - w_0
-        big_d = k * radius * cos_latitude / (k + self.first_eccentricity**2)
-        return k, big_z, big_d
+        big_d = k * radius * coslat / (k + self.eccentricity**2)
+        hypot_dz = np.hypot(big_d, big_z)
+        latitude = np.degrees(2 * np.arctan2(big_z, (big_d + hypot_dz)))
+        height = (k + self.eccentricity**2 - 1) / k * hypot_dz
+        return longitude, latitude, height
 
     def normal_gravity(self, latitude, height, si_units=False):
         r"""

boule/tests/test_ellipsoid.py

@@ -7,6 +7,8 @@
 """
 Test the base Ellipsoid class.
 """
+import warnings
+
 import numpy as np
 import numpy.testing as npt
 import pytest
@@ -18,20 +20,6 @@
 ELLIPSOID_NAMES = [e.name for e in ELLIPSOIDS]
 
 
-def test_coordinate_conversion_deprecations():
-    """
-    Check if warn is raised when using coordinate conversion functions
-    """
-    with pytest.warns(FutureWarning) as warning:
-        WGS84.geodetic_to_spherical(0, 0, 0)
-        assert len(warning) >= 1
-        assert "geodetic2spherical" in warning[0].message.args[0]
-    with pytest.warns(FutureWarning) as warning:
-        WGS84.spherical_to_geodetic(0, 0, 6_000_000)
-        assert len(warning) >= 1
-        assert "spherical2geodetic" in warning[0].message.args[0]
-
-
 def test_check_flattening():
     """
     Check if error/warns is raised after invalid flattening
@@ -68,7 +56,7 @@ def test_check_flattening():
             geocentric_grav_const=0,
             angular_velocity=0,
         )
-    with pytest.warns(UserWarning) as warn:
+    with warnings.catch_warnings(record=True) as warn:
         Ellipsoid(
             name="almost-zero-flattening",
             semimajor_axis=1,
@@ -105,7 +93,7 @@ def test_check_geocentric_grav_const():
     """
     Check if warn is raised after negative geocentric_grav_const
     """
-    with pytest.warns(UserWarning) as warn:
+    with warnings.catch_warnings(record=True) as warn:
         Ellipsoid(
             name="negative_gm",
             semimajor_axis=1,
@@ -346,7 +334,7 @@ def test_normal_gravity_computed_on_internal_point(ellipsoid):
     Check if warn is raised if height is negative
     """
     latitude = np.linspace(-90, 90, 100)
-    with pytest.warns(UserWarning) as warn:
+    with warnings.catch_warnings(record=True) as warn:
         ellipsoid.normal_gravity(latitude, height=-10)
         assert len(warn) >= 1
 

doc/ellipsoids.rst

@@ -10,6 +10,14 @@ All ellipsoids are instances of the :class:`~boule.Ellipsoid`,
 class documentations for a list their derived physical properties (attributes)
 and computations/transformations that they can perform (methods).
 
+.. admonition:: Help!
+    :class: hint
+
+    If an ellipsoid you need isn't in Boule yet, please `reach out
+    <https://www.fatiando.org/contact>`__ to the team and consider adding it
+    yourself. It requires no special knowledge of the code and is a great way
+    to help the project!
+
 .. jupyter-execute::
     :hide-code:
 

doc/index.rst

@@ -135,8 +135,8 @@ Boule is designed for:
 
     user_guide/normal_gravity.rst
     user_guide/gravity_disturbance.rst
-    user_guide/defining_ellipsoids.rst
     user_guide/coordinates.rst
+    user_guide/defining_ellipsoids.rst
 
 .. toctree::
     :maxdepth: 2

doc/user_guide/coordinates.rst

@@ -3,74 +3,92 @@
 Coordinate conversions
 ======================
 
-Boule's :class:`~boule.Ellipsoid` and :class:`~boule.Sphere` classes can be
-used with `pymap3d <https://github.com/geospace-code/pymap3d/>`__ for
-converting between different coordinate systems.
-While pymap3d defines some ellipsoids internally, you may want to use one from
-Boule if:
-
-* You want to be certain that the parameters used for coordinate conversions
-  and gravity calculations are consistent.
-* Need to :ref:`define your own ellipsoid <defining_ellipsoids>` either because
-  you need different parameters than the built-in ones or they aren't available
-  in either Boule or pymap3d.
+Geodetic to geocentric spherical
+--------------------------------
 
-.. admonition:: Help!
-    :class: hint
-
-    If an ellipsoid you need isn't in Boule yet, please `reach out
-    <https://www.fatiando.org/contact>`__ to the team and consider adding it
-    yourself. It requires no special knowledge of the code and is a great way
-    to help the project!
-
-Geodetic to spherical
----------------------
+The :class:`boule.Ellipsoid` class implements coordinate conversions between
+geocentric geodetic coordinates and geocentric spherical coordinates. Both are
+common in geophysical applications when dealing with spherical harmonics or
+spherical modeling of topography.
 
 The example below will show you how to convert geodetic latitude and height
 into geocentric spherical latitude and radius.
 
 .. jupyter-execute::
 
     import boule as bl
-    import pymap3d
     import numpy as np
 
-    latitude = np.linspace(-90, 90, 45)
     longitude = 40
+    latitude = np.linspace(-90, 90, 45)
     height = 481_000  # ICESat-2 orbit height in meters
 
-    latitude_sph, longitude_sph, radius = pymap3d.geodetic2spherical(
-        latitude, longitude, height, ell=bl.WGS84,
+    longitude_sph, latitude_sph, radius = bl.WGS84.geodetic_to_spherical(
+        longitude, latitude, height,
     )
-    print("Geodetic latitude:", latitude)
-    print("Spherical latitude:", latitude_sph)
-    print()
+
     print("Geodetic longitude:", longitude)
     print("Spherical longitude:", longitude_sph)
-    print()
+    print("Geodetic latitude:", latitude)
+    print("Spherical latitude:", latitude_sph)
     print("Height (m):", height)
     print("Radius (m):", radius)
 
 Notice that:
 
-1. The latitude is slightly different except for the poles and equator.
-2. The longitude is the same in both coordinates systems.
+1. The longitude is the same in both coordinates systems.
+2. The latitude is slightly different except for the poles and equator.
 3. The radius (distance from the center of the ellipsoid) varies even though
    the height is constant.
 
 .. tip::
 
     We used the WGS84 ellipsoid here but the workflow is the same for any
-    other oblate ellipsoid or sphere. Checkout :ref:`ellipsoids` for options.
+    other oblate ellipsoid. Checkout :ref:`ellipsoids` for options.
 
-Geodetic to Cartesian
----------------------
+Other conversions using pymap3d
+-------------------------------
+
+Boule's :class:`~boule.Ellipsoid` and :class:`~boule.Sphere` classes can be
+used with `pymap3d <https://github.com/geospace-code/pymap3d/>`__ for
+converting between different coordinate systems.
+While pymap3d defines some ellipsoids internally, you may want to use one from
+Boule if:
+
+* You want to be certain that the parameters used for coordinate conversions
+  and gravity calculations are consistent.
+* Need to :ref:`define your own ellipsoid <defining_ellipsoids>` either because
+  you need different parameters than the built-in ones or they aren't available
+  in either Boule or pymap3d.
+
+The example below converts between geodetic and geocentric spherical using
+``pymap3d.geodetic2spherical`` instead of
+:meth:`boule.Ellipsoid.geodetic_to_spherical` to achieve the same outcome:
+
+.. jupyter-execute::
+
+    import pymap3d
+
+    longitude = 40
+    latitude = np.linspace(-90, 90, 45)
+    height = 481_000  # ICESat-2 orbit height in meters
+
+    latitude_sph, longitude_sph, radius = pymap3d.geodetic2spherical(
+        latitude, longitude, height, ell=bl.WGS84,
+    )
+
+    print("Geodetic longitude:", longitude)
+    print("Spherical longitude:", longitude_sph)
+    print("Geodetic latitude:", latitude)
+    print("Spherical latitude:", latitude_sph)
+    print("Height (m):", height)
+    print("Radius (m):", radius)
 
 Another common coordinate conversion done in global studies is from geodetic
 latitude, longitude, and height to geocentric Cartesian X, Y, and Z.
 The example below performs this conversion for the location of the
 `Insight lander <https://en.wikipedia.org/wiki/InSight>`__ on Mars based on
-[Parker2019]_:
+[Parker2019]_ using the Martian ellipsoid defined in Boule:
 
 .. jupyter-execute::
 

doc/user_guide/defining_ellipsoids.rst

@@ -1,7 +1,7 @@
 .. _defining_ellipsoids:
 
-Defining ellipsoids
-===================
+Defining custom ellipsoids
+==========================
 
 Boule comes with a range of :ref:`built-in ellipsoids <ellipsoids>` already but
 you may want to define your own.

env/requirements-test.txt

@@ -2,4 +2,4 @@
 pytest
 pytest-cov
 coverage
-pymap3d
+pymap3d>=2.9.0

environment.yml

@@ -16,7 +16,7 @@ dependencies:
   - pytest
   - pytest-cov
   - coverage
-  - pymap3d>=2.9.0
+  - pymap3d>=2.9
   # Documentation
   - sphinx==4.5.*
   - sphinx-book-theme==0.3.*