diff --git a/boule/ellipsoid.py b/boule/ellipsoid.py index de6e4e18..441d9d6d 100644 --- a/boule/ellipsoid.py +++ b/boule/ellipsoid.py @@ -146,6 +146,94 @@ def gravity_pole(self): ) 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 @@ def prime_vertical_radius(self, sinlat): 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 @@ def prime_vertical_radius(self, sinlat): ------- 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 diff --git a/boule/tests/test_ellipsoid.py b/boule/tests/test_ellipsoid.py index 7bd6b454..fe07cd4b 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) + )