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')

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')

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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')

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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)

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

[1]

D. F. Crouse, “On measurement-based light-time corrections for bistatic orbital debris tracking,” IEEE Transactions on Aerospace and Electronic Systems, vol. 51, no. 3, pp. 2502-2518, Jul. 2015.

pytcl.coordinate_systems.rotations.rotx(angle)

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)

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)

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')

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')

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)

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)

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)

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)

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')

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')

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)

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)

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)

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)

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)

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)

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)

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)

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)

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

Coordinate-transformation Jacobians, the component gradients of bistatic measurement functions (range, range rate, spherical/polar angles, direction cosines, TDOA) and the full measurement Jacobians (spherical, polar, r-u-v, with range-rate and converted variants) composed from them.

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.range_gradient(x, use_half_range=False, l_tx=None, l_rx=None)

Gradient of bistatic range with respect to Cartesian position.

A transmitter exactly collocated with the target is assumed to move with it, dropping the transmitter term.

Port of rangeGradient.

Parameters:

• x (array_like) – (d,) or (d, N) target positions, d in {1, 2, 3}.

• use_half_range (bool, optional) – Divide the bistatic range by two. Default False.

• l_tx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• l_rx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

Returns:

J – (1, d) gradient, or a (1, d, N) stack.

Return type:

ndarray

Examples

>>> range_gradient([3.0, 4.0])
array([[1.2, 1.6]])

pytcl.coordinate_systems.jacobians.range_rate_gradient(x, use_half_range=False, l_tx=None, l_rx=None)

Gradient of bistatic range rate with respect to the target state.

The state stacks position over velocity (2d elements, d in {1, 2, 3}); l_tx/l_rx are transmitter/receiver states of the same layout, so moving platforms are handled. A transmitter collocated with the target moves with it.

Port of rangeRateGradient.

Parameters:

• x (array_like) – (2d,) or (2d, N) target states.

• use_half_range (bool, optional) – Divide by two for one-way rates. Default False.

• l_tx (array_like, optional) – Transmitter/receiver states. Default: the origin, at rest.

• l_rx (array_like, optional) – Transmitter/receiver states. Default: the origin, at rest.

Returns:

J – (1, 2d) gradient, or a (1, 2d, N) stack.

Return type:

ndarray

Examples

>>> J = range_rate_gradient([0.0, 5e3, 100.0, 0.0])
>>> J.shape
(1, 4)

pytcl.coordinate_systems.jacobians.spher_ang_gradient(x, system_type=0, l_rx=None, m=None)

Gradient of spherical azimuth and elevation w.r.t. position.

Port of spherAngGradient.

Parameters:

• x (array_like) – (3,) or (3, N) global Cartesian positions.

• system_type (int, optional) – Angle convention 0-3 (see the module docstring). Default 0.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

Returns:

J – (2, 3) [d(az); d(el)] rows, or a (2, 3, N) stack.

Return type:

ndarray

Examples

>>> J = spher_ang_gradient([1e3, 2e3, 500.0])
>>> J.shape
(2, 3)

pytcl.coordinate_systems.jacobians.pol_ang_gradient(x, system_type=0, l_rx=None)

Gradient of a 2-D polar angle with respect to position.

system_type 0 measures counterclockwise from +x, 1 clockwise from +y.

Port of polAngGradient.

Parameters:

• x (array_like) – (2,) or (2, N) positions.

• system_type (int, optional) – Angle convention (0 or 1). Default 0.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

Returns:

J – (1, 2) gradient, or a (1, 2, N) stack.

Return type:

ndarray

Examples

>>> pol_ang_gradient([2.0, 0.0])
array([[-0. ,  0.5]])

pytcl.coordinate_systems.jacobians.uv_gradient(x, l_rx=None, m=None, include_w=False)

Gradient of u-v(-w) direction cosines with respect to position.

Port of uvGradient.

Parameters:

• x (array_like) – (3,) or (3, N) global Cartesian positions.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

• include_w (bool, optional) – Also return the gradient of the third direction cosine. Default False.

Returns:

J – (2, 3) [du; dv] rows ((3, 3) with include_w), or the corresponding stack.

Return type:

ndarray

Examples

>>> uv_gradient([0.0, 0.0, 10.0])
array([[0.1, 0. , 0. ],
       [0. , 0.1, 0. ]])

pytcl.coordinate_systems.jacobians.u_gradient_2d(x, l_rx=None, m=None, include_v=False)

Gradient of the 2-D direction cosine u (optionally v) w.r.t. position.

Port of uGradient2D.

Parameters:

• x (array_like) – (2,) or (2, N) positions.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (2, 2) global-to-local rotation. Default: identity.

• include_v (bool, optional) – Also return the second direction cosine’s gradient. Default False.

Returns:

J – (1, 2) ((2, 2) with include_v), or the corresponding stack.

Return type:

ndarray

Examples

>>> u_gradient_2d([0.0, 4.0])
array([[0.25, 0.  ]])

pytcl.coordinate_systems.jacobians.u_gradient_3d(x, l_rx=None, m=None)

Gradient of the first 3-D direction cosine u w.r.t. position.

Port of uGradient3D.

Parameters:

• x (array_like) – (3,) or (3, N) positions.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

Returns:

J – (1, 3) gradient, or a (1, 3, N) stack.

Return type:

ndarray

Examples

>>> u_gradient_3d([0.0, 5.0, 0.0])
array([[0.2, 0. , 0. ]])

pytcl.coordinate_systems.jacobians.tdoa_gradient(x, l_ref=None, l_rx=None, c=299792458.0)

Gradient of a TDOA measurement with respect to position.

The time difference of arrival between the receiver at l_rx and the reference receiver at l_ref, in seconds.

Port of TDOAGradient.

Parameters:

• x (array_like) – (d,) or (d, N) emitter positions.

• l_ref (array_like, optional) – Reference receiver position. Default: the origin.

• l_rx (array_like, optional) – Second receiver position. Default: the origin.

• c (float, optional) – Propagation speed. Default: the vacuum speed of light.

Returns:

J – (1, d) gradient, or a (1, d, N) stack.

Return type:

ndarray

Examples

>>> J = tdoa_gradient([1e3, 2e3], l_ref=[0.0, 0.0], l_rx=[5e3, 0.0])
>>> J.shape
(1, 2)

pytcl.coordinate_systems.jacobians.calc_spher_jacob(x, system_type=0, use_half_range=None, l_tx=None, l_rx=None, m=None)

Jacobian of a bistatic [range; azimuth; elevation] measurement.

Port of calcSpherJacob.

Parameters:

• x (array_like) – (3,) target position in global Cartesian coordinates.

• system_type (int, optional) – Spherical angle convention 0-3. Default 0.

• use_half_range (bool, optional) – One-way range convention. Default: True when no transmitter is given (monostatic), False otherwise – the MATLAB default.

• l_tx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• l_rx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local receiver rotation. Default: identity.

Returns:

J – (3, 3) Jacobian; rows d(r), d(az), d(el).

Return type:

ndarray

Examples

>>> J = calc_spher_jacob([1e3, 2e3, 500.0])
>>> J.shape
(3, 3)

pytcl.coordinate_systems.jacobians.calc_spher_inv_jacob(z, system_type=0)

Jacobian of Cartesian position w.r.t. a monostatic spherical triple.

Port of calcSpherInvJacob.

Parameters:

• z (array_like) – (3,) measurement [r, az, el].

• system_type (int, optional) – Spherical angle convention 0-3. Default 0.

Returns:

J – (3, 3) Jacobian; columns d/dr, d/d(az), d/d(el).

Return type:

ndarray

Examples

>>> J = calc_spher_inv_jacob([1.0, 0.0, 0.0])
>>> float(J[0, 0])
1.0

pytcl.coordinate_systems.jacobians.calc_spher_rr_jacob(x, system_type=0, use_half_range=None, l_tx=None, l_rx=None, m=None)

Jacobian of [range; azimuth; elevation; range rate] w.r.t. state.

The state stacks position over velocity; l_tx/l_rx are six-element states so moving platforms are handled.

Port of calcSpherRRJacob.

Parameters:

• x (array_like) – (6,) target state.

• system_type (int) – As in calc_spher_jacob(), with states for the platforms.

• use_half_range (bool | None) – As in calc_spher_jacob(), with states for the platforms.

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob(), with states for the platforms.

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob(), with states for the platforms.

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob(), with states for the platforms.

Returns:

J – (4, 6) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_spher_rr_jacob([1e3, 2e3, 500.0, 10.0, -20.0, 5.0])
>>> J.shape
(4, 6)

pytcl.coordinate_systems.jacobians.calc_polar_jacob(x, system_type=0, use_half_range=True, l_tx=None, l_rx=None)

Jacobian of a 2-D bistatic [range; angle] measurement.

Port of calcPolarJacob.

Parameters:

• x (array_like) – (2,) target position.

• system_type (int, optional) – Polar angle convention (0 or 1). Default 0.

• use_half_range (bool, optional) – One-way range convention. Default True (the MATLAB default).

• l_tx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• l_rx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

Returns:

J – (2, 2) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_polar_jacob([3.0, 4.0])
>>> J.shape
(2, 2)

pytcl.coordinate_systems.jacobians.calc_polar_rr_jacob(x, system_type=0, use_half_range=True, l_tx=None, l_rx=None)

Jacobian of a 2-D [range; angle; range rate] w.r.t. the state.

Port of calcPolarRRJacob.

Parameters:

• x (array_like) – (4,) target state [x, y, vx, vy].

• system_type (int) – As in calc_polar_jacob(), with four-element states for the platforms.

• use_half_range (bool) – As in calc_polar_jacob(), with four-element states for the platforms.

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_jacob(), with four-element states for the platforms.

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_jacob(), with four-element states for the platforms.

Returns:

J – (3, 4) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_polar_rr_jacob([3e3, 4e3, 10.0, -5.0])
>>> J.shape
(3, 4)

pytcl.coordinate_systems.jacobians.calc_ruv_jacob(x, use_half_range=False, l_tx=None, l_rx=None, m=None, include_w=False)

Jacobian of a bistatic r-u-v(-w) measurement w.r.t. position.

Port of calcRuvJacob.

Parameters:

• x (array_like) – (3,) target position.

• use_half_range (bool, optional) – One-way range convention. Default False (the MATLAB default).

• l_tx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• l_rx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local receiver rotation. Default: identity.

• include_w (bool, optional) – Include the redundant third direction cosine row. Default False.

Returns:

J – (3, 3) Jacobian ((4, 3) with include_w).

Return type:

ndarray

Examples

>>> J = calc_ruv_jacob([1e3, 2e3, 5e3])
>>> J.shape
(3, 3)

pytcl.coordinate_systems.jacobians.calc_ruv_rr_jacob(x, use_half_range=False, l_tx=None, l_rx=None, m=None, include_w=False)

Jacobian of [r; u; v(; w); range rate] w.r.t. the target state.

Port of calcRuvRRJacob.

Parameters:

• x (array_like) – (6,) target state.

• use_half_range (bool) – As in calc_ruv_jacob(), with six-element states for the platforms.

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob(), with six-element states for the platforms.

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob(), with six-element states for the platforms.

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob(), with six-element states for the platforms.

• include_w (bool) – As in calc_ruv_jacob(), with six-element states for the platforms.

Returns:

J – (4, 6) Jacobian ((5, 6) with include_w).

Return type:

ndarray

Examples

>>> J = calc_ruv_rr_jacob([1e3, 2e3, 5e3, 10.0, -20.0, 5.0])
>>> J.shape
(4, 6)

pytcl.coordinate_systems.jacobians.calc_cart_rr_jacob(x, components=0, use_half_range=False, l_tx=None, l_rx=None)

Jacobian of a Cartesian-position(-plus-range-rate) measurement.

With components=0 the measurement is [position; range rate]; any other value drops the range-rate row (position only), matching the MATLAB semantics.

Port of calcCartRRJacob.

Parameters:

• x (array_like) – (2d,) target state, d in {1, 2, 3}.

• components (int, optional) – 0 to include the range-rate row. Default 0.

• use_half_range (bool) – As in range_rate_gradient().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in range_rate_gradient().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in range_rate_gradient().

Returns:

J – (d+1, 2d) Jacobian ((d, 2d) without range rate).

Return type:

ndarray

Examples

>>> J = calc_cart_rr_jacob([1e3, 2e3, 5e3, 10.0, -20.0, 5.0])
>>> J.shape
(4, 6)

pytcl.coordinate_systems.jacobians.calc_spher_conv_jacob(z, system_type=0, use_half_range=None, l_tx=None, l_rx=None, m=None)

Spherical measurement Jacobian evaluated at the measurement itself.

The measurement is inverted to Cartesian and the Jacobian of calc_spher_jacob() is evaluated there – the form needed to transform a measured covariance.

Port of calcSpherConvJacob.

Parameters:

• z (array_like) – (3,) measurement [r, az, el].

• system_type (int) – As in calc_spher_jacob().

• use_half_range (bool | None) – As in calc_spher_jacob().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob().

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob().

Returns:

J – (3, 3) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_spher_conv_jacob([1e3, 0.5, 0.2])
>>> J.shape
(3, 3)

pytcl.coordinate_systems.jacobians.calc_polar_conv_jacob(z, system_type=0, use_half_range=True, l_tx=None, l_rx=None, m=None)

Polar measurement Jacobian evaluated at the measurement itself.

Port of calcPolarConvJacob.

Parameters:

• z (array_like) – (2,) measurement [r, angle].

• system_type (int) – As in calc_polar_jacob().

• use_half_range (bool) – As in calc_polar_jacob().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_jacob().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_jacob().

• m (array_like, optional) – (2, 2) global-to-local receiver rotation used in the measurement inversion. Default: identity.

Returns:

J – (2, 2) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_polar_conv_jacob([2e3, 0.7])
>>> J.shape
(2, 2)

pytcl.coordinate_systems.jacobians.calc_polar_rr_conv_jacob(z, system_type=0, use_half_range=True, l_tx=None, l_rx=None, m=None)

Converted-polar Jacobian with a pass-through range-rate row.

Port of calcPolarRRConvJacob.

Parameters:

• z (array_like) – (3,) measurement [r, angle, range rate].

• system_type (int) – As in calc_polar_conv_jacob().

• use_half_range (bool) – As in calc_polar_conv_jacob().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_conv_jacob().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_conv_jacob().

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_polar_conv_jacob().

Returns:

J – (3, 3) block-diagonal Jacobian.

Return type:

ndarray

Examples

>>> J = calc_polar_rr_conv_jacob([2e3, 0.7, 15.0])
>>> float(J[2, 2])
1.0

pytcl.coordinate_systems.jacobians.calc_ruv_conv_jacob(z, use_half_range=False, l_tx=None, l_rx=None, m=None)

r-u-v measurement Jacobian evaluated at the measurement itself.

Port of calcRuvConvJacob.

Parameters:

• z (array_like) – (3,) measurement [r, u, v].

• use_half_range (bool) – As in calc_ruv_jacob().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob().

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_jacob().

Returns:

J – (3, 3) Jacobian.

Return type:

ndarray

Examples

>>> J = calc_ruv_conv_jacob([2e3, 0.1, 0.2])
>>> J.shape
(3, 3)

pytcl.coordinate_systems.jacobians.calc_ruv_rr_conv_jacob(z, use_half_range=False, l_tx=None, l_rx=None, m=None)

Converted r-u-v Jacobian with a pass-through range-rate row.

Port of calcRuvRRConvJacob.

Parameters:

• z (array_like) – (4,) measurement [r, u, v, range rate].

• use_half_range (bool) – As in calc_ruv_conv_jacob().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_conv_jacob().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_conv_jacob().

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_ruv_conv_jacob().

Returns:

J – (4, 4) block-diagonal Jacobian.

Return type:

ndarray

Examples

>>> J = calc_ruv_rr_conv_jacob([2e3, 0.1, 0.2, 25.0])
>>> float(J[3, 3])
1.0

pytcl.coordinate_systems.jacobians.norm_vec_jacob(x)

Jacobian of the unit vector x/norm(x) with respect to x.

Port of normVecJacob.

Parameters:

x (array_like) – (d,) vector.

Returns:

J – (d, d) Jacobian.

Return type:

ndarray

Examples

>>> norm_vec_jacob([2.0, 0.0])
array([[0. , 0. ],
       [0. , 0.5]])

pytcl.coordinate_systems.jacobians.enu_jacobian(lat, lon)

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)

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)

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)

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
>>> from pytcl.coordinate_systems.jacobians import calc_spher_inv_jacob
>>> P_sph = np.diag([1, 0.01, 0.01])  # [r, az, el] variances
>>> r, az, el = 1000, np.radians(45), np.radians(30)
>>> J = calc_spher_inv_jacob([r, az, el])
>>> P_cart = cross_covariance_transform(J, P_sph)

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

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

Hessians

Second derivatives of measurement components (bistatic range, spherical angles, direction cosines), the composed spherical measurement Hessians, chain-rule helpers and the cross gradients/Hessians between angular parameterizations.

Hessians of measurement functions.

Ports of the MATLAB TCL’s Coordinate_Systems/Hessians: second derivatives of bistatic range, spherical angles and direction cosines with respect to Cartesian position, the composed spherical measurement Hessians (direct, inverse and converted), and the chain-rule/affine helpers. Layout mirrors MATLAB: a Hessian stack has shape (n, n, num_out) with H[:, :, k] the Hessian of the k-th output component.

pytcl.coordinate_systems.hessians.range_hessian(x, use_half_range=False, l_tx=None, l_rx=None)

Hessian of bistatic range with respect to Cartesian position.

A transmitter collocated with the target contributes nothing, as in the gradient. Port of rangeHessian.

Parameters:

• x (array_like) – (d,) target position, d in {1, 2, 3}.

• use_half_range (bool, optional) – Divide by two for one-way ranges. Default False.

• l_tx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

• l_rx (array_like, optional) – Transmitter/receiver positions. Default: the origin.

Returns:

H – (d, d) Hessian.

Return type:

ndarray

Examples

>>> H = range_hessian([3.0, 4.0])
>>> H.shape
(2, 2)

pytcl.coordinate_systems.hessians.spher_ang_hessian(x, system_type=0, l_rx=None, m=None)

Hessians of spherical azimuth and elevation w.r.t. position.

Port of spherAngHessian.

Parameters:

• x (array_like) – (3,) global Cartesian position.

• system_type (int, optional) – Angle convention 0-3. Default 0.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

Returns:

H – (3, 3, 2) stack: H[:, :, 0] the azimuth Hessian, H[:, :, 1] the elevation Hessian.

Return type:

ndarray

Examples

>>> H = spher_ang_hessian([1e3, 2e3, 500.0])
>>> H.shape
(3, 3, 2)

pytcl.coordinate_systems.hessians.uv_hessian(x, l_rx=None, m=None, include_w=False)

Hessians of the u-v(-w) direction cosines w.r.t. position.

Port of uvHessian.

Parameters:

• x (array_like) – (3,) global Cartesian position.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

• include_w (bool, optional) – Also return the third direction cosine’s Hessian. Default False.

Returns:

H – (3, 3, 2) stack ((3, 3, 3) with include_w).

Return type:

ndarray

Examples

>>> H = uv_hessian([1e3, 2e3, 5e3])
>>> H.shape
(3, 3, 2)

pytcl.coordinate_systems.hessians.u_hessian_2d(x, l_rx=None, m=None, include_v=False)

Hessian of the 2-D direction cosine u (optionally v) w.r.t. position.

Port of uHessian2D.

Parameters:

• x (array_like) – (2,) position.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (2, 2) global-to-local rotation. Default: identity.

• include_v (bool, optional) – Also return the second direction cosine’s Hessian. Default False.

Returns:

H – (2, 2, 1) stack ((2, 2, 2) with include_v).

Return type:

ndarray

Examples

>>> H = u_hessian_2d([3.0, 4.0])
>>> H.shape
(2, 2, 1)

pytcl.coordinate_systems.hessians.u_hessian_3d(x, l_rx=None, m=None)

Hessian of the first 3-D direction cosine u w.r.t. position.

Port of uHessian3D.

Parameters:

• x (array_like) – (3,) position.

• l_rx (array_like, optional) – Receiver position. Default: the origin.

• m (array_like, optional) – (3, 3) global-to-local rotation. Default: identity.

Returns:

H – (3, 3, 1) stack.

Return type:

ndarray

Examples

>>> H = u_hessian_3d([1e3, 2e3, 5e3])
>>> H.shape
(3, 3, 1)

pytcl.coordinate_systems.hessians.calc_spher_hessian(x, system_type=0, use_half_range=None, l_tx=None, l_rx=None, m=None)

Hessians of a bistatic [range; azimuth; elevation] measurement.

Port of calcSpherHessian.

Parameters:

• x (array_like) – (3,) target position.

• system_type (int) – As in calc_spher_jacob() (the same MATLAB half-range default applies).

• use_half_range (bool | None) – As in calc_spher_jacob() (the same MATLAB half-range default applies).

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob() (the same MATLAB half-range default applies).

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob() (the same MATLAB half-range default applies).

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_jacob() (the same MATLAB half-range default applies).

Returns:

H – (3, 3, 3) stack: range, azimuth and elevation Hessians.

Return type:

ndarray

Examples

>>> H = calc_spher_hessian([1e3, 2e3, 500.0])
>>> H.shape
(3, 3, 3)

pytcl.coordinate_systems.hessians.calc_spher_inv_hessian(z, system_type=0)

Hessians of Cartesian position w.r.t. a spherical triple.

Port of calcSpherInvHessian.

Parameters:

• z (array_like) – (3,) measurement [r, az, el].

• system_type (int, optional) – Spherical angle convention 0-3. Default 0.

Returns:

H – (3, 3, 3) stack: the Hessians of x, y and z with respect to (r, az, el).

Return type:

ndarray

Examples

>>> H = calc_spher_inv_hessian([1e3, 0.5, 0.2])
>>> H.shape
(3, 3, 3)

pytcl.coordinate_systems.hessians.calc_spher_conv_hessian(z, system_type=0, use_half_range=None, l_tx=None, l_rx=None, m=None)

Spherical measurement Hessian evaluated at the measurement itself.

Port of calcSpherConvHessian.

Parameters:

• z (array_like) – (3,) measurement [r, az, el].

• system_type (int) – As in calc_spher_hessian().

• use_half_range (bool | None) – As in calc_spher_hessian().

• l_tx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_hessian().

• l_rx (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_hessian().

• m (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in calc_spher_hessian().

Returns:

H – (3, 3, 3) stack.

Return type:

ndarray

Examples

>>> H = calc_spher_conv_hessian([1e3, 0.5, 0.2])
>>> H.shape
(3, 3, 3)

pytcl.coordinate_systems.hessians.hessian_of_affine_trans_fun(h, m)

Hessian stack of M @ f(x) given the Hessian stack of f(x).

Port of HessianOfAffineTransFun.

Parameters:

• h (array_like) – (n, n, num_out) Hessian stack of f.

• m (array_like) – (num_out_new, num_out) linear map applied to f’s outputs.

Returns:

H – (n, n, num_out_new) Hessian stack of M @ f.

Return type:

ndarray

Notes

The MATLAB original assigns M*H back into H in place and so accepts only square maps; rectangular output maps are a deliberate widening here.

Examples

>>> import numpy as np
>>> H = np.zeros((2, 2, 2)); H[:, :, 0] = np.eye(2)
>>> hessian_of_affine_trans_fun(H, [[0.0, 1.0], [1.0, 0.0]]).shape
(2, 2, 2)

pytcl.coordinate_systems.hessians.hessian_chain_rule(h_f, h_g, j_f, j_g)

Hessian stack of the composition f(g(x)).

Given f’s Hessians/Jacobian in its own argument and g’s Hessians/Jacobian in x, returns the Hessians of f(g(x)) in x.

Port of HessianChainRule.

Parameters:

• h_f (array_like) – (num_in, num_in, num_out) Hessian stack of f at g(x).

• h_g (array_like) – (num_x, num_x, num_in) Hessian stack of g at x.

• j_f (array_like) – (num_out, num_in) Jacobian of f at g(x).

• j_g (array_like) – (num_in, num_x) Jacobian of g at x.

Returns:

H – (num_x, num_x, num_out) Hessian stack of the composition.

Return type:

ndarray

Examples

>>> import numpy as np
>>> hf = np.zeros((2, 2, 2)); hg = np.zeros((2, 2, 2))
>>> hessian_chain_rule(hf, hg, np.eye(2), 2 * np.eye(2)).shape
(2, 2, 2)

pytcl.coordinate_systems.hessians.uv_spher_ang_cross_grad(az_el, system_type=0, include_w=False, ms=None, muv=None)

Gradient of u-v(-w) direction cosines w.r.t. spherical angles.

The spherical angles are defined in the frame rotated by ms, the direction cosines in the frame rotated by muv.

Port of uvSpherAngCrossGrad.

Parameters:

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

• system_type (int, optional) – Spherical angle convention 0-3. Default 0.

• include_w (bool, optional) – Also return the third direction cosine’s row. Default False.

• ms (array_like, optional) – (3, 3) global-to-local rotations of the spherical and u-v sensor frames. Default: identity.

• muv (array_like, optional) – (3, 3) global-to-local rotations of the spherical and u-v sensor frames. Default: identity.

Returns:

J – (2, 2) d[u; v]/d[az, el] ((3, 2) with include_w).

Return type:

ndarray

Examples

>>> J = uv_spher_ang_cross_grad([0.5, 0.2])
>>> J.shape
(2, 2)

pytcl.coordinate_systems.hessians.spher_ang_uv_cross_grad(uv, system_type=0, ms=None, muv=None)

Gradient of spherical angles w.r.t. u-v direction cosines.

Port of spherAngUvCrossGrad.

Parameters:

• uv (array_like) – (2,) [u, v] (w inferred on the +w hemisphere) or (3,) [u, v, w].

• system_type (int, optional) – Spherical angle convention 0-3. Default 0.

• ms (array_like, optional) – (3, 3) frame rotations as in uv_spher_ang_cross_grad().

• muv (array_like, optional) – (3, 3) frame rotations as in uv_spher_ang_cross_grad().

Returns:

J – (2, 2) d[az; el]/d[u, v].

Return type:

ndarray

Examples

>>> J = spher_ang_uv_cross_grad([0.1, 0.2])
>>> J.shape
(2, 2)

pytcl.coordinate_systems.hessians.u_polar_2d_cross_grad(azimuth, system_type=0, mp=None, mu=None, include_v=False)

Gradient of the 2-D direction cosine(s) w.r.t. the polar angle.

Port of uPolar2DCrossGrad (whose rotation arguments are unusable upstream; see the module docstring).

Parameters:

• azimuth (array_like) – Scalar or (N,) polar angles in radians.

• system_type (int, optional) – Polar convention (0 or 1). Default 0.

• mp (array_like, optional) – (2, 2) rotations of the polar and direction-cosine frames. Default: identity.

• mu (array_like, optional) – (2, 2) rotations of the polar and direction-cosine frames. Default: identity.

• include_v (bool, optional) – Also return the second direction cosine’s row. Default False.

Returns:

J – (1, N) ((2, N) with include_v).

Return type:

ndarray

Examples

>>> u_polar_2d_cross_grad(0.0)
array([[-0.]])

pytcl.coordinate_systems.hessians.polar_u_2d_cross_grad(u_list, system_type=0)

Gradient of the polar angle w.r.t. the 2-D direction cosine(s).

Port of polarU2DCrossGrad.

Parameters:

• u_list (array_like) – (1, N) or (N,) u values (v inferred on the +v half-plane), or (2, N) [u; v] columns.

• system_type (int, optional) – Polar convention (0 or 1). Default 0.

Returns:

J – (1, N) d(az)/du, or (2, N) [d(az)/du; d(az)/dv] when v was given.

Return type:

ndarray

Examples

>>> polar_u_2d_cross_grad(0.0)
array([[-1.]])

pytcl.coordinate_systems.hessians.uv_spher_ang_cross_hessian(az_el, system_type=0, include_w=False, ms=None, muv=None)

Hessians of u-v(-w) direction cosines w.r.t. spherical angles.

Port of uvSpherAngCrossHessian.

Parameters:

• az_el (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – As in uv_spher_ang_cross_grad().

• system_type (int) – As in uv_spher_ang_cross_grad().

• include_w (bool) – As in uv_spher_ang_cross_grad().

• ms (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in uv_spher_ang_cross_grad().

• muv (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in uv_spher_ang_cross_grad().

Returns:

H – (2, 2, 2) stack over (az, el) for u and v ((2, 2, 3) with include_w).

Return type:

ndarray

Examples

>>> H = uv_spher_ang_cross_hessian([0.5, 0.2])
>>> H.shape
(2, 2, 2)

pytcl.coordinate_systems.hessians.spher_ang_uv_cross_hessian(uv, system_type=0, ms=None, muv=None)

Hessians of spherical angles w.r.t. u-v direction cosines.

Port of spherAngUvCrossHessian.

Parameters:

• uv (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – As in spher_ang_uv_cross_grad().

• system_type (int) – As in spher_ang_uv_cross_grad().

• ms (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in spher_ang_uv_cross_grad().

• muv (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str] | None) – As in spher_ang_uv_cross_grad().

Returns:

H – (2, 2, 2) stack over (u, v) for azimuth and elevation.

Return type:

ndarray

Examples

>>> H = spher_ang_uv_cross_hessian([0.1, 0.2])
>>> H.shape
(2, 2, 2)

pytcl.coordinate_systems.hessians.u_polar_2d_cross_hessian(azimuth, system_type=0, include_v=False)

Second derivative of the 2-D direction cosine(s) w.r.t. the angle.

Port of uPolar2DCrossHessian.

Parameters:

• azimuth (array_like) – Scalar or (N,) polar angles in radians.

• system_type (int, optional) – Polar convention (0 or 1). Default 0.

• include_v (bool, optional) – Also return the second direction cosine’s row. Default False.

Returns:

H – (1, N) ((2, N) with include_v).

Return type:

ndarray

Examples

>>> u_polar_2d_cross_hessian(0.0)
array([[-1.]])

pytcl.coordinate_systems.hessians.polar_u_2d_cross_hessian(u_list, system_type=0)

Second derivatives of the polar angle w.r.t. the direction cosines.

Port of polarU2DCrossHessian.

Parameters:

• u_list (array_like) – (1, N) or (N,) u values (v inferred on the +v half-plane), or (2, N) [u; v] columns.

• system_type (int, optional) – Polar convention (0 or 1). Default 0.

Returns:

H – (1, 1, N), or (2, 2, N) when v was given, matching the MATLAB layout.

Return type:

ndarray

Examples

>>> polar_u_2d_cross_hessian(0.0)[..., 0]
array([[-0.]])

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)

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)

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))

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)

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))

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))

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)

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)

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)

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')

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)

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))

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

Do not use lat0 = +-pi/2 for polar work: this function’s conformal latitude substitution diverges from PROJ’s polar stereographic by 5.7 km at 85 degrees latitude and 34.7 km at 60 degrees (measured against +proj=stere on WGS84). Use polar_stereographic() instead, which matches PROJ/UPS to sub-nanometer precision. For the ellipsoidal oblique case (any other lat0), use oblique_stereographic() (EPSG method 9809), which matches PROJ’s +proj=sterea to sub-micrometer precision; this function’s own oblique accuracy is documented below and is generally worse.

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 coordinates that match PROJ (to sub-micrometer), use oblique_stereographic() – the EPSG method 9809 implementation.

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)

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.oblique_stereographic(lat, lon, lat0, lon0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))

Oblique (double) stereographic projection, EPSG method 9809 (forward).

The ellipsoid is first mapped conformally onto a sphere of Gaussian radius at the origin using the full EPSG latitude mapping (constants n and c), then that sphere is projected stereographically. This is the algorithm PROJ implements as +proj=sterea (the projection of the Dutch RD grid, EPSG:28992), so unlike stereographic() – whose simpler conformal-latitude substitution diverges from sterea by kilometers far from the origin (gh-25) – coordinates from this function are interchangeable with PROJ output.

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• lat0 (float) – Origin latitude in radians.

• lon0 (float) – Origin longitude in radians.

• k0 (float, optional) – Scale factor at the origin. 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

References

• IOGP Publication 373-7-2, “Coordinate Conversions and Transformations including Formulas”, EPSG method code 9809.

Examples

>>> import numpy as np
>>> res = oblique_stereographic(np.radians(53.0), np.radians(6.0),
...                             np.radians(52.156161), np.radians(5.387639),
...                             k0=0.9999079)
>>> bool(abs(res.scale - 0.9999079) < 1e-3)
True

pytcl.coordinate_systems.projections.projections.oblique_stereographic_inverse(x, y, lat0, lon0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))

Oblique (double) stereographic projection, EPSG method 9809 (inverse).

Parameters:

• x (float) – Easting in meters.

• y (float) – Northing in meters.

• lat0 (float) – Origin latitude in radians.

• lon0 (float) – Origin longitude in radians.

• k0 (float, optional) – Scale factor at the origin. 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]

Examples

>>> import numpy as np
>>> lat0, lon0 = np.radians(52.156161), np.radians(5.387639)
>>> fwd = oblique_stereographic(np.radians(53.0), np.radians(6.0),
...                             lat0, lon0, k0=0.9999079)
>>> lat, lon = oblique_stereographic_inverse(fwd.x, fwd.y, lat0, lon0,
...                                          k0=0.9999079)
>>> bool(abs(np.degrees(lat) - 53.0) < 1e-9)
True

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

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))

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)

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))

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 exact ellipsoidal distances (matching PROJ’s +proj=aeqd), use azimuthal_equidistant_exact().

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))

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]

pytcl.coordinate_systems.projections.projections.azimuthal_equidistant_exact(lat, lon, lat0, lon0, a=6378137.0, e2=np.float64(0.006694379990141316))

Exact ellipsoidal azimuthal equidistant projection (forward).

The map coordinates are the true geodesic distance and azimuth from the projection centre: x = s12 sin(az1), y = s12 cos(az1), computed with Karney’s algorithm (geographiclib). This is what PROJ’s +proj=aeqd computes on an ellipsoid, so unlike azimuthal_equidistant() – a spherical approximation on the authalic sphere – distances from the centre are exact everywhere, antipodes included.

Requires the geodesy extra (geographiclib).

Parameters:

• lat (float) – Geodetic latitude in radians.

• lon (float) – Longitude in radians.

• lat0 (float) – Centre latitude in radians.

• lon0 (float) – Centre 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. scale is the radial scale, exactly 1 by construction; convergence is the true grid convergence, computed from the geodesic azimuths and reduced length.

Return type:

ProjectionResult

Raises:

DependencyError – If geographiclib is not installed.

Examples

>>> import numpy as np
>>> res = azimuthal_equidistant_exact(np.radians(40), np.radians(-75),
...                                   np.radians(38.9), np.radians(-77))
>>> bool(np.hypot(res.x, res.y) > 0)
True
>>> res.scale
1.0

pytcl.coordinate_systems.projections.projections.azimuthal_equidistant_exact_inverse(x, y, lat0, lon0, a=6378137.0, e2=np.float64(0.006694379990141316))

Exact ellipsoidal azimuthal equidistant projection (inverse).

Solves the direct geodesic problem from the projection centre with azimuth atan2(x, y) and distance hypot(x, y) (Karney’s algorithm, geographiclib). Requires the geodesy extra.

Parameters:

• x (float) – Easting in meters.

• y (float) – Northing in meters.

• lat0 (float) – Centre latitude in radians.

• lon0 (float) – Centre 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]

Raises:

DependencyError – If geographiclib is not installed.

Examples

>>> import numpy as np
>>> lat0, lon0 = np.radians(38.9), np.radians(-77)
>>> fwd = azimuthal_equidistant_exact(np.radians(40), np.radians(-75),
...                                   lat0, lon0)
>>> lat, lon = azimuthal_equidistant_exact_inverse(fwd.x, fwd.y,
...                                                lat0, lon0)
>>> bool(abs(np.degrees(lat) - 40.0) < 1e-9)
True