Coordinate Systems

Coordinate system conversions and transformations.

This module provides functions for converting between different coordinate systems commonly used in tracking applications:

Cartesian coordinates (x, y, z)
Spherical coordinates (range, azimuth, elevation)
Polar and cylindrical coordinates
Geodetic coordinates (latitude, longitude, altitude)
Various local tangent plane frames (ENU, NED)
Direction cosine representations (r-u-v)
Rotation representations (matrices, quaternions, Euler angles)
Jacobian matrices for error propagation

Conversions

Coordinate conversions.

This module provides: - Spherical/polar coordinate conversions - Geodetic (lat/lon/alt) to ECEF conversions - Local tangent plane frames (ENU, NED) - Direction cosine representations (r-u-v)

Spherical Coordinates

Spherical and polar coordinate conversions.

This module provides functions for converting between Cartesian and spherical/polar coordinate systems, following tracking conventions.

pytcl.coordinate_systems.conversions.spherical.cart2sphere(cart_points, system_type='standard')[source]

Convert Cartesian coordinates to spherical coordinates.

Parameters:

cart_points (array_like) – Cartesian coordinates. Can be: - Shape (3,) for a single point [x, y, z] - Shape (3, n) for n points (each column is a point) - Shape (n, 3) will be transposed automatically
system_type ({'standard', 'az-el', 'range-az-el'}, optional) – Spherical coordinate convention: - ‘standard’: Physics convention (r, θ polar from +z, φ azimuth from +x) - ‘az-el’: Tracking convention (r, azimuth from +x, elevation from xy-plane) - ‘range-az-el’: Same as ‘az-el’ (alias) Default is ‘standard’.

Returns:

r (ndarray) – Range (radial distance from origin).
az (ndarray) – Azimuth angle in radians. - ‘standard’: Angle in xy-plane from +x axis [0, 2π) - ‘az-el’: Angle in xy-plane from +x axis [-π, π]
el (ndarray) – Elevation/polar angle in radians. - ‘standard’: Polar angle from +z axis [0, π] - ‘az-el’: Elevation from xy-plane [-π/2, π/2]

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> x, y, z = 1.0, 1.0, 1.0
>>> r, az, el = cart2sphere([x, y, z], system_type='az-el')
>>> r
1.7320508075688772
>>> np.degrees(az)
45.0
>>> round(float(np.degrees(el)), 6)
35.26439

See also

sphere2cart: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.sphere2cart(r, az, el, system_type='standard')[source]

Convert spherical coordinates to Cartesian coordinates.

Parameters:

r (array_like) – Range (radial distance).
az (array_like) – Azimuth angle in radians.
el (array_like) – Elevation/polar angle in radians.
system_type ({'standard', 'az-el', 'range-az-el'}, optional) – Spherical coordinate convention (see cart2sphere).

Returns:

cart_points – Cartesian coordinates of shape (3,) or (3, n).

Return type:

ndarray

Examples

>>> r, az, el = 1.732, np.radians(45), np.radians(35.26)
>>> cart2sphere(sphere2cart(r, az, el, 'az-el'), 'az-el')
(1.732..., 0.785..., 0.615...)

See also

cart2sphere: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.cart2pol(cart_points)[source]

Convert 2D Cartesian coordinates to polar coordinates.

Parameters:

cart_points (array_like) – Cartesian coordinates. Can be: - Shape (2,) for a single point [x, y] - Shape (2, n) for n points - Shape (n, 2) will be transposed

Returns:

r (ndarray) – Radial distance from origin.
theta (ndarray) – Angle in radians from +x axis, in range [-π, π].

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, theta = cart2pol([1, 1])
>>> r
1.4142135623730951
>>> np.degrees(theta)
45.0

See also

pol2cart: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.pol2cart(r, theta)[source]

Convert polar coordinates to 2D Cartesian coordinates.

Parameters:

r (array_like) – Radial distance.
theta (array_like) – Angle in radians from +x axis.

Returns:

cart_points – Cartesian coordinates of shape (2,) or (2, n).

Return type:

ndarray

Examples

>>> x, y = pol2cart(1.414, np.radians(45))
>>> x, y
(0.999..., 0.999...)

See also

cart2pol: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.cart2cyl(cart_points)[source]

Convert 3D Cartesian coordinates to cylindrical coordinates.

Parameters:

cart_points (array_like) – Cartesian coordinates [x, y, z].

Returns:

rho (ndarray) – Radial distance in xy-plane.
phi (ndarray) – Azimuth angle in radians from +x axis.
z (ndarray) – Height (same as Cartesian z).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> rho, phi, z = cart2cyl([1, 1, 5])
>>> rho
1.4142135623730951
>>> np.degrees(phi)
45.0
>>> z
5.0

See also

cyl2cart: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.cyl2cart(rho, phi, z)[source]

Convert cylindrical coordinates to 3D Cartesian coordinates.

Parameters:

rho (array_like) – Radial distance in xy-plane.
phi (array_like) – Azimuth angle in radians from +x axis.
z (array_like) – Height.

Returns:

cart_points – Cartesian coordinates of shape (3,) or (3, n).

Return type:

ndarray

Examples

>>> cart = cyl2cart(np.sqrt(2), np.radians(45), 5.0)
>>> cart
array([1., 1., 5.])

See also

cart2cyl: Inverse conversion.

pytcl.coordinate_systems.conversions.spherical.ruv2cart(r, u, v)[source]

Convert r-u-v (range, direction cosines) to Cartesian coordinates.

The r-u-v system uses direction cosines where: - u = cos(az) * cos(el) = x/r - v = sin(az) * cos(el) = y/r - w = sin(el) = z/r (derived from u, v)

Parameters:

r (array_like) – Range.
u (array_like) – Direction cosine along x-axis.
v (array_like) – Direction cosine along y-axis.

Returns:

cart_points – Cartesian coordinates.

Return type:

ndarray

Examples

>>> # Target at 45 deg azimuth, 30 deg elevation, range 100
>>> az, el = np.radians(45), np.radians(30)
>>> u = np.cos(az) * np.cos(el)
>>> v = np.sin(az) * np.cos(el)
>>> cart = ruv2cart(100, u, v)
>>> cart
array([61.23724357, 61.23724357, 50.        ])

Notes

This representation is common in radar tracking systems.

The upper hemisphere is assumed. u and v fix only the x and y direction cosines; the third is recovered as w = sqrt(1 - u**2 - v**2), which is never negative. A target at z < 0 therefore comes back with z mirrored, so ruv2cart(*cart2ruv(p)) round-trips only for p[2] >= 0. The r-u-v triple does not carry the hemisphere; keep it alongside if the geometry spans both, or use cart2sphere(), whose elevation is signed.

pytcl.coordinate_systems.conversions.spherical.cart2ruv(cart_points)[source]

Convert Cartesian coordinates to r-u-v (range, direction cosines).

Parameters:

cart_points (array_like) – Cartesian coordinates [x, y, z].

Returns:

r (ndarray) – Range.
u (ndarray) – Direction cosine along x-axis (x/r).
v (ndarray) – Direction cosine along y-axis (y/r).

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> r, u, v = cart2ruv([100, 0, 0])
>>> r, u, v
(100.0, 1.0, 0.0)
>>> r, u, v = cart2ruv([50, 50, 50])
>>> r
86.602...

See also

ruv2cart: Inverse conversion.

Notes

The returned (r, u, v) does not encode which hemisphere the point is in: u and v are the x and y direction cosines only. ruv2cart() assumes z >= 0 when inverting, so a point below the x-y plane does not survive the round trip.

u-v Direction Cosines

The angle-only u-v(-w) measurement system of planar phased arrays, and the full bistatic r-u-v conversions.

Direction-cosine u-v measurement coordinates.

The u-v(-w) system is the natural measurement space of a planar phased array: u and v are the first two components of a unit direction vector in the sensor’s coordinate system. This module ports the MATLAB TCL angle-only u-v conversions and the full bistatic r-u-v conversions with sensor offsets and pointing rotations (the simplified aligned-monostatic ruv2cart/cart2ruv live in pytcl.coordinate_systems.conversions.spherical).

All conventions follow D. F. Crouse, “Basic tracking using nonlinear 3D monostatic and bistatic measurements,” IEEE Aerospace and Electronic Systems Magazine, vol. 29, no. 8, Part II, pp. 4-53, Aug. 2014.

pytcl.coordinate_systems.conversions.uv.camera_coords2uv(z_cam, a, m=None, include_w=True)[source]

Convert camera pixel coordinates to u-v(-w) direction cosines.

Parameters:

z_cam (array_like) – (2, N) camera coordinates [x, y] in the image plane.
a (array_like) – (3, 3) camera intrinsics-style matrix; its third row must be [0, 0, a33].
m (array_like, optional) – (3, 3) rotation from the global frame to the camera’s frame; the returned directions are rotated back into the global frame. Default identity (no rotation applied).
include_w (bool, optional) – Include the third component. Default True.

Returns:

dir_vecs – (3, N) u-v-w unit vectors, or (2, N) u-v pairs when include_w is False.

Return type:

ndarray

Examples

>>> import numpy as np
>>> a = np.diag([500.0, 500.0, 1.0])
>>> d = camera_coords2uv(np.array([0.0, 0.0]), a)
>>> np.round(d.ravel(), 12)
array([0., 0., 1.])

Notes

Port of cameraCoords2UVCoords.m.

pytcl.coordinate_systems.conversions.uv.cart2ruv_bistatic(z_c, use_half_range=False, z_tx=None, z_rx=None, m=None, include_w=False)[source]

Convert global Cartesian points to bistatic r-u-v(-w) measurements.

Parameters:

z_c (array_like) – (3, N) global Cartesian positions.
use_half_range (bool, optional) – Halve the returned bistatic range (monostatic convention). Default False.
z_tx (array_like, optional) – Transmitter/receiver positions and receiver rotations as in ruv2cart_bistatic().
z_rx (array_like, optional) – Transmitter/receiver positions and receiver rotations as in ruv2cart_bistatic().
m (array_like, optional) – Transmitter/receiver positions and receiver rotations as in ruv2cart_bistatic().
include_w (bool, optional) – Also return the third direction cosine w. Default False.

Returns:

z – (3, N) r-u-v or (4, N) r-u-v-w measurements.

Return type:

ndarray

Examples

>>> import numpy as np
>>> z_c = np.array([0.0, 0.0, 5.0])
>>> np.round(cart2ruv_bistatic(z_c).ravel(), 12)
array([10.,  0.,  0.])

Notes

Port of Cart2Ruv.m (Crouse 2014, see the module docstring).

pytcl.coordinate_systems.conversions.uv.ruv2cart_bistatic(z, use_half_range=False, z_tx=None, z_rx=None, m=None)[source]

Convert bistatic r-u-v(-w) measurements to global Cartesian points.

The full conversion with transmitter/receiver offsets and receiver pointing rotations; the aligned-monostatic simplification is pytcl.coordinate_systems.conversions.spherical.ruv2cart().

Parameters:

z (array_like) – (3, N) r-u-v or (4, N) r-u-v-w measurements. The range is the bistatic range (transmitter to target to receiver).
use_half_range (bool, optional) – True if the ranges are one-way (monostatic convention). Default False.
z_tx (array_like, optional) – (3, N) transmitter positions, or a single (3,) position shared by all measurements. Default: the origin.
z_rx (array_like, optional) – Receiver positions, like z_tx. Default: the origin.
m (array_like, optional) – (3, 3, N) rotations from the global frame to each receiver’s local frame (the local z-axis is the pointing direction), or a single (3, 3) shared rotation. Default identity.

Returns:

z_c – (3, N) global Cartesian positions.

Return type:

ndarray

Examples

>>> import numpy as np
>>> z = np.array([10.0, 0.0, 0.0])  # r-u-v along the local z-axis
>>> np.round(ruv2cart_bistatic(z).ravel(), 12)
array([0., 0., 5.])

Notes

Port of ruv2Cart.m (Crouse 2014, see the module docstring).

pytcl.coordinate_systems.conversions.uv.ruv2ruv(z, use_half_range=False, z_tx1=None, z_rx1=None, m1=None, z_tx2=None, z_rx2=None, m2=None, include_w=None)[source]

Convert bistatic r-u-v(-w) measurements between two bistatic pairs.

Converts measurements taken by one transmitter/receiver pair into the coordinate system of another (possibly rotated, displaced) pair, via Cartesian coordinates.

Parameters:

z (array_like) – (3, N) r-u-v or (4, N) r-u-v-w measurements.
use_half_range (bool or tuple of (bool, bool), optional) – One-way-range convention for the input and output systems; a scalar applies to both. Default False.
z_tx1 (array_like, optional) – Transmitter/receiver positions and receiver rotation of the system the measurements come from.
z_rx1 (array_like, optional) – Transmitter/receiver positions and receiver rotation of the system the measurements come from.
m1 (array_like, optional) – Transmitter/receiver positions and receiver rotation of the system the measurements come from.
z_tx2 (array_like, optional) – The same for the system converted into.
z_rx2 (array_like, optional) – The same for the system converted into.
m2 (array_like, optional) – The same for the system converted into.
include_w (bool, optional) – Include w in the output. Default: True when the input is (4, N), else False.

Returns:

z_new – The measurements in the second system.

Return type:

ndarray

Examples

>>> import numpy as np
>>> z = np.array([100.0, 0.3, -0.2])
>>> rx2 = np.array([10.0, -5.0, 2.0])
>>> out = ruv2ruv(z, False, None, None, None, None, rx2)
>>> back = ruv2ruv(out, False, None, rx2)
>>> np.allclose(back.ravel(), z)
True

Notes

Port of ruv2Ruv.m: composes ruv2cart_bistatic() and cart2ruv_bistatic().

pytcl.coordinate_systems.conversions.uv.spher_ang2uv(az_el, system_type=0, include_w=False, m_s=None, m_uv=None)[source]

Convert spherical azimuth and elevation to u-v direction cosines.

Parameters:

az_el (array_like) – (2, N) azimuth and elevation in radians.
system_type (int, optional) – Spherical axis convention; see uv2spher_ang().
include_w (bool, optional) – Also return the third unit-vector component w. Default False.
m_s (array_like, optional) – Rotations as in uv2spher_ang().
m_uv (array_like, optional) – Rotations as in uv2spher_ang().

Returns:

uv – (2, N) u-v pairs, or (3, N) u-v-w unit vectors when include_w is True.

Return type:

ndarray

Examples

>>> import numpy as np
>>> uv = spher_ang2uv(np.array([0.4, 0.7]), include_w=True)
>>> az_el = uv2spher_ang(uv)
>>> np.allclose(az_el.ravel(), [0.4, 0.7])
True

Notes

Port of spherAng2Uv.m (Crouse 2014, see the module docstring).

pytcl.coordinate_systems.conversions.uv.state_ruv2cart(x)[source]

Convert an r-u-v state with derivatives to a Cartesian state.

Handles 6-element states [r, u, v, rdot, udot, vdot] and 9-element states with second derivatives appended, converting position, velocity and (when present) acceleration into 3D Cartesian components. The measurement is monostatic with w > 0.

Parameters:: x (array_like) – (6, N) or (9, N) r-u-v states.
Returns:: cart_states – (6, N) or (9, N) Cartesian states [x, y, z, xdot, ydot, zdot (, xddot, yddot, zddot)].
Return type:: ndarray

Examples

>>> import numpy as np
>>> x = np.array([100.0, 0.0, 0.0, 5.0, 0.0, 0.0])
>>> np.round(state_ruv2cart(x).ravel(), 12)
array([  0.,   0., 100.,   0.,   0.,   5.])

Notes

Port of stateRuv2Cart.m. The local basis vectors u1 (radial), u2 and u3 and the coefficients follow the original expressions verbatim.

pytcl.coordinate_systems.conversions.uv.uv2spher_ang(uv, system_type=0, m_s=None, m_uv=None)[source]

Convert u-v direction cosines to spherical azimuth and elevation.

Parameters:

uv (array_like) – (2, N) u-v pairs, or (3, N) u-v-w unit vectors. With only u-v given, w is taken positive (in front of the sensor).
system_type (int, optional) –
Spherical axis convention:
- 0 (default): azimuth counterclockwise from x in the x-y plane, elevation up from the x-y plane.
- 1: azimuth counterclockwise from z in the z-x plane, elevation up from the z-x plane (z-axis boresight).
- 2: like 0, but the second angle is measured down from the z-axis (pi/2 - elevation).
- 3: like 0, but azimuth is measured clockwise from the y-axis (East-of-North bearings in an ENU frame).
m_s (array_like, optional) – (3, 3) rotation from the global frame to the frame the spherical angles are expressed in. Default identity.
m_uv (array_like, optional) – (3, 3) rotation from the global frame to the frame the u-v coordinates are expressed in. Default identity.

Returns:

az_el – (2, N) azimuth and elevation in radians.

Return type:

ndarray

Examples

>>> import numpy as np
>>> az_el = uv2spher_ang(np.array([0.0, 0.0]))  # the +z boresight
>>> np.round(az_el.ravel(), 6).tolist()
[0.0, 1.570796]

Notes

Port of uv2SpherAng.m (Crouse 2014, see the module docstring).

Geodetic Coordinates

Geodetic coordinate conversions.

This module provides functions for converting between geodetic (latitude, longitude, altitude) and Earth-centered coordinate systems (ECEF), as well as local tangent plane coordinates (ENU, NED).

pytcl.coordinate_systems.conversions.geodetic.geodetic2ecef(lat, lon, alt, a=6378137.0, f=0.0033528106647474805)[source]

Convert geodetic coordinates to ECEF (Earth-Centered Earth-Fixed).

Parameters:

lat (array_like) – Geodetic latitude in radians.
lon (array_like) – Geodetic longitude in radians.
alt (array_like) – Altitude above the reference ellipsoid in meters.
a (float, optional) – Semi-major axis of the reference ellipsoid in meters. Default is WGS84 value.
f (float, optional) – Flattening of the reference ellipsoid. Default is WGS84 value.

Returns:

ecef – ECEF coordinates [x, y, z] in meters. Shape is (3,) for single point or (3, n) for multiple points.

Return type:

ndarray

Examples

>>> lat, lon, alt = np.radians(45), np.radians(-75), 100.0
>>> ecef = geodetic2ecef(lat, lon, alt)
>>> np.round(ecef / 1e6, 3)  # In millions of meters
array([ 1.169, -4.364,  4.487])

See also

ecef2geodetic: Inverse conversion.

pytcl.coordinate_systems.conversions.geodetic.ecef2geodetic(ecef, a=6378137.0, f=0.0033528106647474805, method='iterative')[source]

Convert ECEF coordinates to geodetic coordinates.

Parameters:

ecef (array_like) – ECEF coordinates [x, y, z] in meters.
a (float, optional) – Semi-major axis of the reference ellipsoid.
f (float, optional) – Flattening of the reference ellipsoid.
method (str, optional) – Algorithm to use: - ‘iterative’: Bowring’s iterative method (default) - ‘direct’: Closed-form solution (Vermeille’s method)

Returns:

lat (ndarray) – Geodetic latitude in radians.
lon (ndarray) – Geodetic longitude in radians.
alt (ndarray) – Altitude above the ellipsoid in meters.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> ecef = geodetic2ecef(np.radians(45), np.radians(-75), 100.0)
>>> lat, lon, alt = ecef2geodetic(ecef)
>>> round(float(np.degrees(lat)), 6), round(float(np.degrees(lon)), 6)
(45.0, -75.0)
>>> round(float(alt), 3)
100.0

See also

geodetic2ecef: Inverse conversion.

pytcl.coordinate_systems.conversions.geodetic.geodetic2enu(lat, lon, alt, lat_ref, lon_ref, alt_ref, a=6378137.0, f=0.0033528106647474805)[source]

Convert geodetic coordinates to local ENU (East-North-Up) coordinates.

Parameters:

lat (array_like) – Geodetic latitude in radians.
lon (array_like) – Geodetic longitude in radians.
alt (array_like) – Altitude in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
alt_ref (float) – Reference point altitude in meters.
a (float, optional) – Semi-major axis of the reference ellipsoid.
f (float, optional) – Flattening of the reference ellipsoid.

Returns:

enu – Local ENU coordinates [east, north, up] in meters.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from pytcl.coordinate_systems import geodetic2enu
>>> # Reference point: San Francisco Airport (-122.375°, 37.615°)
>>> lat_ref = np.radians(37.615)
>>> lon_ref = np.radians(-122.375)
>>> alt_ref = 0
>>> # Target point: 2 km north, 1 km east of reference
>>> # Approximate location using small offsets
>>> lat_target = lat_ref + np.radians(0.01)
>>> lon_target = lon_ref + np.radians(0.01)
>>> alt_target = 100  # 100 m elevation
>>> enu = geodetic2enu(lat_target, lon_target, alt_target,
...                      lat_ref, lon_ref, alt_ref)
>>> # ENU coordinates should show positive north and east offsets
>>> enu[0] > 0 and enu[1] > 0  # East > 0, North > 0
True
>>> bool(abs(enu[2] - 100) < 1)  # Up approximately the altitude difference
True

See also

enu2geodetic: Inverse conversion.
ecef2enu: ECEF to ENU conversion.

pytcl.coordinate_systems.conversions.geodetic.ecef2enu(ecef, lat_ref, lon_ref, ecef_ref=None)[source]

Convert ECEF coordinates to local ENU coordinates.

Parameters:

ecef (array_like) – ECEF coordinates [x, y, z] in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference point ECEF coordinates. If None, computed from lat_ref, lon_ref.

Returns:

enu – Local ENU coordinates [east, north, up] in meters.

Return type:

ndarray

Examples

Convert an aircraft position from ECEF to local ENU frame:

>>> import numpy as np
>>> from pytcl.coordinate_systems.conversions import ecef2enu, geodetic2ecef
>>> # Reference point (airport at 38.9°N, 77.0°W)
>>> lat_ref = np.radians(38.9)
>>> lon_ref = np.radians(-77.0)
>>> # Aircraft ECEF position
>>> ecef_aircraft = np.array([1130000.0, -4830000.0, 3990000.0])
>>> enu = ecef2enu(ecef_aircraft, lat_ref, lon_ref)
>>> enu.shape
(3,)

See also

enu2ecef: Inverse conversion.
ecef2ned: Convert to NED (North-East-Down) frame.

pytcl.coordinate_systems.conversions.geodetic.enu2ecef(enu, lat_ref, lon_ref, ecef_ref=None)[source]

Convert local ENU coordinates to ECEF coordinates.

Parameters:

enu (array_like) – Local ENU coordinates [east, north, up] in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

ecef – ECEF coordinates [x, y, z] in meters.

Return type:

ndarray

Examples

Convert local ENU offset to ECEF coordinates:

>>> import numpy as np
>>> from pytcl.coordinate_systems.conversions import enu2ecef
>>> # Reference point (airport at 38.9°N, 77.0°W)
>>> lat_ref = np.radians(38.9)
>>> lon_ref = np.radians(-77.0)
>>> # Aircraft 1km east, 2km north, 500m up
>>> enu = np.array([1000.0, 2000.0, 500.0])
>>> ecef = enu2ecef(enu, lat_ref, lon_ref)
>>> ecef.shape
(3,)

See also

ecef2enu: Inverse conversion.

pytcl.coordinate_systems.conversions.geodetic.ecef2ned(ecef, lat_ref, lon_ref, ecef_ref=None)[source]

Convert ECEF coordinates to local NED (North-East-Down) coordinates.

Parameters:

ecef (array_like) – ECEF coordinates [x, y, z] in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

ned – Local NED coordinates [north, east, down] in meters.

Return type:

ndarray

Examples

Convert aircraft ECEF position to local NED frame:

>>> import numpy as np
>>> from pytcl.coordinate_systems.conversions import ecef2ned, geodetic2ecef
>>> # Reference point (airport)
>>> lat_ref = np.radians(38.9)
>>> lon_ref = np.radians(-77.0)
>>> # Aircraft ECEF position
>>> ecef_aircraft = np.array([1130000.0, -4830000.0, 3990000.0])
>>> ned = ecef2ned(ecef_aircraft, lat_ref, lon_ref)
>>> ned.shape  # [north, east, down]
(3,)

See also

ned2ecef: Inverse conversion.
ecef2enu: Similar, but ENU frame.

pytcl.coordinate_systems.conversions.geodetic.ned2ecef(ned, lat_ref, lon_ref, ecef_ref=None)[source]

Convert local NED coordinates to ECEF coordinates.

Parameters:

ned (array_like) – Local NED coordinates [north, east, down] in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

ecef – ECEF coordinates [x, y, z] in meters.

Return type:

ndarray

Examples

>>> import numpy as np
>>> from pytcl.coordinate_systems import ned2ecef, geodetic2ecef
>>> # Reference point: Kennedy Space Center (28.5°N, 80.65°W)
>>> lat_ref = np.radians(28.5)
>>> lon_ref = np.radians(-80.65)
>>> # NED offset from reference: 5 km north, 3 km east, 1 km down
>>> ned = np.array([5000.0, 3000.0, 1000.0])
>>> # Convert to ECEF coordinates
>>> ecef = ned2ecef(ned, lat_ref, lon_ref)
>>> ecef.shape
(3,)
>>> # Verify roundtrip conversion
>>> from pytcl.coordinate_systems import ecef2ned
>>> ned_back = ecef2ned(ecef, lat_ref, lon_ref)
>>> np.allclose(ned, ned_back, atol=0.1)  # Allow small numerical error
True

See also

ecef2ned: Inverse conversion.

pytcl.coordinate_systems.conversions.geodetic.enu2ned(enu)[source]

Convert ENU coordinates to NED coordinates.

Parameters:: enu (array_like) – ENU coordinates [east, north, up].
Returns:: ned – NED coordinates [north, east, down].
Return type:: ndarray

pytcl.coordinate_systems.conversions.geodetic.ned2enu(ned)[source]

Convert NED coordinates to ENU coordinates.

Parameters:: ned (array_like) – NED coordinates [north, east, down].
Returns:: enu – ENU coordinates [east, north, up].
Return type:: ndarray

pytcl.coordinate_systems.conversions.geodetic.geodetic2sez(lat, lon, alt, lat_ref, lon_ref, alt_ref, a=6378137.0, f=0.0033528106647474805)[source]

Convert geodetic coordinates to local SEZ (South-East-Zenith) coordinates.

SEZ is a horizon-relative coordinate frame where: - S (South) points in the southward direction - E (East) points in the eastward direction - Z (Zenith) points upward (away from Earth center)

Parameters:

lat (array_like) – Geodetic latitude in radians.
lon (array_like) – Geodetic longitude in radians.
alt (array_like) – Altitude in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
alt_ref (float) – Reference point altitude in meters.
a (float, optional) – Semi-major axis of the reference ellipsoid.
f (float, optional) – Flattening of the reference ellipsoid.

Returns:

sez – Local SEZ coordinates [south, east, zenith] in meters.

Return type:

ndarray

See also

sez2geodetic: Inverse conversion.
ecef2sez: ECEF to SEZ conversion.

Notes

Standard SEZ convention (Vallado): S points south, E east, Z up. Relative to NED: SEZ = [S, E, Z] = [-NED[0], NED[1], -NED[2]]

Examples

>>> import numpy as np
>>> # A target 0.1 deg NORTH of the reference has a NEGATIVE S component
>>> sez = geodetic2sez(
...     np.radians(45.1), np.radians(-75.0), 0.0,
...     np.radians(45.0), np.radians(-75.0), 0.0)
>>> bool(sez[0] < 0)
True
>>> round(float(sez[1]), 6)  # Due north: no east component
0.0

pytcl.coordinate_systems.conversions.geodetic.ecef2sez(ecef, lat_ref, lon_ref, ecef_ref=None)[source]

Convert ECEF coordinates to local SEZ coordinates.

Parameters:

ecef (array_like) – ECEF coordinates [X, Y, Z] in meters, shape (3,) or (3, N).
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference ECEF position. If None, the reference point is at (lat_ref, lon_ref) with zero altitude.

Returns:

sez – SEZ coordinates [south, east, zenith] in meters.

Return type:

ndarray

See also

sez2ecef: Inverse conversion.

pytcl.coordinate_systems.conversions.geodetic.sez2ecef(sez, lat_ref, lon_ref, ecef_ref=None)[source]

Convert local SEZ coordinates to ECEF coordinates.

Parameters:

sez (array_like) – SEZ coordinates [south, east, zenith] in meters, shape (3,) or (3, N).
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
ecef_ref (array_like, optional) – Reference ECEF position. If None, the reference point is at (lat_ref, lon_ref) with zero altitude.

Returns:

ecef – ECEF coordinates [X, Y, Z] in meters.

Return type:

ndarray

See also

ecef2sez: Forward conversion.

pytcl.coordinate_systems.conversions.geodetic.sez2geodetic(sez, lat_ref, lon_ref, alt_ref, a=6378137.0, f=0.0033528106647474805)[source]

Convert local SEZ coordinates to geodetic coordinates.

Parameters:

sez (array_like) – SEZ coordinates [south, east, zenith] in meters.
lat_ref (float) – Reference point latitude in radians.
lon_ref (float) – Reference point longitude in radians.
alt_ref (float) – Reference point altitude in meters.
a (float, optional) – Semi-major axis.
f (float, optional) – Flattening.

Returns:

lat (ndarray) – Geodetic latitude in radians.
lon (ndarray) – Geodetic longitude in radians.
alt (ndarray) – Altitude in meters.

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]]]

Examples

>>> import numpy as np
>>> from pytcl.coordinate_systems.conversions.geodetic import (
...     sez2geodetic, geodetic2sez)
>>> # Observer location: Arecibo Observatory (18.3°N, 66.75°W)
>>> lat_ref = np.radians(18.3)
>>> lon_ref = np.radians(-66.75)
>>> alt_ref = 500  # m
>>> # Observed satellite: 30 km south, 20 km east, 35 km zenith
>>> sez = np.array([-30000.0, 20000.0, 35000.0])
>>> # Convert to geodetic coordinates
>>> lat, lon, alt = sez2geodetic(sez, lat_ref, lon_ref, alt_ref)
>>> # Verify roundtrip conversion
>>> sez_back = geodetic2sez(lat, lon, alt, lat_ref, lon_ref, alt_ref)
>>> bool(np.allclose(sez, sez_back, atol=1.0))  # Allow ~1m numerical error
True

See also

geodetic2sez: Forward conversion.

pytcl.coordinate_systems.conversions.geodetic.geocentric_radius(lat, a=6378137.0, f=0.0033528106647474805)[source]

Compute the geocentric radius at a given geodetic latitude.

Parameters:

lat (array_like) – Geodetic latitude in radians.
a (float, optional) – Semi-major axis.
f (float, optional) – Flattening.

Returns:

r – Geocentric radius in meters.

Return type:

ndarray

pytcl.coordinate_systems.conversions.geodetic.prime_vertical_radius(lat, a=6378137.0, f=0.0033528106647474805)[source]

Compute the prime vertical radius of curvature.

Parameters:

lat (array_like) – Geodetic latitude in radians.
a (float, optional) – Semi-major axis.
f (float, optional) – Flattening.

Returns:

N – Prime vertical radius of curvature in meters.

Return type:

ndarray

pytcl.coordinate_systems.conversions.geodetic.meridional_radius(lat, a=6378137.0, f=0.0033528106647474805)[source]

Compute the meridional radius of curvature.

Parameters:

lat (array_like) – Geodetic latitude in radians.
a (float, optional) – Semi-major axis.
f (float, optional) – Flattening.

Returns:

M – Meridional radius of curvature in meters.

Return type:

ndarray

Rotations

Rotation representations and conversions.

This module provides: - Basic rotation matrices (rotx, roty, rotz) - Euler angle conversions - Quaternion operations - Axis-angle and Rodrigues representations - Rotation interpolation (SLERP)

pytcl.coordinate_systems.rotations.rot_axis_to_vec(u, axis=None)[source]

Rotation matrix taking a coordinate axis into the direction of u.

Returns R such that u/||u|| = R @ e_axis where e_axis is the selected coordinate axis unit vector. Works in any number of dimensions >= 1.

Parameters:

u (array_like) – Direction vector of shape (n,), or (n, N) for a batch of N vectors. Non-unit vectors are normalized.
axis (int or str, optional) – The axis to rotate: an index from 0 to n-1, or ‘x’, ‘y’, ‘z’ for 0, 1, 2. Default: ‘z’ (index 2) in 3D, index 0 otherwise.

Returns:

R – (n, n) rotation matrix, or (n, n, N) for batched input, with det(R) = +1.

Return type:

ndarray

Examples

>>> u = np.array([1.0, 2.0, 3.0])
>>> R = rot_axis_to_vec(u, 'x')
>>> bool(np.allclose(R @ [1, 0, 0], u / np.linalg.norm(u)))
True
>>> bool(np.allclose(R.T @ u, [np.linalg.norm(u), 0, 0]))
True

Notes

Port of rotAxis2Vec.m (method 0, the default), used in [1] for bistatic orbital debris tracking: a Householder reflection whose possible reflection component is removed by a mirror plus a 180-degree Givens rotation, which is numerically stabler than the shortest-rotation quaternion construction (method 1, not ported — no caller here needs it). MATLAB’s 1-based axis numbers and letters map to 0-based indices here.

References

pytcl.coordinate_systems.rotations.rotx(angle)[source]

Create rotation matrix for rotation about the x-axis.

Parameters:: angle (float) – Rotation angle in radians.
Returns:: R – 3x3 rotation matrix.
Return type:: ndarray

Examples

>>> R = rotx(np.pi/2)  # 90 degree rotation about x
>>> bool(np.allclose(R @ [0, 1, 0], [0, 0, 1]))  # y-axis maps to z-axis
True

pytcl.coordinate_systems.rotations.roty(angle)[source]

Create rotation matrix for rotation about the y-axis.

Parameters:: angle (float) – Rotation angle in radians.
Returns:: R – 3x3 rotation matrix.
Return type:: ndarray

Examples

>>> R = roty(np.pi/2)  # 90 degree rotation about y
>>> bool(np.allclose(R @ [1, 0, 0], [0, 0, -1]))  # x-axis maps to -z-axis
True

pytcl.coordinate_systems.rotations.rotz(angle)[source]

Create rotation matrix for rotation about the z-axis.

Parameters:: angle (float) – Rotation angle in radians.
Returns:: R – 3x3 rotation matrix.
Return type:: ndarray

Examples

>>> R = rotz(np.pi/2)  # 90 degree rotation about z
>>> bool(np.allclose(R @ [1, 0, 0], [0, 1, 0]))  # x-axis maps to y-axis
True

pytcl.coordinate_systems.rotations.euler2rotmat(angles, sequence='ZYX')[source]

Convert Euler angles to rotation matrix.

Parameters:

angles (array_like) – Three Euler angles in radians [angle1, angle2, angle3].
sequence (str, optional) – Rotation sequence (e.g., ‘ZYX’, ‘XYZ’, ‘ZXZ’). Default is ‘ZYX’ (aerospace convention: yaw-pitch-roll).

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

>>> yaw, pitch, roll = np.radians([45, 30, 15])
>>> R = euler2rotmat([yaw, pitch, roll], 'ZYX')

Notes

The rotation is applied in the order specified by the sequence, from right to left. For ‘ZYX’: R = Rz(yaw) @ Ry(pitch) @ Rx(roll).

pytcl.coordinate_systems.rotations.rotmat2euler(R, sequence='ZYX')[source]

Convert rotation matrix to Euler angles.

Parameters:

R (array_like) – 3x3 rotation matrix.
sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

angles – Three Euler angles in radians.

Return type:

ndarray

Examples

>>> R = rotz(np.radians(45)) @ roty(np.radians(30)) @ rotx(np.radians(15))
>>> angles = rotmat2euler(R, 'ZYX')
>>> np.degrees(angles)
array([45., 30., 15.])

Notes

May have singularities (gimbal lock) at certain angles.

pytcl.coordinate_systems.rotations.axisangle2rotmat(axis, angle)[source]

Convert axis-angle representation to rotation matrix.

Parameters:

axis (array_like) – Unit vector defining the rotation axis [ax, ay, az].
angle (float) – Rotation angle in radians.

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

>>> axis = [0, 0, 1]  # Z-axis
>>> R = axisangle2rotmat(axis, np.pi/2)  # 90 deg about Z
>>> bool(np.allclose(R @ [1, 0, 0], [0, 1, 0]))  # x-axis maps to y-axis
True

Notes

Uses Rodrigues’ rotation formula.

pytcl.coordinate_systems.rotations.rotmat2axisangle(R)[source]

Convert rotation matrix to axis-angle representation.

Parameters:

R (array_like) – 3x3 rotation matrix.

Returns:

axis (ndarray) – Unit vector rotation axis.
angle (float) – Rotation angle in radians [0, π].

Return type:

Tuple[ndarray[tuple[Any, …], dtype[floating]], float]

Examples

>>> R = rotz(np.pi/2)  # 90 deg about Z
>>> axis, angle = rotmat2axisangle(R)
>>> axis  # Z-axis
array([0., 0., 1.])
>>> np.degrees(angle)
90.0

pytcl.coordinate_systems.rotations.quat2rotmat(q)[source]

Convert quaternion to rotation matrix.

Parameters:: q (array_like) – Quaternion [qw, qx, qy, qz] (scalar-first convention).
Returns:: R – 3x3 rotation matrix.
Return type:: ndarray

Examples

>>> q = [1, 0, 0, 0]  # Identity quaternion
>>> quat2rotmat(q)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])

pytcl.coordinate_systems.rotations.rotmat2quat(R)[source]

Convert rotation matrix to quaternion.

Parameters:: R (array_like) – 3x3 rotation matrix.
Returns:: q – Quaternion [qw, qx, qy, qz] (scalar-first, positive qw).
Return type:: ndarray

Examples

>>> R = np.eye(3)  # Identity rotation
>>> rotmat2quat(R)
array([1., 0., 0., 0.])

Notes

Uses Shepperd’s method for numerical stability.

pytcl.coordinate_systems.rotations.euler2quat(angles, sequence='ZYX')[source]

Convert Euler angles to quaternion.

Parameters:

angles (array_like) – Three Euler angles in radians.
sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

q – Quaternion [qw, qx, qy, qz].

Return type:

ndarray

Examples

Convert yaw-pitch-roll angles to quaternion:

>>> import numpy as np
>>> from pytcl.coordinate_systems.rotations import euler2quat
>>> # 45° yaw, 30° pitch, 0° roll
>>> angles = np.radians([45, 30, 0])
>>> q = euler2quat(angles, sequence='ZYX')
>>> q.shape
(4,)
>>> np.abs(q[0]) > 0.5  # scalar part should be significant
True

See also

quat2euler: Inverse conversion.
euler2rotmat: Convert to rotation matrix instead.

pytcl.coordinate_systems.rotations.quat2euler(q, sequence='ZYX')[source]

Convert quaternion to Euler angles.

Parameters:

q (array_like) – Quaternion [qw, qx, qy, qz].
sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

angles – Three Euler angles in radians.

Return type:

ndarray

Examples

>>> q = [1, 0, 0, 0]  # Identity quaternion
>>> angles = quat2euler(q, 'ZYX')
>>> np.allclose(angles, [0, 0, 0])
True

pytcl.coordinate_systems.rotations.quat_multiply(q1, q2)[source]

Multiply two quaternions.

Parameters:

q1 (array_like) – First quaternion [qw, qx, qy, qz].
q2 (array_like) – Second quaternion.

Returns:

q – Product quaternion q1 * q2.

Return type:

ndarray

Notes

Quaternion multiplication represents composition of rotations. q1 * q2 applies q2 first, then q1.

Examples

Combine two rotations using quaternion multiplication:

>>> import numpy as np
>>> from pytcl.coordinate_systems.rotations import quat_multiply, euler2quat
>>> # 90° rotation about Z, then 45° about X
>>> q_z90 = euler2quat(np.radians([90, 0, 0]), 'ZYX')
>>> q_x45 = euler2quat(np.radians([0, 0, 45]), 'ZYX')
>>> q_combined = quat_multiply(q_z90, q_x45)
>>> q_combined.shape
(4,)
>>> np.isclose(np.linalg.norm(q_combined), 1.0)  # unit quaternion
True

See also

quat_inverse: Compute quaternion inverse.
quat_rotate: Rotate a vector by a quaternion.

pytcl.coordinate_systems.rotations.quat_conjugate(q)[source]

Compute quaternion conjugate.

Parameters:: q (array_like) – Quaternion [qw, qx, qy, qz].
Returns:: q_conj – Conjugate quaternion [qw, -qx, -qy, -qz].
Return type:: ndarray

Examples

>>> quat_conjugate([0.707, 0.707, 0, 0])
array([ 0.707, -0.707, -0.   , -0.   ])

pytcl.coordinate_systems.rotations.quat_inverse(q)[source]

Compute quaternion inverse.

Parameters:: q (array_like) – Quaternion [qw, qx, qy, qz].
Returns:: q_inv – Inverse quaternion.
Return type:: ndarray

Examples

>>> q = euler2quat(np.radians([45, 0, 0]), 'ZYX')
>>> q_inv = quat_inverse(q)
>>> quat_multiply(q, q_inv)  # Should be identity
array([1., 0., 0., 0.])

Notes

For unit quaternions, inverse equals conjugate.

pytcl.coordinate_systems.rotations.quat_rotate(q, v)[source]

Rotate a vector using a quaternion.

Parameters:

q (array_like) – Quaternion [qw, qx, qy, qz].
v (array_like) – Vector to rotate [x, y, z].

Returns:

v_rot – Rotated vector.

Return type:

ndarray

Examples

>>> # 90 degree rotation about z-axis
>>> q = euler2quat(np.radians([90, 0, 0]), 'ZYX')
>>> v = np.array([1.0, 0.0, 0.0])
>>> v_rot = quat_rotate(q, v)
>>> bool(np.allclose(v_rot, [0, 1, 0]))  # x-axis becomes y-axis
True

Notes

Computes q * v * q^(-1) where v is treated as a pure quaternion.

pytcl.coordinate_systems.rotations.slerp(q1, q2, t)[source]

Spherical linear interpolation between two quaternions.

Parameters:

q1 (array_like) – Start quaternion.
q2 (array_like) – End quaternion.
t (float) – Interpolation parameter in [0, 1].

Returns:

q – Interpolated quaternion.

Return type:

ndarray

Examples

>>> q1 = np.array([1, 0, 0, 0])  # identity
>>> q2 = euler2quat(np.radians([90, 0, 0]), 'ZYX')  # 90 deg about z
>>> q_mid = slerp(q1, q2, 0.5)  # halfway = 45 deg
>>> angles = quat2euler(q_mid, 'ZYX')
>>> np.degrees(angles[0])  # yaw should be ~45
45.0...

pytcl.coordinate_systems.rotations.rodrigues2rotmat(rvec)[source]

Convert Rodrigues vector to rotation matrix.

Parameters:: rvec (array_like) – Rodrigues vector (axis * angle).
Returns:: R – 3x3 rotation matrix.
Return type:: ndarray

Examples

>>> rvec = [0, 0, np.pi/2]  # 90 deg about Z
>>> R = rodrigues2rotmat(rvec)
>>> bool(np.allclose(R @ [1, 0, 0], [0, 1, 0]))  # x-axis maps to y-axis
True

Notes

The Rodrigues vector encodes both the rotation axis and angle:

rvec = axis * angle, where |rvec| = angle

pytcl.coordinate_systems.rotations.rotmat2rodrigues(R)[source]

Convert rotation matrix to Rodrigues vector.

Parameters:: R (array_like) – 3x3 rotation matrix.
Returns:: rvec – Rodrigues vector (axis * angle).
Return type:: ndarray

Examples

>>> R = rotz(np.pi/2)  # 90 deg about Z
>>> rvec = rotmat2rodrigues(R)
>>> np.linalg.norm(rvec)  # magnitude is the angle
1.5707...

pytcl.coordinate_systems.rotations.dcm_rate(R, omega)[source]

Compute the time derivative of a rotation matrix.

Parameters:

R (array_like) – Current rotation matrix.
omega (array_like) – Angular velocity vector [wx, wy, wz] in body frame.

Returns:

R_dot – Time derivative of R.

Return type:

ndarray

Examples

>>> R = np.eye(3)
>>> omega = [0, 0, 1]  # 1 rad/s about Z
>>> R_dot = dcm_rate(R, omega)
>>> R_dot[0, 1]  # Off-diagonal elements show rotation
-1.0

Notes

R_dot = R @ skew(omega)

pytcl.coordinate_systems.rotations.is_rotation_matrix(R, tol=1e-06)[source]

Check if a matrix is a valid rotation matrix.

Parameters:

R (array_like) – Matrix to check.
tol (float, optional) – Tolerance for numerical checks.

Returns:

valid – True if R is a valid rotation matrix.

Return type:

bool

Examples

>>> R = rotx(np.pi/4)
>>> is_rotation_matrix(R)
True
>>> is_rotation_matrix(np.eye(3) * 2)  # not orthonormal
False

Jacobians

Jacobian matrices for coordinate transformations.

This module provides: - Jacobians for spherical/Cartesian transformations - Jacobians for polar transformations - Jacobians for r-u-v direction cosines - Jacobians for ECEF/ENU/NED transformations - Jacobians for geodetic transformations - Covariance transformation utilities

pytcl.coordinate_systems.jacobians.spherical_jacobian(cart_point, system_type='az-el')[source]

Compute Jacobian of Cartesian to spherical transformation.

Returns the Jacobian matrix J where d[r, az, el] = J @ d[x, y, z].

Parameters:

cart_point (array_like) – Cartesian coordinates [x, y, z].
system_type ({'standard', 'az-el'}, optional) – Spherical coordinate convention. Default is ‘az-el’.

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

Notes

For tracking convention (‘az-el’): - r = sqrt(x² + y² + z²) - az = atan2(y, x) - el = atan2(z, sqrt(x² + y²))

Examples

>>> J = spherical_jacobian([1, 1, 1])
>>> J.shape
(3, 3)

pytcl.coordinate_systems.jacobians.spherical_jacobian_inv(r, az, el, system_type='az-el')[source]

Compute Jacobian of spherical to Cartesian transformation.

Returns the Jacobian matrix J where d[x, y, z] = J @ d[r, az, el].

Parameters:

r (float) – Range.
az (float) – Azimuth in radians.
el (float) – Elevation in radians.
system_type ({'standard', 'az-el'}, optional) – Spherical coordinate convention.

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

pytcl.coordinate_systems.jacobians.polar_jacobian(cart_point)[source]

Compute Jacobian of 2D Cartesian to polar transformation.

Returns J where d[r, theta] = J @ d[x, y].

Parameters:: cart_point (array_like) – Cartesian coordinates [x, y].
Returns:: J – 2x2 Jacobian matrix.
Return type:: ndarray

pytcl.coordinate_systems.jacobians.polar_jacobian_inv(r, theta)[source]

Compute Jacobian of polar to 2D Cartesian transformation.

Returns J where d[x, y] = J @ d[r, theta].

Parameters:

r (float) – Radial distance.
theta (float) – Angle in radians.

Returns:

J – 2x2 Jacobian matrix.

Return type:

ndarray

pytcl.coordinate_systems.jacobians.ruv_jacobian(cart_point)[source]

Compute Jacobian of Cartesian to r-u-v transformation.

The r-u-v representation uses range and direction cosines: - u = x/r (direction cosine along x) - v = y/r (direction cosine along y)

Parameters:: cart_point (array_like) – Cartesian coordinates [x, y, z].
Returns:: J – 3x3 Jacobian matrix where d[r, u, v] = J @ d[x, y, z].
Return type:: ndarray

pytcl.coordinate_systems.jacobians.enu_jacobian(lat, lon)[source]

Compute Jacobian of ECEF to ENU transformation.

Returns J where d[e, n, u] = J @ d[x, y, z].

Parameters:

lat (float) – Reference latitude in radians.
lon (float) – Reference longitude in radians.

Returns:

J – 3x3 rotation matrix (Jacobian is constant for this linear transformation).

Return type:

ndarray

Notes

Uses cached computation with quantized inputs for performance.

pytcl.coordinate_systems.jacobians.ned_jacobian(lat, lon)[source]

Compute Jacobian of ECEF to NED transformation.

Parameters:

lat (float) – Reference latitude in radians.
lon (float) – Reference longitude in radians.

Returns:

J – 3x3 rotation matrix.

Return type:

ndarray

Notes

Uses cached computation with quantized inputs for performance.

pytcl.coordinate_systems.jacobians.geodetic_jacobian(lat, lon, alt, a=6378137.0, f=0.0033528106647474805)[source]

Compute Jacobian of geodetic to ECEF transformation.

Returns J where d[x, y, z] = J @ d[lat, lon, alt].

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Geodetic longitude in radians.
alt (float) – Altitude above ellipsoid in meters.
a (float, optional) – Semi-major axis (default: WGS84).
f (float, optional) – Flattening (default: WGS84).

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

pytcl.coordinate_systems.jacobians.cross_covariance_transform(J, P)[source]

Transform a covariance matrix through a Jacobian.

Computes P_new = J @ P @ J.T for error propagation.

Parameters:

J (array_like) – Jacobian matrix of the transformation.
P (array_like) – Original covariance matrix.

Returns:

P_new – Transformed covariance matrix.

Return type:

ndarray

Examples

>>> # Transform spherical covariance to Cartesian
>>> P_sph = np.diag([1, 0.01, 0.01])  # [r, az, el] variances
>>> r, az, el = 1000, np.radians(45), np.radians(30)
>>> J = spherical_jacobian_inv(r, az, el)
>>> P_cart = cross_covariance_transform(J, P_sph)

pytcl.coordinate_systems.jacobians.numerical_jacobian(func, x, dx=1e-07)[source]

Compute Jacobian numerically using central differences.

Parameters:

func (callable) – Function f(x) -> y.
x (array_like) – Point at which to compute Jacobian.
dx (float, optional) – Step size for finite differences.

Returns:

J – Jacobian matrix.

Return type:

ndarray

Projections

Map Projections for Tracking and Navigation.

This module provides map projection functions commonly used in tracking, navigation, and geospatial applications. All projections convert between geodetic coordinates (latitude, longitude) and planar map coordinates (x, y).

Projections included: - Mercator: Conformal cylindrical projection, good for small-scale maps - Transverse Mercator: Conformal projection for UTM zones - UTM: Universal Transverse Mercator with zone handling - Stereographic: Conformal azimuthal projection, good for polar regions - Lambert Conformal Conic: Conformal conic projection for mid-latitudes - Azimuthal Equidistant: Preserves distances from center point

All angles are in radians unless otherwise noted.

References

Snyder, J. P. “Map Projections: A Working Manual.” U.S. Geological Survey Professional Paper 1395, 1987.
NIMA Technical Report 8350.2, “Department of Defense World Geodetic System 1984,” Third Edition, 2000.
Karney, C. F. F. “Transverse Mercator with an accuracy of a few nanometers.” Journal of Geodesy 85.8 (2011): 475-485.

class pytcl.coordinate_systems.projections.projections.ProjectionResult(x, y, scale, convergence)[source]

Bases: NamedTuple

Result of a map projection.

Parameters:

x (float) – Easting coordinate in meters.
y (float) – Northing coordinate in meters.
scale (float) – Scale factor at the point (ratio of map distance to ground distance).
convergence (float) – Grid convergence angle in radians (angle from true north to grid north).

x: float: Alias for field number 0

y: float: Alias for field number 1

scale: float: Alias for field number 2

convergence: float: Alias for field number 3

class pytcl.coordinate_systems.projections.projections.UTMResult(easting, northing, zone, hemisphere, scale, convergence)[source]

Bases: NamedTuple

Result of UTM projection.

Parameters:

easting (float) – UTM easting coordinate in meters.
northing (float) – UTM northing coordinate in meters.
zone (int) – UTM zone number (1-60).
hemisphere (str) – ‘N’ for northern hemisphere, ‘S’ for southern.
scale (float) – Scale factor at the point.
convergence (float) – Grid convergence angle in radians.

easting: float: Alias for field number 0

northing: float: Alias for field number 1

zone: int: Alias for field number 2

hemisphere: str: Alias for field number 3

scale: float: Alias for field number 4

convergence: float: Alias for field number 5

pytcl.coordinate_systems.projections.projections.mercator(lat, lon, lon0=0.0, a=6378137.0, e=np.float64(0.08181919084262149))[source]

Ellipsoidal Mercator projection (forward).

The Mercator projection is a conformal cylindrical projection where rhumb lines (lines of constant bearing) appear as straight lines. Scale increases with latitude, becoming infinite at the poles.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
lon0 (float, optional) – Central meridian in radians. Default is 0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e (float, optional) – First eccentricity. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

The Mercator projection is conformal, meaning it preserves local angles. However, it distorts areas, especially at high latitudes.

Examples

>>> import numpy as np
>>> result = mercator(np.radians(45), np.radians(-75))
>>> print(f"x={result.x:.1f}, y={result.y:.1f}")
x=-8348961.8, y=5591295.9

pytcl.coordinate_systems.projections.projections.mercator_inverse(x, y, lon0=0.0, a=6378137.0, e=np.float64(0.08181919084262149), tol=1e-12, max_iter=10)[source]

Ellipsoidal Mercator projection (inverse).

Parameters:

x (float) – Easting coordinate in meters.
y (float) – Northing coordinate in meters.
lon0 (float, optional) – Central meridian in radians. Default is 0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e (float, optional) – First eccentricity. Default is WGS84.
tol (float, optional) – Convergence tolerance. Default is 1e-12.
max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

Examples

>>> import numpy as np
>>> lat, lon = mercator_inverse(1000000, 5000000)

pytcl.coordinate_systems.projections.projections.transverse_mercator(lat, lon, lat0=0.0, lon0=0.0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Transverse Mercator projection (forward).

The Transverse Mercator is a conformal projection that rotates the cylinder to be tangent along a meridian. It’s the basis for UTM.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
lat0 (float, optional) – Origin latitude in radians. Default is 0.
lon0 (float, optional) – Central meridian in radians. Default is 0.
k0 (float, optional) – Scale factor at central meridian. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

This implementation uses the Redfearn series expansion. Measured against PROJ’s +proj=tmerc on WGS84, worst case over latitudes 0-75 degrees:

From the meridian	Difference
2 degrees	< 0.1 mm
3 degrees	0.1 mm
4 degrees	0.7 mm
5 degrees	3.3 mm
6 degrees	11.5 mm

The docstring previously said “about 1 meter within 4 degrees”, which understated it by three orders of magnitude and could have led a caller to reject the function as too coarse (gh-25 follow-up).

Examples

>>> import numpy as np
>>> result = transverse_mercator(np.radians(45), np.radians(-75),
...                              lon0=np.radians(-75))

pytcl.coordinate_systems.projections.projections.transverse_mercator_inverse(x, y, lat0=0.0, lon0=0.0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Transverse Mercator projection (inverse).

Parameters:

x (float) – Easting coordinate in meters.
y (float) – Northing coordinate in meters.
lat0 (float, optional) – Origin latitude in radians. Default is 0.
lon0 (float, optional) – Central meridian in radians. Default is 0.
k0 (float, optional) – Scale factor at central meridian. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

pytcl.coordinate_systems.projections.projections.utm_zone(lon, lat=0.0)[source]

Determine UTM zone number from longitude.

Parameters:

lon (float) – Longitude in radians.
lat (float, optional) – Latitude in radians (used for Norway/Svalbard exceptions).

Returns:

UTM zone number (1-60).

Return type:

int

Notes

Standard zones are 6 degrees wide. Special zones exist for Norway and Svalbard.

pytcl.coordinate_systems.projections.projections.utm_central_meridian(zone)[source]

Get central meridian for UTM zone.

Parameters:: zone (int) – UTM zone number (1-60).
Returns:: Central meridian in radians.
Return type:: float

pytcl.coordinate_systems.projections.projections.geodetic2utm(lat, lon, zone=None)[source]

Convert geodetic coordinates to UTM.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
zone (int, optional) – Force specific UTM zone. If None, computed from longitude.

Returns:

UTM coordinates with zone information.

Return type:

UTMResult

Examples

>>> import numpy as np
>>> result = geodetic2utm(np.radians(45.0), np.radians(-75.5))
>>> print(f"Zone {result.zone}{result.hemisphere}: "
...       f"E={result.easting:.1f}, N={result.northing:.1f}")
Zone 18N: E=460592.4, N=4983072.0

pytcl.coordinate_systems.projections.projections.utm2geodetic(easting, northing, zone, hemisphere='N')[source]

Convert UTM coordinates to geodetic.

Parameters:

easting (float) – UTM easting in meters.
northing (float) – UTM northing in meters.
zone (int) – UTM zone number (1-60).
hemisphere (str, optional) – ‘N’ for northern, ‘S’ for southern hemisphere.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

Examples

>>> lat, lon = utm2geodetic(500000, 5000000, 18, 'N')
>>> print(f"Lat: {np.degrees(lat):.4f}, Lon: {np.degrees(lon):.4f}")
Lat: 45.1535, Lon: -75.0000

pytcl.coordinate_systems.projections.projections.geodetic2utm_batch(lats, lons, zone=None)[source]

Batch convert geodetic coordinates to UTM.

Parameters:

lats (ndarray) – Geodetic latitudes in radians.
lons (ndarray) – Longitudes in radians.
zone (int, optional) – Force specific UTM zone. If None, computed per point.

Returns:

eastings (ndarray) – UTM eastings in meters.
northings (ndarray) – UTM northings in meters.
zones (ndarray) – UTM zone numbers.
hemispheres (ndarray) – Hemisphere characters (‘N’ or ‘S’).

Return type:

tuple[ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[floating]], ndarray[tuple[Any, …], dtype[int64]], ndarray[tuple[Any, …], dtype[Any]]]

pytcl.coordinate_systems.projections.projections.stereographic(lat, lon, lat0, lon0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Oblique stereographic projection (forward).

The stereographic projection is conformal and azimuthal. It’s commonly used for polar regions and local surveys.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
lat0 (float) – Center latitude in radians.
lon0 (float) – Center longitude in radians.
k0 (float, optional) – Scale factor at center. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

For polar stereographic, use lat0 = +-pi/2.

Examples

>>> import numpy as np
>>> # Polar stereographic centered at North Pole
>>> result = stereographic(np.radians(85), np.radians(45),
...                        np.radians(90), 0)

Notes

The oblique case maps onto a conformal sphere of Gaussian radius, sqrt(M0 * N0) at the origin latitude. PROJ’s +proj=sterea uses a sphere of radius N0 instead, so the two disagree away from the center (gh-25). How much depends strongly on the origin latitude, because the two radii coincide where the meridional and prime-vertical curvatures do. Measured against +proj=sterea on WGS84:

`lat0`	At 400 km	At 1,738 km
0 deg	1.34 km	5.9 km
30 deg	0.03 km	0.73 km
45 deg	1.38 km	6.7 km
60 deg	2.73 km	12.5 km
75 deg	3.71 km	16.5 km

An earlier version of this note quoted a single pair of figures, 2.5 km and 8.7 km, as though they were properties of the projection. They are one origin latitude: the disagreement varies eightyfold across the table above.

Both are conformal and each is self-consistent with its own inverse. But coordinates from this function are not interchangeable with +proj=sterea output, and it is not survey-grade far from the projection center. For work that must match PROJ, use pyproj.

pytcl.coordinate_systems.projections.projections.stereographic_inverse(x, y, lat0, lon0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316), tol=1e-12, max_iter=10)[source]

Oblique stereographic projection (inverse).

Parameters:

x (float) – Easting in meters.
y (float) – Northing in meters.
lat0 (float) – Center latitude in radians.
lon0 (float) – Center longitude in radians.
k0 (float, optional) – Scale factor at center. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.
tol (float, optional) – Convergence tolerance. Default is 1e-12.
max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

pytcl.coordinate_systems.projections.projections.polar_stereographic(lat, lon, north=True, k0=0.994, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Polar stereographic projection (forward).

Standard polar stereographic used for polar regions.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
north (bool, optional) – True for North Pole center, False for South Pole. Default is True.
k0 (float, optional) – Scale factor at pole. Default is 0.994 (UPS standard).
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

Uses the exact ellipsoidal polar stereographic equations (Snyder 1987, eqs. 15-9 and 21-33/21-34), matching standard UPS/EPSG implementations.

Examples

>>> import numpy as np
>>> # Arctic location
>>> result = polar_stereographic(np.radians(80), np.radians(45))

pytcl.coordinate_systems.projections.projections.lambert_conformal_conic(lat, lon, lat0, lon0, lat1, lat2, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Lambert Conformal Conic projection (forward).

A conformal conic projection with two standard parallels where scale is exact. Good for mid-latitude regions with east-west extent.

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
lat0 (float) – Origin latitude in radians.
lon0 (float) – Central meridian in radians.
lat1 (float) – First standard parallel in radians.
lat2 (float) – Second standard parallel in radians.
k0 (float, optional) – Scale factor. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Examples

>>> import numpy as np
>>> # Continental US projection
>>> result = lambert_conformal_conic(
...     np.radians(40), np.radians(-100),
...     lat0=np.radians(39), lon0=np.radians(-96),
...     lat1=np.radians(33), lat2=np.radians(45)
... )

pytcl.coordinate_systems.projections.projections.lambert_conformal_conic_inverse(x, y, lat0, lon0, lat1, lat2, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316), tol=1e-12, max_iter=10)[source]

Lambert Conformal Conic projection (inverse).

Parameters:

x (float) – Easting in meters.
y (float) – Northing in meters.
lat0 (float) – Origin latitude in radians.
lon0 (float) – Central meridian in radians.
lat1 (float) – First standard parallel in radians.
lat2 (float) – Second standard parallel in radians.
k0 (float, optional) – Scale factor. Default is 1.0.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.
tol (float, optional) – Convergence tolerance. Default is 1e-12.
max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

pytcl.coordinate_systems.projections.projections.azimuthal_equidistant(lat, lon, lat0, lon0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Azimuthal equidistant projection (forward).

Distances from the center point are preserved. Useful for showing distances from a specific location (e.g., radio coverage).

Parameters:

lat (float) – Geodetic latitude in radians.
lon (float) – Longitude in radians.
lat0 (float) – Center latitude in radians.
lon0 (float) – Center longitude in radians.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

This implementation uses the spherical approximation for simplicity. For high accuracy over long distances, geodesic calculations should be used.

This projection is not conformal, so there is no single scale factor. The radial scale is exactly 1 by construction and is what ProjectionResult.scale reports; the tangential scale is c / sin(c) for angular distance c from the centre, reaching about 1.005 at 1,000 km and 1.11 at 5,000 km. convergence accounts for this and is the true grid convergence, not the conformal approximation.

Examples

>>> import numpy as np
>>> result = azimuthal_equidistant(np.radians(40), np.radians(-75),
...                                np.radians(38.9), np.radians(-77))

pytcl.coordinate_systems.projections.projections.azimuthal_equidistant_inverse(x, y, lat0, lon0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Azimuthal equidistant projection (inverse).

Parameters:

x (float) – Easting in meters.
y (float) – Northing in meters.
lat0 (float) – Center latitude in radians.
lon0 (float) – Center longitude in radians.
a (float, optional) – Semi-major axis in meters. Default is WGS84.
e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]