add geod2geoc, geocentric_radius

README.md

@@ -1,3 +1,5 @@
+# Python 3-D coordinate conversions
+
 [![image](https://zenodo.org/badge/DOI/10.5281/zenodo.213676.svg)](https://doi.org/10.5281/zenodo.213676)
 [![image](http://joss.theoj.org/papers/10.21105/joss.00580/status.svg)](https://doi.org/10.21105/joss.00580)
 [![astronomer](https://img.shields.io/endpoint.svg?url=https%3A%2F%2Fastronomer.ullaakut.eu%2Fshields%3Fowner%3Dscivision%26name%3Dpymap3d)](https://github.com/Ullaakut/astronomer/)
@@ -10,8 +12,6 @@
 [![image](https://img.shields.io/pypi/pyversions/pymap3d.svg)](https://pypi.python.org/pypi/pymap3d)
 [![PyPi Download stats](http://pepy.tech/badge/pymap3d)](http://pepy.tech/project/pymap3d)
 
-# Python 3-D coordinate conversions
-
 Pure Python (no prerequistes beyond Python itself) 3-D geographic coordinate conversions and geodesy.
 API similar to popular $1000 Matlab Mapping Toolbox routines for:
 
@@ -32,7 +32,6 @@ Thanks to our [contributors](./contributors.md).
 Pymap3d is compatible with Python ≥ 3.5 including PyPy.
 Numpy and AstroPy are optional; algorithms from Vallado and Meeus are used if AstroPy is not present.
 
-
 ## Install
 
 ```sh
@@ -48,6 +47,7 @@ pip install -e pymap3d
 ```
 
 One can verify Python functionality after installation by:
+
 ```sh
 pytest pymap3d -r a -v
 ```
@@ -79,7 +79,6 @@ lla = pm.aer2geodetic(*aer,*obslla)
 
 where tuple `lla` is comprised of scalar or N-D arrays `(lat,lon,alt)`.
 
-
 Example scripts are in the [examples](./examples) directory.
 
 ### Functions
@@ -97,17 +96,16 @@ converted to the desired coordinate system:
     azel2radec radec2azel
     vreckon vdist
     lookAtSpheroid
-    track2 departure meanm 
-    rcurve rsphere 
-
+    track2 departure meanm
+    rcurve rsphere
+    geod2geoc geoc2geod
 
 Additional functions:
 
 * loxodrome_inverse: rhumb line distance and azimuth between ellipsoid points (lat,lon)  akin to Matlab `distance('rh', ...)` and `azimuth('rh', ...)`
 * loxodrome_direct
 * geodetic latitude transforms to/from: parametric, authalic, isometric, and more in pymap3d.latitude
 
-
 Abbreviations:
 
 * [AER: Azimuth, Elevation, Range](https://en.wikipedia.org/wiki/Spherical_coordinate_system)
@@ -128,7 +126,6 @@ pymap3d seamlessly falls back to Python's math module if Numpy isn't present.
 To keep the code clean, only scalar data can be used without Numpy.
 As noted above, use list comprehension if you need vector data without Numpy.
 
-
 ### Caveats
 
 * Atmospheric effects neglected in all functions not invoking AstroPy.

pymap3d/__init__.py

@@ -38,6 +38,8 @@
 from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu
 from .sidereal import datetime2sidereal
 from .latitude import (
+    geod2geoc,
+    geoc2geod,
     geodetic2geocentric,
     geocentric2geodetic,
     geodetic2isometric,
@@ -51,7 +53,7 @@
     geodetic2parametric,
     parametric2geodetic,
 )
-from .rcurve import rcurve_parallel, rcurve_meridian, rcurve_transverse
+from .rcurve import geocentric_radius, rcurve_parallel, rcurve_meridian, rcurve_transverse
 from .rsphere import (
     rsphere_eqavol,
     rsphere_authalic,

pymap3d/latitude.py

@@ -3,6 +3,7 @@
 
 from .ellipsoid import Ellipsoid
 from .utils import sanitize
+from .rcurve import rcurve_transverse
 
 try:
     from numpy import radians, degrees, tan, sin, exp, pi, sqrt, inf, vectorize
@@ -27,22 +28,69 @@
     "geocentric2geodetic",
     "geodetic2authalic",
     "authalic2geodetic",
+    "geod2geoc",
+    "geoc2geod",
 ]
 
 if typing.TYPE_CHECKING:
     from numpy import ndarray
 
 
-def geodetic2geocentric(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
+def geoc2geod(geocentric_lat: "ndarray", geocentric_distance: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
     """
-    convert geodetic latitude to geocentric latitude.
+    convert geocentric latitude to geodetic latitude, consider mean sea level altitude
 
-    like Matlab geocentricLatitude()
+    like Matlab geoc2geod()
+
+    Parameters
+    ----------
+    geocentric_lat : "ndarray"
+        geocentric latitude
+    geocentric_distance: "ndarray"
+        distance from planet center, meters (NOT altitude above ground!)
+    ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+    deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+    geodetic_lat : "ndarray"
+         geodetic latiude
+
+
+    References
+    ----------
+    Long, S.A.T. "General-Altitude Transformation between Geocentric
+        and Geodetic Coordinates. Celestial Mechanics (12), 2, p. 225-230 (1975)
+        doi: 10.1007/BF01230214"
+    """
+    geocentric_lat, ell = sanitize(geocentric_lat, ell, deg)
+
+    r = geocentric_distance / ell.semimajor_axis
+
+    geodetic_lat = (
+        geocentric_lat
+        + (sin(2 * geocentric_lat) / r) * ell.flattening
+        + ((1 / r ** 2 + 1 / (4 * r)) * sin(4 * geocentric_lat)) * ell.flattening ** 2
+    )
+
+    return degrees(geodetic_lat) if deg else geodetic_lat
+
+
+def geodetic2geocentric(geodetic_lat: "ndarray", alt_m: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
+    """
+    convert geodetic latitude to geocentric latitude on spheroid surface
+
+    like Matlab geocentricLatitude() with alt_m = 0
+    like Matlab geod2geoc()
 
     Parameters
     ----------
     geodetic_lat : "ndarray"
         geodetic latitude
+    alt_m: "ndarray"
+        altitude above ellipsoid
     ell : Ellipsoid, optional
          reference ellipsoid (default WGS84)
     deg : bool, optional
@@ -60,22 +108,28 @@ def geodetic2geocentric(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: boo
     Office, Washington, DC, 1987, pp. 13-18.
     """
     geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
-
-    geocentric_lat = atan((1 - (ell.eccentricity) ** 2) * tan(geodetic_lat))
+    r = rcurve_transverse(geodetic_lat, ell, deg=False)
+    geocentric_lat = atan((1 - ell.eccentricity ** 2 * (r / (r + alt_m))) * tan(geodetic_lat))
 
     return degrees(geocentric_lat) if deg else geocentric_lat
 
 
-def geocentric2geodetic(geocentric_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
+geod2geoc = geodetic2geocentric
+
+
+def geocentric2geodetic(geocentric_lat: "ndarray", alt_m: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
     """
     converts from geocentric latitude to geodetic latitude
 
-    like Matlab geodeticLatitudeFromGeocentric()
+    like Matlab geodeticLatitudeFromGeocentric() when alt_m = 0
+    like Matlab geod2geoc() but with sea level altitude rather than planet center distance
 
     Parameters
     ----------
     geocentric_lat : "ndarray"
          geocentric latitude
+    alt_m: "ndarray"
+        altitude above ellipsoid
     ell : Ellipsoid, optional
          reference ellipsoid (default WGS84)
     deg : bool, optional
@@ -93,8 +147,8 @@ def geocentric2geodetic(geocentric_lat: "ndarray", ell: Ellipsoid = None, deg: b
     Office, Washington, DC, 1987, pp. 13-18.
     """
     geocentric_lat, ell = sanitize(geocentric_lat, ell, deg)
-
-    geodetic_lat = atan(tan(geocentric_lat) / (1 - (ell.eccentricity) ** 2))
+    r = rcurve_transverse(geocentric_lat, ell, deg=False)
+    geodetic_lat = atan(tan(geocentric_lat) / (1 - ell.eccentricity ** 2 * (r / (r + alt_m))))
 
     return degrees(geodetic_lat) if deg else geodetic_lat
 

pymap3d/rcurve.py

@@ -6,14 +6,30 @@
     from numpy import radians, sin, cos, sqrt
 except ImportError:
     from math import radians, sin, cos, sqrt
+
 from .ellipsoid import Ellipsoid
+from .utils import sanitize
 
-__all__ = ["rcurve_parallel", "rcurve_meridian", "rcurve_transverse"]
+__all__ = ["rcurve_parallel", "rcurve_meridian", "rcurve_transverse", "geocentric_radius"]
 
 if typing.TYPE_CHECKING:
     from numpy import ndarray
 
 
+def geocentric_radius(geodetic_lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
+    """
+    compute geocentric radius at geodetic latitude
+
+    https://en.wikipedia.org/wiki/Earth_radius#Geocentric_radius
+    """
+    geodetic_lat, ell = sanitize(geodetic_lat, ell, deg)
+
+    return sqrt(
+        ((ell.semimajor_axis ** 2 * cos(geodetic_lat)) ** 2 + (ell.semiminor_axis ** 2 * sin(geodetic_lat)) ** 2)
+        / ((ell.semimajor_axis * cos(geodetic_lat)) ** 2 + (ell.semiminor_axis * sin(geodetic_lat)) ** 2)
+    )
+
+
 def rcurve_parallel(lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
     """
     computes the radius of the small circle encompassing the globe at the specified latitude
@@ -42,6 +58,8 @@ def rcurve_parallel(lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) ->
 def rcurve_meridian(lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -> "ndarray":
     """computes the meridional radius of curvature for the ellipsoid
 
+    like Matlab rcurve('meridian', ...)
+
     Parameters
     ----------
     lat : "ndarray"
@@ -72,10 +90,12 @@ def rcurve_transverse(lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -
     intersecting the ellipsoid at the latitude which is
     normal to the surface of the ellipsoid
 
+    like Matlab rcurve('transverse', ...)
+
     Parameters
     ----------
     lat : "ndarray"
-        geodetic latitude (degrees)
+        latitude (degrees)
     ell : Ellipsoid, optional
           reference ellipsoid
     deg : bool, optional
@@ -87,9 +107,6 @@ def rcurve_transverse(lat: "ndarray", ell: Ellipsoid = None, deg: bool = True) -
         radius of ellipsoid
     """
 
-    if ell is None:
-        ell = Ellipsoid()
-    if deg:
-        lat = radians(lat)
+    lat, ell = sanitize(lat, ell, deg)
 
     return ell.semimajor_axis / sqrt(1 - (ell.eccentricity * sin(lat)) ** 2)

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pymap3d
-version = 2.1.1
+version = 2.2.0
 author = Michael Hirsch, Ph.D.
 author_email = [email protected]
 description = pure Python (no prereqs) coordinate conversions, following convention of several popular Matlab routines.

tests/test_aer.py

@@ -69,4 +69,4 @@ def test_aer_ned(aer, ned):
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])

tests/test_astropy.py

@@ -67,4 +67,4 @@ def test_eci_aer(useastropy):
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])

tests/test_eci.py

@@ -33,4 +33,4 @@ def test_eci_times(useastropy):
 
 
 if __name__ == "__main__":
-    pytest.main(["-x", __file__])
+    pytest.main([__file__])

tests/test_enu.py

@@ -53,4 +53,4 @@ def test_enu_ecef(enu, lla, xyz):
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])

tests/test_geodetic.py

@@ -184,4 +184,4 @@ def test_somenan():
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])

tests/test_latitude.py

@@ -6,20 +6,33 @@
 import pymap3d as pm
 
 
+@pytest.mark.parametrize(
+    "geodetic_lat,alt_m,geocentric_lat", [(0, 0, 0), (90, 0, 90), (-90, 0, -90), (45, 0, 44.80757678), (-45, 0, -44.80757678)]
+)
+def test_geodetic_alt_geocentric(geodetic_lat, alt_m, geocentric_lat):
+    assert pm.geod2geoc(geodetic_lat, alt_m) == approx(geocentric_lat)
+
+    r = pm.geocentric_radius(geodetic_lat)
+    assert pm.geoc2geod(geocentric_lat, r) == approx(geodetic_lat)
+    assert pm.geoc2geod(geocentric_lat, 1e5 + r) == approx(pm.geocentric2geodetic(geocentric_lat, 1e5 + alt_m))
+
+    assert pm.geod2geoc(geodetic_lat, 1e5 + alt_m) == approx(pm.geodetic2geocentric(geodetic_lat, 1e5 + alt_m))
+
+
 @pytest.mark.parametrize("geodetic_lat,geocentric_lat", [(0, 0), (90, 90), (-90, -90), (45, 44.80757678), (-45, -44.80757678)])
 def test_geodetic_geocentric(geodetic_lat, geocentric_lat):
 
-    assert pm.geodetic2geocentric(geodetic_lat) == approx(geocentric_lat)
-    assert pm.geodetic2geocentric(radians(geodetic_lat), deg=False) == approx(radians(geocentric_lat))
+    assert pm.geodetic2geocentric(geodetic_lat, 0) == approx(geocentric_lat)
+    assert pm.geodetic2geocentric(radians(geodetic_lat), 0, deg=False) == approx(radians(geocentric_lat))
 
-    assert pm.geocentric2geodetic(geocentric_lat) == approx(geodetic_lat)
-    assert pm.geocentric2geodetic(radians(geocentric_lat), deg=False) == approx(radians(geodetic_lat))
+    assert pm.geocentric2geodetic(geocentric_lat, 0) == approx(geodetic_lat)
+    assert pm.geocentric2geodetic(radians(geocentric_lat), 0, deg=False) == approx(radians(geodetic_lat))
 
 
 def test_numpy_geodetic_geocentric():
     pytest.importorskip("numpy")
-    assert pm.geodetic2geocentric([45, 0]) == approx([44.80757678, 0])
-    assert pm.geocentric2geodetic([44.80757678, 0]) == approx([45, 0])
+    assert pm.geodetic2geocentric([45, 0], 0) == approx([44.80757678, 0])
+    assert pm.geocentric2geodetic([44.80757678, 0], 0) == approx([45, 0])
 
 
 @pytest.mark.parametrize(
@@ -105,9 +118,9 @@ def test_numpy_geodetic_parametric():
 def test_badvals(lat):
     # geodetic_isometric is not included on purpose
     with pytest.raises(ValueError):
-        pm.geodetic2geocentric(lat)
+        pm.geodetic2geocentric(lat, 0)
     with pytest.raises(ValueError):
-        pm.geocentric2geodetic(lat)
+        pm.geocentric2geodetic(lat, 0)
     with pytest.raises(ValueError):
         pm.geodetic2conformal(lat)
     with pytest.raises(ValueError):
@@ -127,4 +140,4 @@ def test_badvals(lat):
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])

tests/test_look_spheroid.py

@@ -26,3 +26,7 @@ def test_array():
     truth = np.array([[42.00103959, lla0[1], 230.9413173], [42.00177328, -81.9995808, 282.84715651], [nan, nan, nan]])
 
     assert np.column_stack((lat, lon, sr)) == approx(truth, nan_ok=True)
+
+
+if __name__ == "__main__":
+    pytest.main([__file__])

tests/test_ned.py

@@ -52,4 +52,4 @@ def test_ned_geodetic():
 
 
 if __name__ == "__main__":
-    pytest.main(["-v", __file__])
+    pytest.main([__file__])