21944858113

Committed 12 Feb 2026 11:31AM UTC coverage: 73.219%. Remained the same

Build # 21944858113

Build Type

push

github

Committed by

ioshchepkov

Commit Message

Upgrade CI

Coverage Stats

73 of 119 branches covered (61.34%)

Branch coverage included in aggregate %.

172 of 295 new or added lines in 20 files covered. (58.31%)

13 existing lines in 5 files now uncovered.

1078 of 1453 relevant lines covered (74.19%)

0.74 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.65

/src/pygeoid/coordinates/transform.py

"""This module contains functions for coordinate transformations"""

import functools as _functools

import astropy.units as u
import numpy as _np

##############################################################################
# 3D coordinates
##############################################################################


@u.quantity_input
def geodetic_to_cartesian(lat: u.deg, lon: u.deg, height: u.m, ell):
    """Convert geodetic to 3D cartesian coordinates.

    Convert geodetic coordinates (`lat`, `lon`, `height`) given w.r.t.
    ellipsoid `ell` to 3D cartesian coordinates (`x`, `y`, `z`).

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.
    """

    N = ell.prime_vertical_curvature_radius(lat)
    x = (N + height) * _np.cos(lat) * _np.cos(lon)
    y = (N + height) * _np.cos(lat) * _np.sin(lon)
    z = (N + height - N * ell.e2) * _np.sin(lat)

    return x, y, z


@_functools.partial(
    _np.vectorize, otypes=(_np.float64, _np.float64, _np.float64), excluded=[3, 4]
)
def _cartesian_to_geodetic(x, y, z, ell, degrees=True):
    """Convert 3D cartesian to geodetic coordinates.

    Parameters
    ----------
    x, y, z : float or array_like of floats
        Cartesian coordinates, in metres.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.
    degrees : bool, optional
        If True, the output geodetic latitude and longitude will be in degrees,
        otherwise radians.

    Returns
    -------
    lat, lon : float or array_like of floats
        Geodetic latitude and longitude.
    height : float or array_like of floats
        Geodetic height, in metres.

    Notes
    -----
    The algorithm of H. Vermeille is used for this transformation [1]_.

    References
    ----------
    .. [1] Vermeille, H., 2011. An analytical method to transform geocentric
    into geodetic coordinates. Journal of Geodesy, 85(2), pp.105-117.
    """
    e2 = ell.e2.value
    a = ell.a.value

    e4 = e2**2

    # Step 1
    p = (x**2 + y**2) / a**2
    q = (1 - e2) * z**2 / a**2
    r = (p + q - e4) / 6

    # Step 2 - 3
    e4pq = e4 * p * q
    t = 8 * r**3 + e4pq

    if (t > 0) or (t <= 0 and q != 0):
        if t > 0:
            li = _np.power(_np.sqrt(t) + _np.sqrt(e4pq), 1 / 3)
            u = 3 / 2 * r**2 / li**2 + 0.5 * (li + r / li) ** 2
        elif t <= 0 and q != 0:
            u_aux = (
                2 / 3 * _np.arctan2(_np.sqrt(e4pq), _np.sqrt(-t) + _np.sqrt(-8 * r**3))
            )

            u = -4 * r * _np.sin(u_aux) * _np.cos(_np.pi / 6 + u_aux)

        v = _np.sqrt(u**2 + e4 * q)
        w = e2 * (u + v - q) / (2 * v)
        k = (u + v) / (_np.sqrt(w**2 + u + v) + w)
        D = (k * _np.sqrt(x**2 + y**2)) / (k + e2)

        height = (k + e2 - 1) * _np.sqrt(D**2 + z**2) / k
        lat = 2 * _np.arctan2(z, D + _np.sqrt(D**2 + z**2))
    # Step 4
    elif q == 0 and p <= e4:
        e2p = _np.sqrt(e2 - p)
        me2 = _np.sqrt(1 - e2)

        height = -(a * me2 * e2p) / _np.sqrt(e2)
        lat = 2 * _np.arctan2(_np.sqrt(e4 - p), _np.sqrt(e2) * e2p + me2 * _np.sqrt(p))

    lon = _np.arctan2(y, x)

    if degrees:
        lat = _np.degrees(lat)
        lon = _np.degrees(lon)

    return lat, lon, height


@u.quantity_input
def cartesian_to_geodetic(x: u.m, y: u.m, z: u.m, ell):
    """Convert 3D cartesian to geodetic coordinates.

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : float or array_like of floats
        Geodetic height.

    Notes
    -----
    The algorithm of H. Vermeille is used for this transformation [1]_.

    References
    ----------
    .. [1] Vermeille, H., 2011. An analytical method to transform geocentric
    into geodetic coordinates. Journal of Geodesy, 85(2), pp.105-117.
    """

    lat, lon, height = _cartesian_to_geodetic(
        x.to("m").value, y.to("m").value, z.to("m").value, ell=ell, degrees=True
    )

    return lat * u.deg, lon * u.deg, height * u.m


@u.quantity_input
def cartesian_to_spherical(x: u.m, y: u.m, z: u.m):
    """Convert 3D cartesian to spherical coordinates.

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.

    Returns
    -------
    lat : ~astropy.units.Quantity
        Spherical latitude.
    lon : ~astropy.units.Quantity
        Spherical longitude.
    r : ~astropy.units.Quantity
        Radius.
    """

    radius = _np.sqrt(x**2 + y**2 + z**2)
    lat = _np.arctan2(z, _np.sqrt(x**2 + y**2))
    lon = _np.arctan2(y, x)

    return lat, lon, radius


@u.quantity_input
def spherical_to_cartesian(lat: u.deg, lon: u.deg, radius: u.m):
    """Convert spherical to 3D cartesian coordinates.

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Spherical latitude.
    lon : ~astropy.units.Quantity
        Spherical longitude.
    r : ~astropy.units.Quantity
        Radius.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.
    """

    x = radius * _np.cos(lat) * _np.cos(lon)
    y = radius * _np.cos(lat) * _np.sin(lon)
    z = radius * _np.sin(lat)

    return x, y, z


@u.quantity_input
def cartesian_to_ellipsoidal(x: u.m, y: u.m, z: u.m, ell):
    """Convert 3D cartesian to ellipsoidal-harmonic coordinates.

    Note that point (x, y, z) must be on or outside of the sphere with the
    radius equals to the linear eccentricity of the reference ellipsoid `ell`,
    i. e. (x**2 + y**2 + z**2) >= E**2.

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which ellipsoidal coordinates are referenced to.

    Returns
    -------
    rlat : ~astropy.units.Quantity
        Reduced latitude.
    lon : ~astropy.units.Quantity
        Longitude.
    u_ax : ~astropy.units.Quantity
        Polar axis of the ellipsoid passing through the given point.
    """

    le2 = ell.linear_eccentricity**2

    k = x**2 + y**2 + z**2 - le2

    if _np.any(k < 0):
        raise ValueError(
            "x**2 + y**2 + z**2 must be grater or equal to "
            + "the linear eccentricity of the reference ellipsoid."
        )

    u_ax = k * (0.5 + 0.5 * _np.sqrt(1 + (4 * le2 * z**2) / k**2))

    u_ax = _np.sqrt(u_ax)
    rlat = _np.arctan2(z * _np.sqrt(u_ax**2 + le2), u_ax * _np.sqrt(x**2 + y**2))

    lon = _np.arctan2(y, x)

    return rlat, lon, u_ax


@u.quantity_input
def ellipsoidal_to_cartesian(rlat: u.deg, lon: u.deg, u_ax: u.m, ell):
    """Convert ellipsoidal-harmonic coordinates to 3D cartesian coordinates.

    Parameters
    ----------
    rlat : ~astropy.units.Quantity
        Reduced latitude.
    lon : ~astropy.units.Quantity
        Longitude.
    u_ax : ~astropy.units.Quantity
        Polar axis of the ellipsoid passing through the given point.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Cartesian coordinates.
    """

    k = _np.sqrt(u_ax**2 + ell.linear_eccentricity**2)

    x = k * _np.cos(rlat) * _np.cos(lon)
    y = k * _np.cos(rlat) * _np.sin(lon)
    z = u_ax * _np.sin(rlat)

    return x, y, z


@u.quantity_input
def geodetic_to_spherical(lat: u.deg, lon: u.deg, height: u.m, ell):
    """Convert from geodetic to spherical coordinates.

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    lat, lon : ~astropy.units.Quantity
        Spherical latitude and longitude.
    r : ~astropy.units.Quantity
        Radius.
    """
    return cartesian_to_spherical(*geodetic_to_cartesian(lat, lon, height, ell=ell))


@u.quantity_input
def spherical_to_geodetic(lat: u.deg, lon: u.deg, radius: u.m, ell):
    """Convert spherical to geodetic coordinates.

    Parameters
    ----------
    lat, lon : ~astropy.units.Quantity
        Spherical latitude and longitude.
    r : ~astropy.units.Quantity
        Radius.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    """
    return cartesian_to_geodetic(*spherical_to_cartesian(lat, lon, radius), ell=ell)


@u.quantity_input
def geodetic_to_ellipsoidal(lat: u.deg, lon: u.deg, height: u.m, ell):
    """Convert from geodetic to ellipsoidal-harmonic coordinates.

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    rlat : ~astropy.units.Quantity
        Reduced latitude.
    lon : ~astropy.units.Quantity
        Longitude.
    u_ax : ~astropy.units.Quantity
        Polar axis of the ellipsoid passing through the given point.

    """
    return cartesian_to_ellipsoidal(
        *geodetic_to_cartesian(lat, lon, height, ell=ell), ell=ell
    )


@u.quantity_input
def ellipsoidal_to_geodetic(rlat: u.deg, lon: u.deg, u_ax: u.m, ell):
    """Convert from ellipsoidal-harmonic to geodetic coordinates.

    Parameters
    ----------
    rlat : ~astropy.units.Quantity
        Reduced latitude.
    lon : ~astropy.units.Quantity
        Longitude.
    u_ax : ~astropy.units.Quantity
        Polar axis of the ellipsoid passing through the given point.
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    lat, lon : ~astropy.units.Quantity
        Geodetic latitude and longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    """
    return cartesian_to_geodetic(
        *ellipsoidal_to_cartesian(rlat, lon, u_ax, ell=ell), ell=ell
    )


@u.quantity_input
def _ecef_to_enu_rotation_matrix(lat: u.deg, lon: u.deg):
    """Return ECEF to ENU rotation matrix.

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Geodetic or spherical latitude.
    lon : ~astropy.units.Quantity
        Geodetic or spherical longitude.

    """
    clat = _np.cos(lat)
    slat = _np.sin(lat)
    clon = _np.cos(lon)
    slon = _np.sin(lon)

    rotation_matrix = _np.array(
        [
            [-slon, clon, 0],
            [-slat * clon, -slat * slon, clat],
            [clat * clon, clat * slon, slat],
        ]
    )

    return rotation_matrix


@u.quantity_input
def ecef_to_enu(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell=None):
    """Convert geocentric cartesian to local cartesian coordinates.

    Convert Earth Centered Earth Fixed (ECEF) cartesian coordinates
    (`x`,`y`,`z`) to the local east-north-up (ENU) local cartesian
    coordinate system with the origin in (`lat0`, `lon0`, `height0`)
    or in (`lat0`, `lon0`, `r0`).

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Geocentric cartesian coordinates.
    origin : tuple of ~astropy.units.Quantity
        Ggeocentric (spherical) or geodetic coordinates of the origin
        (`lat0`, `lon0`, `r0`) or (`lat0`, `lon0`, `h0`).
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`, optional
        Reference ellipsoid to which geodetic coordinates are referenced to.
        Default is None, meaning spherical coordinates instead of geodetic.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Local east-north-up cartesian coordinates.
    """
    rotation_matrix = _ecef_to_enu_rotation_matrix(origin[0], origin[1])
    if ell is None:
        x0, y0, z0 = spherical_to_cartesian(*origin)
    else:
        x0, y0, z0 = geodetic_to_cartesian(*origin, ell=ell)

    out_shape = x.shape

    xyz_shifted = _np.array(
        [
            _np.asarray((x - x0).to("m").value).flatten(),
            _np.asarray((y - y0).to("m").value).flatten(),
            _np.asarray((z - z0).to("m").value).flatten(),
        ]
    )

    out = _np.dot(rotation_matrix, xyz_shifted) * u.m

    return (
        out[0].reshape(out_shape),
        out[1].reshape(out_shape),
        out[2].reshape(out_shape),
    )


@u.quantity_input
def enu_to_ecef(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell=None):
    """Convert local cartesian to geocentric cartesian coordinates.

    Convert local east-north-up (ENU) local cartesian
    coordinate system with the origin in (`lat0`, `lon0`, `height0`)
    to the Earth Centered Earth Fixed (ECEF) cartesian coordinates.

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Local east-north-uo cartesian coordinates.
    origin : tuple of ~astropy.units.Quantity
        Ggeocentric (spherical) or geodetic coordinates of the origin
        (`lat0`, `lon0`, `r0`) or (`lat0`, `lon0`, `h0`).
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.
        Default is None, meaning spherical coordinates instead of geodetic.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Geocentric cartesian coordinates, in metres.
    """
    rotation_matrix = _ecef_to_enu_rotation_matrix(origin[0], origin[1]).T

    out_shape = x.shape

    x, y, z = (
        _np.dot(
            rotation_matrix,
            _np.array(
                [
                    _np.asarray(x.to("m").value).flatten(),
                    _np.asarray(y.to("m").value).flatten(),
                    _np.asarray(z.to("m").value).flatten(),
                ]
            ),
        )
        * u.m
    )

    if ell is None:
        x0, y0, z0 = spherical_to_cartesian(*origin)
    else:
        x0, y0, z0 = geodetic_to_cartesian(*origin, ell=ell)

    return (
        (x + x0).reshape(out_shape),
        (y + y0).reshape(out_shape),
        (z + z0).reshape(out_shape),
    )


@u.quantity_input
def geodetic_to_enu(
    lat: u.deg, lon: u.deg, height: u.m, origin: tuple[u.deg, u.deg, u.m], ell
):
    """Convert geodetic coordinates to local cartesian coordinates.

    Convert geodetic coordinates
    (`lat`,`lon`,`height`) to the local east-north-up (ENU) local cartesian
    coordinate system with the origin in (`lat0`, `lon0`, `height0`).

    Parameters
    ----------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    origin : tuple of ~astropy.units.Quantity
        Geodetic coordinates of the origin (`lat0`, `lon0`, `h0`).
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    x, y, z : ~astropy.units.Quantity
        Local east-north-up cartesian coordinates.
    """
    return ecef_to_enu(
        *geodetic_to_cartesian(lat, lon, height, ell=ell), origin=origin, ell=ell
    )


@u.quantity_input
def enu_to_geodetic(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell):
    """Convert local cartesian coordinates to geodetic coordinates.

    Convert the local east-north-up (ENU) local cartesian
    coordinate system with the origin in (`lat0`, `lon0`, `height0`)
    to the geodetic coordinates.

    Parameters
    ----------
    x, y, z : ~astropy.units.Quantity
        Local east-north-uo cartesian coordinates.
    origin : tuple of ~astropy.units.Quantity
        Geodetic coordinates of the origin (`lat0`, `lon0`, `h0`).
    ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
        Reference ellipsoid to which geodetic coordinates are referenced to.

    Returns
    -------
    lat : ~astropy.units.Quantity
        Geodetic latitude.
    lon : ~astropy.units.Quantity
        Geodetic longitude.
    height : ~astropy.units.Quantity
        Geodetic height.
    """
    return cartesian_to_geodetic(*enu_to_ecef(x, y, z, origin=origin, ell=ell), ell=ell)


##############################################################################
# 2D coordinates
##############################################################################


@u.quantity_input
def polar_to_cartesian(theta: u.deg, radius: u.m):
    """Convert polar coordinates to 2D cartesian.

    Parameters
    ----------
    theta : ~astropy.units.Quantity
        Polar angle.
    radius : ~astropy.units.Quantity
        Radius.

    Returns
    -------
    x, y : ~astropy.units.Quantity
        Local east-north-uo cartesian coordinates.
    """
    return radius * _np.cos(theta), radius * _np.sin(theta)


@u.quantity_input
def cartesian_to_polar(x: u.m, y: u.m):
    """Convert 2D cartesian coordinates to polar coordinates.

    Parameters
    ----------
    x, y : ~astropy.units.Quantity
        Cartesian coordinates.

    Returns
    -------
    theta : ~astropy.units.Quantity
        Polar angle.
    radius : ~astropy.units.Quantity
        Radius.
    """

    radius = _np.sqrt(x**2 + y**2)
    theta = _np.arctan2(y, x)

    return theta, radius

1	"""This module contains functions for coordinate transformations"""
2
3	import functools as _functools	1✔
4
5	import astropy.units as u	1✔
6	import numpy as _np	1✔
7
8	##############################################################################
9	# 3D coordinates
10	##############################################################################
11
12
13	@u.quantity_input	1✔
14	def geodetic_to_cartesian(lat: u.deg, lon: u.deg, height: u.m, ell):	1✔
15	"""Convert geodetic to 3D cartesian coordinates.
16
17	Convert geodetic coordinates (`lat`, `lon`, `height`) given w.r.t.
18	ellipsoid `ell` to 3D cartesian coordinates (`x`, `y`, `z`).
19
20	Parameters
21	----------
22	lat : ~astropy.units.Quantity
23	Geodetic latitude.
24	lon : ~astropy.units.Quantity
25	Geodetic longitude.
26	height : ~astropy.units.Quantity
27	Geodetic height.
28	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
29	Reference ellipsoid to which geodetic coordinates are referenced to.
30
31	Returns
32	-------
33	x, y, z : ~astropy.units.Quantity
34	Cartesian coordinates.
35	"""
36
37	N = ell.prime_vertical_curvature_radius(lat)	1✔
38	x = (N + height) * _np.cos(lat) * _np.cos(lon)	1✔
39	y = (N + height) * _np.cos(lat) * _np.sin(lon)	1✔
40	z = (N + height - N * ell.e2) * _np.sin(lat)	1✔
41
42	return x, y, z	1✔
43
44
45	@_functools.partial(	1✔
46	_np.vectorize, otypes=(_np.float64, _np.float64, _np.float64), excluded=[3, 4]
47	)
48	def _cartesian_to_geodetic(x, y, z, ell, degrees=True):	1✔
49	"""Convert 3D cartesian to geodetic coordinates.
50
51	Parameters
52	----------
53	x, y, z : float or array_like of floats
54	Cartesian coordinates, in metres.
55	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
56	Reference ellipsoid to which geodetic coordinates are referenced to.
57	degrees : bool, optional
58	If True, the output geodetic latitude and longitude will be in degrees,
59	otherwise radians.
60
61	Returns
62	-------
63	lat, lon : float or array_like of floats
64	Geodetic latitude and longitude.
65	height : float or array_like of floats
66	Geodetic height, in metres.
67
68	Notes
69	-----
70	The algorithm of H. Vermeille is used for this transformation [1]_.
71
72	References
73	----------
74	.. [1] Vermeille, H., 2011. An analytical method to transform geocentric
75	into geodetic coordinates. Journal of Geodesy, 85(2), pp.105-117.
76	"""
77	e2 = ell.e2.value	1✔
78	a = ell.a.value	1✔
79
80	e4 = e2**2	1✔
81
82	# Step 1
83	p = (x2 + y2) / a**2	1✔
84	q = (1 - e2) * z2 / a2	1✔
85	r = (p + q - e4) / 6	1✔
86
87	# Step 2 - 3
88	e4pq = e4 * p * q	1✔
89	t = 8 * r**3 + e4pq	1✔
90
91	if (t > 0) or (t <= 0 and q != 0):	1✔
92	if t > 0:	1✔
93	li = _np.power(_np.sqrt(t) + _np.sqrt(e4pq), 1 / 3)	1✔
94	u = 3 / 2 * r2 / li2 + 0.5 * (li + r / li) ** 2	1✔
95	elif t <= 0 and q != 0:	1✔
96	u_aux = (	1✔
97	2 / 3 * _np.arctan2(_np.sqrt(e4pq), _np.sqrt(-t) + _np.sqrt(-8 * r**3))
98	)
99
100	u = -4 * r * _np.sin(u_aux) * _np.cos(_np.pi / 6 + u_aux)	1✔
101
102	v = _np.sqrt(u*2 + e4 q)	1✔
103	w = e2 * (u + v - q) / (2 * v)	1✔
104	k = (u + v) / (_np.sqrt(w**2 + u + v) + w)	1✔
105	D = (k * _np.sqrt(x2 + y2)) / (k + e2)	1✔
106
107	height = (k + e2 - 1) * _np.sqrt(D2 + z2) / k	1✔
108	lat = 2 * _np.arctan2(z, D + _np.sqrt(D2 + z2))	1✔
109	# Step 4
110	elif q == 0 and p <= e4:	×
111	e2p = _np.sqrt(e2 - p)	×
112	me2 = _np.sqrt(1 - e2)	×
113
NEW 114	height = -(a * me2 * e2p) / _np.sqrt(e2)	×
NEW 115	lat = 2 * _np.arctan2(_np.sqrt(e4 - p), _np.sqrt(e2) * e2p + me2 * _np.sqrt(p))	×
116
117	lon = _np.arctan2(y, x)	1✔
118
119	if degrees:	1✔
120	lat = _np.degrees(lat)	1✔
121	lon = _np.degrees(lon)	1✔
122
123	return lat, lon, height	1✔
124
125
126	@u.quantity_input	1✔
127	def cartesian_to_geodetic(x: u.m, y: u.m, z: u.m, ell):	1✔
128	"""Convert 3D cartesian to geodetic coordinates.
129
130	Parameters
131	----------
132	x, y, z : ~astropy.units.Quantity
133	Cartesian coordinates.
134	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
135	Reference ellipsoid to which geodetic coordinates are referenced to.
136
137	Returns
138	-------
139	lat : ~astropy.units.Quantity
140	Geodetic latitude.
141	lon : ~astropy.units.Quantity
142	Geodetic longitude.
143	height : float or array_like of floats
144	Geodetic height.
145
146	Notes
147	-----
148	The algorithm of H. Vermeille is used for this transformation [1]_.
149
150	References
151	----------
152	.. [1] Vermeille, H., 2011. An analytical method to transform geocentric
153	into geodetic coordinates. Journal of Geodesy, 85(2), pp.105-117.
154	"""
155
156	lat, lon, height = _cartesian_to_geodetic(	1✔
157	x.to("m").value, y.to("m").value, z.to("m").value, ell=ell, degrees=True
158	)
159
160	return lat * u.deg, lon * u.deg, height * u.m	1✔
161
162
163	@u.quantity_input	1✔
164	def cartesian_to_spherical(x: u.m, y: u.m, z: u.m):	1✔
165	"""Convert 3D cartesian to spherical coordinates.
166
167	Parameters
168	----------
169	x, y, z : ~astropy.units.Quantity
170	Cartesian coordinates.
171
172	Returns
173	-------
174	lat : ~astropy.units.Quantity
175	Spherical latitude.
176	lon : ~astropy.units.Quantity
177	Spherical longitude.
178	r : ~astropy.units.Quantity
179	Radius.
180	"""
181
182	radius = _np.sqrt(x2 + y2 + z**2)	1✔
183	lat = _np.arctan2(z, _np.sqrt(x2 + y2))	1✔
184	lon = _np.arctan2(y, x)	1✔
185
186	return lat, lon, radius	1✔
187
188
189	@u.quantity_input	1✔
190	def spherical_to_cartesian(lat: u.deg, lon: u.deg, radius: u.m):	1✔
191	"""Convert spherical to 3D cartesian coordinates.
192
193	Parameters
194	----------
195	lat : ~astropy.units.Quantity
196	Spherical latitude.
197	lon : ~astropy.units.Quantity
198	Spherical longitude.
199	r : ~astropy.units.Quantity
200	Radius.
201
202	Returns
203	-------
204	x, y, z : ~astropy.units.Quantity
205	Cartesian coordinates.
206	"""
207
208	x = radius * _np.cos(lat) * _np.cos(lon)	1✔
209	y = radius * _np.cos(lat) * _np.sin(lon)	1✔
210	z = radius * _np.sin(lat)	1✔
211
212	return x, y, z	1✔
213
214
215	@u.quantity_input	1✔
216	def cartesian_to_ellipsoidal(x: u.m, y: u.m, z: u.m, ell):	1✔
217	"""Convert 3D cartesian to ellipsoidal-harmonic coordinates.
218
219	Note that point (x, y, z) must be on or outside of the sphere with the
220	radius equals to the linear eccentricity of the reference ellipsoid `ell`,
221	i. e. (x2 + y2 + z2) >= E2.
222
223	Parameters
224	----------
225	x, y, z : ~astropy.units.Quantity
226	Cartesian coordinates.
227	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
228	Reference ellipsoid to which ellipsoidal coordinates are referenced to.
229
230	Returns
231	-------
232	rlat : ~astropy.units.Quantity
233	Reduced latitude.
234	lon : ~astropy.units.Quantity
235	Longitude.
236	u_ax : ~astropy.units.Quantity
237	Polar axis of the ellipsoid passing through the given point.
238	"""
239
240	le2 = ell.linear_eccentricity**2	1✔
241
242	k = x2 + y2 + z**2 - le2	1✔
243
244	if _np.any(k < 0):	1✔
245	raise ValueError(	1✔
246	"x2 + y2 + z**2 must be grater or equal to "
247	+ "the linear eccentricity of the reference ellipsoid."
248	)
249
250	u_ax = k * (0.5 + 0.5 * _np.sqrt(1 + (4 * le2 * z2) / k2))	1✔
251
252	u_ax = _np.sqrt(u_ax)	1✔
253	rlat = _np.arctan2(z * _np.sqrt(u_ax*2 + le2), u_ax _np.sqrt(x2 + y2))	1✔
254
255	lon = _np.arctan2(y, x)	1✔
256
257	return rlat, lon, u_ax	1✔
258
259
260	@u.quantity_input	1✔
261	def ellipsoidal_to_cartesian(rlat: u.deg, lon: u.deg, u_ax: u.m, ell):	1✔
262	"""Convert ellipsoidal-harmonic coordinates to 3D cartesian coordinates.
263
264	Parameters
265	----------
266	rlat : ~astropy.units.Quantity
267	Reduced latitude.
268	lon : ~astropy.units.Quantity
269	Longitude.
270	u_ax : ~astropy.units.Quantity
271	Polar axis of the ellipsoid passing through the given point.
272	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
273	Reference ellipsoid to which geodetic coordinates are referenced to.
274
275	Returns
276	-------
277	x, y, z : ~astropy.units.Quantity
278	Cartesian coordinates.
279	"""
280
281	k = _np.sqrt(u_ax2 + ell.linear_eccentricity2)	1✔
282
283	x = k * _np.cos(rlat) * _np.cos(lon)	1✔
284	y = k * _np.cos(rlat) * _np.sin(lon)	1✔
285	z = u_ax * _np.sin(rlat)	1✔
286
287	return x, y, z	1✔
288
289
290	@u.quantity_input	1✔
291	def geodetic_to_spherical(lat: u.deg, lon: u.deg, height: u.m, ell):	1✔
292	"""Convert from geodetic to spherical coordinates.
293
294	Parameters
295	----------
296	lat : ~astropy.units.Quantity
297	Geodetic latitude.
298	lon : ~astropy.units.Quantity
299	Geodetic longitude.
300	height : ~astropy.units.Quantity
301	Geodetic height.
302	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
303	Reference ellipsoid to which geodetic coordinates are referenced to.
304
305	Returns
306	-------
307	lat, lon : ~astropy.units.Quantity
308	Spherical latitude and longitude.
309	r : ~astropy.units.Quantity
310	Radius.
311	"""
312	return cartesian_to_spherical(*geodetic_to_cartesian(lat, lon, height, ell=ell))	1✔
313
314
315	@u.quantity_input	1✔
316	def spherical_to_geodetic(lat: u.deg, lon: u.deg, radius: u.m, ell):	1✔
317	"""Convert spherical to geodetic coordinates.
318
319	Parameters
320	----------
321	lat, lon : ~astropy.units.Quantity
322	Spherical latitude and longitude.
323	r : ~astropy.units.Quantity
324	Radius.
325	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
326	Reference ellipsoid to which geodetic coordinates are referenced to.
327
328	Returns
329	-------
330	lat : ~astropy.units.Quantity
331	Geodetic latitude.
332	lon : ~astropy.units.Quantity
333	Geodetic longitude.
334	height : ~astropy.units.Quantity
335	Geodetic height.
336	"""
337	return cartesian_to_geodetic(*spherical_to_cartesian(lat, lon, radius), ell=ell)	1✔
338
339
340	@u.quantity_input	1✔
341	def geodetic_to_ellipsoidal(lat: u.deg, lon: u.deg, height: u.m, ell):	1✔
342	"""Convert from geodetic to ellipsoidal-harmonic coordinates.
343
344	Parameters
345	----------
346	lat : ~astropy.units.Quantity
347	Geodetic latitude.
348	lon : ~astropy.units.Quantity
349	Geodetic longitude.
350	height : ~astropy.units.Quantity
351	Geodetic height.
352	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
353	Reference ellipsoid to which geodetic coordinates are referenced to.
354
355	Returns
356	-------
357	rlat : ~astropy.units.Quantity
358	Reduced latitude.
359	lon : ~astropy.units.Quantity
360	Longitude.
361	u_ax : ~astropy.units.Quantity
362	Polar axis of the ellipsoid passing through the given point.
363
364	"""
365	return cartesian_to_ellipsoidal(	1✔
366	*geodetic_to_cartesian(lat, lon, height, ell=ell), ell=ell
367	)
368
369
370	@u.quantity_input	1✔
371	def ellipsoidal_to_geodetic(rlat: u.deg, lon: u.deg, u_ax: u.m, ell):	1✔
372	"""Convert from ellipsoidal-harmonic to geodetic coordinates.
373
374	Parameters
375	----------
376	rlat : ~astropy.units.Quantity
377	Reduced latitude.
378	lon : ~astropy.units.Quantity
379	Longitude.
380	u_ax : ~astropy.units.Quantity
381	Polar axis of the ellipsoid passing through the given point.
382	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
383	Reference ellipsoid to which geodetic coordinates are referenced to.
384
385	Returns
386	-------
387	lat, lon : ~astropy.units.Quantity
388	Geodetic latitude and longitude.
389	height : ~astropy.units.Quantity
390	Geodetic height.
391	"""
392	return cartesian_to_geodetic(	1✔
393	*ellipsoidal_to_cartesian(rlat, lon, u_ax, ell=ell), ell=ell
394	)
395
396
397	@u.quantity_input	1✔
398	def _ecef_to_enu_rotation_matrix(lat: u.deg, lon: u.deg):	1✔
399	"""Return ECEF to ENU rotation matrix.
400
401	Parameters
402	----------
403	lat : ~astropy.units.Quantity
404	Geodetic or spherical latitude.
405	lon : ~astropy.units.Quantity
406	Geodetic or spherical longitude.
407
408	"""
409	clat = _np.cos(lat)	1✔
410	slat = _np.sin(lat)	1✔
411	clon = _np.cos(lon)	1✔
412	slon = _np.sin(lon)	1✔
413
414	rotation_matrix = _np.array(	1✔
415	[
416	[-slon, clon, 0],
417	[-slat * clon, -slat * slon, clat],
418	[clat * clon, clat * slon, slat],
419	]
420	)
421
422	return rotation_matrix	1✔
423
424
425	@u.quantity_input	1✔
426	def ecef_to_enu(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell=None):	1✔
427	"""Convert geocentric cartesian to local cartesian coordinates.
428
429	Convert Earth Centered Earth Fixed (ECEF) cartesian coordinates
430	(`x`,`y`,`z`) to the local east-north-up (ENU) local cartesian
431	coordinate system with the origin in (`lat0`, `lon0`, `height0`)
432	or in (`lat0`, `lon0`, `r0`).
433
434	Parameters
435	----------
436	x, y, z : ~astropy.units.Quantity
437	Geocentric cartesian coordinates.
438	origin : tuple of ~astropy.units.Quantity
439	Ggeocentric (spherical) or geodetic coordinates of the origin
440	(`lat0`, `lon0`, `r0`) or (`lat0`, `lon0`, `h0`).
441	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`, optional
442	Reference ellipsoid to which geodetic coordinates are referenced to.
443	Default is None, meaning spherical coordinates instead of geodetic.
444
445	Returns
446	-------
447	x, y, z : ~astropy.units.Quantity
448	Local east-north-up cartesian coordinates.
449	"""
450	rotation_matrix = _ecef_to_enu_rotation_matrix(origin[0], origin[1])	1✔
451	if ell is None:	1✔
452	x0, y0, z0 = spherical_to_cartesian(*origin)	1✔
453	else:
454	x0, y0, z0 = geodetic_to_cartesian(*origin, ell=ell)	1✔
455
456	out_shape = x.shape	1✔
457
458	xyz_shifted = _np.array(	1✔
459	[
460	_np.asarray((x - x0).to("m").value).flatten(),
461	_np.asarray((y - y0).to("m").value).flatten(),
462	_np.asarray((z - z0).to("m").value).flatten(),
463	]
464	)
465
466	out = _np.dot(rotation_matrix, xyz_shifted) * u.m	1✔
467
468	return (	1✔
469	out[0].reshape(out_shape),
470	out[1].reshape(out_shape),
471	out[2].reshape(out_shape),
472	)
473
474
475	@u.quantity_input	1✔
476	def enu_to_ecef(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell=None):	1✔
477	"""Convert local cartesian to geocentric cartesian coordinates.
478
479	Convert local east-north-up (ENU) local cartesian
480	coordinate system with the origin in (`lat0`, `lon0`, `height0`)
481	to the Earth Centered Earth Fixed (ECEF) cartesian coordinates.
482
483	Parameters
484	----------
485	x, y, z : ~astropy.units.Quantity
486	Local east-north-uo cartesian coordinates.
487	origin : tuple of ~astropy.units.Quantity
488	Ggeocentric (spherical) or geodetic coordinates of the origin
489	(`lat0`, `lon0`, `r0`) or (`lat0`, `lon0`, `h0`).
490	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
491	Reference ellipsoid to which geodetic coordinates are referenced to.
492	Default is None, meaning spherical coordinates instead of geodetic.
493
494	Returns
495	-------
496	x, y, z : ~astropy.units.Quantity
497	Geocentric cartesian coordinates, in metres.
498	"""
499	rotation_matrix = _ecef_to_enu_rotation_matrix(origin[0], origin[1]).T	1✔
500
501	out_shape = x.shape	1✔
502
503	x, y, z = (	1✔
504	_np.dot(
505	rotation_matrix,
506	_np.array(
507	[
508	_np.asarray(x.to("m").value).flatten(),
509	_np.asarray(y.to("m").value).flatten(),
510	_np.asarray(z.to("m").value).flatten(),
511	]
512	),
513	)
514	* u.m
515	)
516
517	if ell is None:	1✔
518	x0, y0, z0 = spherical_to_cartesian(*origin)	1✔
519	else:
520	x0, y0, z0 = geodetic_to_cartesian(*origin, ell=ell)	1✔
521
522	return (	1✔
523	(x + x0).reshape(out_shape),
524	(y + y0).reshape(out_shape),
525	(z + z0).reshape(out_shape),
526	)
527
528
529	@u.quantity_input	1✔
530	def geodetic_to_enu(	1✔
531	lat: u.deg, lon: u.deg, height: u.m, origin: tuple[u.deg, u.deg, u.m], ell
532	):
533	"""Convert geodetic coordinates to local cartesian coordinates.
534
535	Convert geodetic coordinates
536	(`lat`,`lon`,`height`) to the local east-north-up (ENU) local cartesian
537	coordinate system with the origin in (`lat0`, `lon0`, `height0`).
538
539	Parameters
540	----------
541	lat : ~astropy.units.Quantity
542	Geodetic latitude.
543	lon : ~astropy.units.Quantity
544	Geodetic longitude.
545	height : ~astropy.units.Quantity
546	Geodetic height.
547	origin : tuple of ~astropy.units.Quantity
548	Geodetic coordinates of the origin (`lat0`, `lon0`, `h0`).
549	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
550	Reference ellipsoid to which geodetic coordinates are referenced to.
551
552	Returns
553	-------
554	x, y, z : ~astropy.units.Quantity
555	Local east-north-up cartesian coordinates.
556	"""
557	return ecef_to_enu(	1✔
558	*geodetic_to_cartesian(lat, lon, height, ell=ell), origin=origin, ell=ell
559	)
560
561
562	@u.quantity_input	1✔
563	def enu_to_geodetic(x: u.m, y: u.m, z: u.m, origin: tuple[u.deg, u.deg, u.m], ell):	1✔
564	"""Convert local cartesian coordinates to geodetic coordinates.
565
566	Convert the local east-north-up (ENU) local cartesian
567	coordinate system with the origin in (`lat0`, `lon0`, `height0`)
568	to the geodetic coordinates.
569
570	Parameters
571	----------
572	x, y, z : ~astropy.units.Quantity
573	Local east-north-uo cartesian coordinates.
574	origin : tuple of ~astropy.units.Quantity
575	Geodetic coordinates of the origin (`lat0`, `lon0`, `h0`).
576	ell : instance of the `pygeoid.coordinates.ellipsoid.Ellipsoid`
577	Reference ellipsoid to which geodetic coordinates are referenced to.
578
579	Returns
580	-------
581	lat : ~astropy.units.Quantity
582	Geodetic latitude.
583	lon : ~astropy.units.Quantity
584	Geodetic longitude.
585	height : ~astropy.units.Quantity
586	Geodetic height.
587	"""
588	return cartesian_to_geodetic(*enu_to_ecef(x, y, z, origin=origin, ell=ell), ell=ell)	1✔
589
590
591	##############################################################################
592	# 2D coordinates
593	##############################################################################
594
595
596	@u.quantity_input	1✔
597	def polar_to_cartesian(theta: u.deg, radius: u.m):	1✔
598	"""Convert polar coordinates to 2D cartesian.
599
600	Parameters
601	----------
602	theta : ~astropy.units.Quantity
603	Polar angle.
604	radius : ~astropy.units.Quantity
605	Radius.
606
607	Returns
608	-------
609	x, y : ~astropy.units.Quantity
610	Local east-north-uo cartesian coordinates.
611	"""
612	return radius * _np.cos(theta), radius * _np.sin(theta)	1✔
613
614
615	@u.quantity_input	1✔
616	def cartesian_to_polar(x: u.m, y: u.m):	1✔
617	"""Convert 2D cartesian coordinates to polar coordinates.
618
619	Parameters
620	----------
621	x, y : ~astropy.units.Quantity
622	Cartesian coordinates.
623
624	Returns
625	-------
626	theta : ~astropy.units.Quantity
627	Polar angle.
628	radius : ~astropy.units.Quantity
629	Radius.
630	"""
631
632	radius = _np.sqrt(x2 + y2)	1✔
633	theta = _np.arctan2(y, x)	1✔
634
635	return theta, radius	1✔

ioshchepkov / pygeoid / 21944858113

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous