geodetic <=> geocentric added

vers

Author: scivision

pymap3d/__init__.py

@@ -37,6 +37,7 @@
 from .ned import ned2ecef, ned2geodetic, geodetic2ned, ecef2nedv, ned2aer, aer2ned, ecef2ned
 from .ecef import geodetic2ecef, ecef2geodetic, eci2geodetic, ecef2enuv, enu2ecef, ecef2enu, enu2uvw, uvw2enu
 from .sidereal import datetime2sidereal
+from .latitude import geodetic2geocentric, geocentric2geodetic
 from .lox import isometric, meridian_dist, loxodrome_inverse
 from .los import lookAtSpheroid
 from .timeconv import str2dt

pymap3d/latitude.py

@@ -0,0 +1,113 @@
+from .ellipsoid import Ellipsoid
+from math import atan, radians, degrees, tan
+
+try:
+    import numpy
+except ImportError:
+    numpy = None
+
+
+def geodetic2geocentric(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    if numpy is not None:
+        fun = numpy.vectorize(geodetic2geocentric_point)
+        return fun(geodetic_lat, ell, deg)
+    else:
+        return geodetic2geocentric_point(geodetic_lat, ell, deg)
+
+
+def geodetic2geocentric_point(geodetic_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    """
+    convert geodetic latitude to geocentric latitude.
+
+    this is like Matlab geocentricLatitude()
+    https://www.mathworks.com/help/map/ref/geocentriclatitude.html
+
+    Parameters
+    ----------
+    geodetic_lat : float
+        geodetic latitude
+    ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+    deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+    geocentric_lat : float
+         geocentric latiude
+
+    Notes
+    -----
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    """
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if abs(geodetic_lat) > 90:
+        raise ValueError("-90 <= latitude <= 90")
+
+    if deg is True:
+        geodetic_lat = radians(geodetic_lat)
+
+    geocentric_lat = atan((1 - (ell.eccentricity) ** 2) * tan(geodetic_lat))
+
+    if deg is True:
+        geocentric_lat = degrees(geocentric_lat)
+
+    return geocentric_lat
+
+
+def geocentric2geodetic(geocentric_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    if numpy is not None:
+        fun = numpy.vectorize(geocentric2geodetic_point)
+        return fun(geocentric_lat, ell, deg)
+    else:
+        return geocentric2geodetic_point(geocentric_lat, ell, deg)
+
+
+def geocentric2geodetic_point(geocentric_lat: float, ell: Ellipsoid = None, deg: bool = True) -> float:
+    """
+    converts from geocentric latitude to geodetic latitude
+
+    like Matlab geodeticLatitudeFromGeocentric
+    https://www.mathworks.com/help/map/ref/geodeticlatitudefromgeocentric.html
+
+    Parameters
+    ----------
+    geocentric_lat : float or numpy.ndarray of float
+         geocentric latitude
+    ell : Ellipsoid, optional
+         reference ellipsoid (default WGS84)
+    deg : bool, optional
+         degrees input/output  (False: radians in/out)
+
+    Returns
+    -------
+    geodetic_lat : float or numpy.ndarray of float
+         geodetic latiude
+
+    Notes
+    -----
+    Equations from J. P. Snyder, "Map Projections - A Working Manual",
+    US Geological Survey Professional Paper 1395, US Government Printing
+    Office, Washington, DC, 1987, pp. 13-18.
+    """
+
+    if ell is None:
+        ell = Ellipsoid()
+
+    if abs(geocentric_lat) > 90:
+        raise ValueError("-90 <= latitude <= 90")
+
+    if deg is True:
+        geocentric_lat = radians(geocentric_lat)
+
+    geodetic_lat = atan(tan(geocentric_lat) / (1 - (ell.eccentricity) ** 2))
+
+    if deg is True:
+        geodetic_lat = degrees(geodetic_lat)
+
+    return geodetic_lat

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pymap3d
-version = 2.0.0
+version = 2.0.1
 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_elliposid.py

@@ -6,11 +6,16 @@
 xyz0 = (660e3, -4700e3, 4247e3)
 
 
-@pytest.mark.parametrize("model,f", [("wgs84", 3.352810664747480e-03),
-                                     ("wgs72", 3.352779454167505e-03),
-                                     ("grs80", 3.352810681182319e-03),
-                                     ("clarke1866", 3.390075303928791e-03),
-                                     ("moon", 0.)])
+@pytest.mark.parametrize(
+    "model,f",
+    [
+        ("wgs84", 3.352810664747480e-03),
+        ("wgs72", 3.352779454167505e-03),
+        ("grs80", 3.352810681182319e-03),
+        ("clarke1866", 3.390075303928791e-03),
+        ("moon", 0.0),
+    ],
+)
 def test_reference(model, f):
     assert pm.Ellipsoid(model).flattening == approx(f)
 

tests/test_latitude.py

@@ -0,0 +1,33 @@
+#!/usr/bin/env python
+import pytest
+from pytest import approx
+from math import radians
+
+import pymap3d as pm
+
+
+@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.geocentric2geodetic(geocentric_lat) == approx(geodetic_lat)
+    assert pm.geocentric2geodetic(radians(geocentric_lat), 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])
+
+
+@pytest.mark.parametrize("lat", [91, -91])
+def test_badvals(lat):
+    with pytest.raises(ValueError):
+        pm.geodetic2geocentric(lat)
+        pm.geocentric2geodetic(lat)
+
+
+if __name__ == "__main__":
+    pytest.main(["-v", __file__])