Repository: fatiando/boule

Issue:
Provide function that returns the radius as a function of latitude
It would be very useful to be able to compute the radius of the reference ellipsoid as a function of latitude (geocentric or geodetic). 

The use case I am thinking of is the following. Most spherical harmonic models of planetary topography provide the absolute radius of the surface. In many cases, it makes more sense to remove either a reference ellipse (which could be a sphere), or to remove the geoid. In pyshtools, I would like to create a grid of the radius of the reference ellipsoid like `grid = pysh.SHGrid.from_ellipsoid()`, and for this it would be most easy to input a boule ellipsoid and to use a boule function to compute the radius as function of latitude.

Discussion:
leouieda: @MarkWieczorek just to check, this is what you had in mind right?
MarkWieczorek: That is about what I was looking for: But would it be possible to allow inputting geodetic or geocentric latitude? Perhaps you could just add a default parameter like `geodetic=True`. If False, the geocentric latitude would first be converted to geodetic.

I have almost never used geodetic coordinates, and when dealing with spherical harmonics, it is natural to use the gecentric latitude :)

leouieda: Hum, I'll need to think about this a bit since we need the radius in the first place to convert geocentric to geodetic. So it doesn't seem to be straight forward. But there should be a way.
leouieda: Alright, I worked out the equation and implemented the option for `geodetic=False` and some tests. I added a note saying that this is an exception since we generally favor solutions where the conversions are done by the user prior to input to simplify the code. But this wouldn't be an option here since there is no way to do the conversion without knowing the radius already.  

Seems like I got my math right :slightly_smiling_face: 

Patch:
diff --git a/boule/ellipsoid.py b/boule/ellipsoid.py
index de6e4e1..441d9d6 100644
--- a/boule/ellipsoid.py
+++ b/boule/ellipsoid.py
@@ -146,6 +146,94 @@ class Ellipsoid:
         )
         return result
 
+    def geocentric_radius(self, latitude, geodetic=True):
+        r"""
+        Distance from the center of the ellipsoid to its surface.
+
+        The geocentric radius and is a function of the geodetic latitude
+        :math:`\phi` and the semi-major and semi-minor axis, a and b:
+
+        .. math::
+
+            R(\phi) = \sqrt{\dfrac{
+                (a^2\cos\phi)^2 + (b^2\sin\phi)^2}{
+                (a\cos\phi)^2 + (b\sin\phi)^2 }
+            }
+
+        See https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius
+
+        The same could be achieved with
+        :meth:`boule.Ellipsoid.geodetic_to_spherical` by passing any value for
+        the longitudes and heights equal to zero. This method provides a
+        simpler and possibly faster alternative.
+
+        Alternatively, the geocentric radius can also be expressed in terms of
+        the geocentric (spherical) latitude :math:`\theta`:
+
+        .. math::
+
+            R(\theta) = \sqrt{\dfrac{1}{
+                (\frac{\cos\theta}{a})^2 + (\frac{\sin\theta}{b})^2 }
+            }
+
+        This can be useful if you already have the geocentric latitudes and
+        need the geocentric radius of the ellipsoid (for example, in spherical
+        harmonic analysis). In these cases, the coordinate conversion route is
+        not possible since we need a radius to do that in the first place.
+
+        Boule generally tries to avoid doing coordinate conversions inside
+        functions in favour of the user handling the conversions prior to
+        input. This simplifies the code and makes sure that users know
+        precisely which conversions are taking place. This method is an
+        exception since a coordinate conversion route would not be possible.
+
+        .. note::
+
+            No elevation is taken into account (the height is zero). If you
+            need the geocentric radius at a height other than zero, use
+            :meth:`boule.Ellipsoid.geodetic_to_spherical` instead.
+
+        Parameters
+        ----------
+        latitude : float or array
+            Latitude coordinates on geodetic coordinate system in degrees.
+        geodetic : bool
+            If True (default), will assume that latitudes are geodetic
+            latitudes. Otherwise, will that they are geocentric spherical
+            latitudes.
+
+        Returns
+        -------
+        geocentric_radius : float or array
+            The geocentric radius for the given latitude(s) in the same units
+            as the ellipsoid axis.
+
+        """
+        latitude_rad = np.radians(latitude)
+        coslat, sinlat = np.cos(latitude_rad), np.sin(latitude_rad)
+        # Avoid doing this in favour of having the user do the conversions when
+        # possible. It's not the case here, so we made an exception.
+        if geodetic:
+            radius = np.sqrt(
+                (
+                    (self.semimajor_axis ** 2 * coslat) ** 2
+                    + (self.semiminor_axis ** 2 * sinlat) ** 2
+                )
+                / (
+                    (self.semimajor_axis * coslat) ** 2
+                    + (self.semiminor_axis * sinlat) ** 2
+                )
+            )
+        else:
+            radius = np.sqrt(
+                1
+                / (
+                    (coslat / self.semimajor_axis) ** 2
+                    + (sinlat / self.semiminor_axis) ** 2
+                )
+            )
+        return radius
+
     def prime_vertical_radius(self, sinlat):
         r"""
         Calculate the prime vertical radius for a given geodetic latitude
@@ -156,10 +244,11 @@ class Ellipsoid:
 
             N(\phi) = \frac{a}{\sqrt{1 - e^2 \sin^2(\phi)}}
 
-        Where :math:`a` is the semimajor axis and :math:`e` is the first eccentricity.
+        Where :math:`a` is the semimajor axis and :math:`e` is the first
+        eccentricity.
 
-        This function receives the sine of the latitude as input to avoid repeated
-        computations of trigonometric functions.
+        This function receives the sine of the latitude as input to avoid
+        repeated computations of trigonometric functions.
 
         Parameters
         ----------
@@ -170,6 +259,7 @@ class Ellipsoid:
         -------
         prime_vertical_radius : float or array-like
             Prime vertical radius given in the same units as the semimajor axis
+
         """
         return self.semimajor_axis / np.sqrt(
             1 - self.first_eccentricity ** 2 * sinlat ** 2


Test patch:
diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py
index 7bd6b45..fe07cd4 100644
--- a/boule/tests/test_ellipsoid.py
+++ b/boule/tests/test_ellipsoid.py
@@ -208,3 +208,50 @@ def test_prime_vertical_radius(ellipsoid):
     )
     # Compare calculated vs expected values
     npt.assert_allclose(prime_vertical_radii, expected_pvr)
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_geocentric_radius(ellipsoid):
+    "Check against geocentric coordinate conversion results"
+    latitude = np.linspace(-80, 80, 100)
+    longitude = np.linspace(-180, 180, latitude.size)
+    height = np.zeros(latitude.size)
+    radius_conversion = ellipsoid.geodetic_to_spherical(longitude, latitude, height)[2]
+    npt.assert_allclose(radius_conversion, ellipsoid.geocentric_radius(latitude))
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_geocentric_radius_pole_equator(ellipsoid):
+    "Check against values at the pole and equator"
+    latitude = np.array([-90, 90, 0])
+    radius_true = np.array(
+        [ellipsoid.semiminor_axis, ellipsoid.semiminor_axis, ellipsoid.semimajor_axis]
+    )
+    npt.assert_allclose(radius_true, ellipsoid.geocentric_radius(latitude))
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_geocentric_radius_geocentric(ellipsoid):
+    "Check against coordinate conversion results with geocentric latitude"
+    latitude = np.linspace(-80, 80, 100)
+    longitude = np.linspace(-180, 180, latitude.size)
+    height = np.zeros(latitude.size)
+    latitude_spherical, radius_conversion = ellipsoid.geodetic_to_spherical(
+        longitude, latitude, height
+    )[1:]
+    npt.assert_allclose(
+        radius_conversion,
+        ellipsoid.geocentric_radius(latitude_spherical, geodetic=False),
+    )
+
+
+@pytest.mark.parametrize("ellipsoid", ELLIPSOIDS, ids=ELLIPSOID_NAMES)
+def test_geocentric_radius_geocentric_pole_equator(ellipsoid):
+    "Check against values at the pole and equator with geocentric latitude"
+    latitude = np.array([-90, 90, 0])
+    radius_true = np.array(
+        [ellipsoid.semiminor_axis, ellipsoid.semiminor_axis, ellipsoid.semimajor_axis]
+    )
+    npt.assert_allclose(
+        radius_true, ellipsoid.geocentric_radius(latitude, geodetic=False)
+    )