Coordinate Systems

Coordinate system conversions and transformations.

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

  • Cartesian coordinates (x, y, z)

  • Spherical coordinates (range, azimuth, elevation)

  • Polar and cylindrical coordinates

  • Geodetic coordinates (latitude, longitude, altitude)

  • Various local tangent plane frames (ENU, NED)

  • Direction cosine representations (r-u-v)

  • Rotation representations (matrices, quaternions, Euler angles)

  • Jacobian matrices for error propagation

Conversions

Coordinate conversions.

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

Spherical Coordinates

Spherical and polar coordinate conversions.

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

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

Convert Cartesian coordinates to spherical coordinates.

Parameters:
  • cart_points (array_like) – Cartesian coordinates. Can be: - Shape (3,) for a single point [x, y, z] - Shape (3, n) for n points (each column is a point) - Shape (n, 3) will be transposed automatically

  • system_type ({'standard', 'az-el', 'range-az-el'}, optional) – Spherical coordinate convention: - ‘standard’: Physics convention (r, θ polar from +z, φ azimuth from +x) - ‘az-el’: Tracking convention (r, azimuth from +x, elevation from xy-plane) - ‘range-az-el’: Same as ‘az-el’ (alias) Default is ‘standard’.

Returns:

  • r (ndarray) – Range (radial distance from origin).

  • az (ndarray) – Azimuth angle in radians. - ‘standard’: Angle in xy-plane from +x axis [0, 2π) - ‘az-el’: Angle in xy-plane from +x axis [-π, π]

  • el (ndarray) – Elevation/polar angle in radians. - ‘standard’: Polar angle from +z axis [0, π] - ‘az-el’: Elevation from xy-plane [-π/2, π/2]

Return type:

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

Examples

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

See also

sphere2cart

Inverse conversion.

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

Convert spherical coordinates to Cartesian coordinates.

Parameters:
  • r (array_like) – Range (radial distance).

  • az (array_like) – Azimuth angle in radians.

  • el (array_like) – Elevation/polar angle in radians.

  • system_type ({'standard', 'az-el', 'range-az-el'}, optional) – Spherical coordinate convention (see cart2sphere).

Returns:

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

Return type:

ndarray

Examples

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

See also

cart2sphere

Inverse conversion.

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

Convert 2D Cartesian coordinates to polar coordinates.

Parameters:

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

Returns:

  • r (ndarray) – Radial distance from origin.

  • theta (ndarray) – Angle in radians from +x axis, in range [-π, π].

Return type:

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

Examples

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

See also

pol2cart

Inverse conversion.

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

Convert polar coordinates to 2D Cartesian coordinates.

Parameters:
  • r (array_like) – Radial distance.

  • theta (array_like) – Angle in radians from +x axis.

Returns:

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

Return type:

ndarray

Examples

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

See also

cart2pol

Inverse conversion.

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

Convert 3D Cartesian coordinates to cylindrical coordinates.

Parameters:

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

Returns:

  • rho (ndarray) – Radial distance in xy-plane.

  • phi (ndarray) – Azimuth angle in radians from +x axis.

  • z (ndarray) – Height (same as Cartesian z).

Return type:

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

Examples

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

See also

cyl2cart

Inverse conversion.

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

Convert cylindrical coordinates to 3D Cartesian coordinates.

Parameters:
  • rho (array_like) – Radial distance in xy-plane.

  • phi (array_like) – Azimuth angle in radians from +x axis.

  • z (array_like) – Height.

Returns:

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

Return type:

ndarray

Examples

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

See also

cart2cyl

Inverse conversion.

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

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

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

Parameters:
  • r (array_like) – Range.

  • u (array_like) – Direction cosine along x-axis.

  • v (array_like) – Direction cosine along y-axis.

Returns:

cart_points – Cartesian coordinates.

Return type:

ndarray

Examples

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

Notes

This representation is common in radar tracking systems.

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

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

Parameters:

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

Returns:

  • r (ndarray) – Range.

  • u (ndarray) – Direction cosine along x-axis (x/r).

  • v (ndarray) – Direction cosine along y-axis (y/r).

Return type:

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

Examples

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

See also

ruv2cart

Inverse conversion.

Geodetic Coordinates

Geodetic coordinate conversions.

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

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

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

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • lon (array_like) – Geodetic longitude in radians.

  • alt (array_like) – Altitude above the reference ellipsoid in meters.

  • a (float, optional) – Semi-major axis of the reference ellipsoid in meters. Default is WGS84 value.

  • f (float, optional) – Flattening of the reference ellipsoid. Default is WGS84 value.

Returns:

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

Return type:

ndarray

Examples

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

See also

ecef2geodetic

Inverse conversion.

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

Convert ECEF coordinates to geodetic coordinates.

Parameters:
  • ecef (array_like) – ECEF coordinates [x, y, z] in meters.

  • a (float, optional) – Semi-major axis of the reference ellipsoid.

  • f (float, optional) – Flattening of the reference ellipsoid.

  • method (str, optional) – Algorithm to use: - ‘iterative’: Bowring’s iterative method (default) - ‘direct’: Closed-form solution (Vermeille’s method)

Returns:

  • lat (ndarray) – Geodetic latitude in radians.

  • lon (ndarray) – Geodetic longitude in radians.

  • alt (ndarray) – Altitude above the ellipsoid in meters.

Return type:

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

Examples

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

See also

geodetic2ecef

Inverse conversion.

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

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

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • lon (array_like) – Geodetic longitude in radians.

  • alt (array_like) – Altitude in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • alt_ref (float) – Reference point altitude in meters.

  • a (float, optional) – Semi-major axis of the reference ellipsoid.

  • f (float, optional) – Flattening of the reference ellipsoid.

Returns:

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

Return type:

ndarray

Examples

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

See also

enu2geodetic

Inverse conversion.

ecef2enu

ECEF to ENU conversion.

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

Convert ECEF coordinates to local ENU coordinates.

Parameters:
  • ecef (array_like) – ECEF coordinates [x, y, z] in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference point ECEF coordinates. If None, computed from lat_ref, lon_ref.

Returns:

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

Return type:

ndarray

Examples

Convert an aircraft position from ECEF to local ENU frame:

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

See also

enu2ecef

Inverse conversion.

ecef2ned

Convert to NED (North-East-Down) frame.

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

Convert local ENU coordinates to ECEF coordinates.

Parameters:
  • enu (array_like) – Local ENU coordinates [east, north, up] in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

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

Return type:

ndarray

Examples

Convert local ENU offset to ECEF coordinates:

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

See also

ecef2enu

Inverse conversion.

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

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

Parameters:
  • ecef (array_like) – ECEF coordinates [x, y, z] in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

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

Return type:

ndarray

Examples

Convert aircraft ECEF position to local NED frame:

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

See also

ned2ecef

Inverse conversion.

ecef2enu

Similar, but ENU frame.

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

Convert local NED coordinates to ECEF coordinates.

Parameters:
  • ned (array_like) – Local NED coordinates [north, east, down] in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference point ECEF coordinates.

Returns:

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

Return type:

ndarray

Examples

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

See also

ecef2ned

Inverse conversion.

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

Convert ENU coordinates to NED coordinates.

Parameters:

enu (array_like) – ENU coordinates [east, north, up].

Returns:

ned – NED coordinates [north, east, down].

Return type:

ndarray

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

Convert NED coordinates to ENU coordinates.

Parameters:

ned (array_like) – NED coordinates [north, east, down].

Returns:

enu – ENU coordinates [east, north, up].

Return type:

ndarray

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

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

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

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • lon (array_like) – Geodetic longitude in radians.

  • alt (array_like) – Altitude in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • alt_ref (float) – Reference point altitude in meters.

  • a (float, optional) – Semi-major axis of the reference ellipsoid.

  • f (float, optional) – Flattening of the reference ellipsoid.

Returns:

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

Return type:

ndarray

See also

sez2geodetic

Inverse conversion.

ecef2sez

ECEF to SEZ conversion.

Notes

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

Examples

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

Convert ECEF coordinates to local SEZ coordinates.

Parameters:
  • ecef (array_like) – ECEF coordinates [X, Y, Z] in meters, shape (3,) or (3, N).

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference ECEF position. If None, the reference point is at (lat_ref, lon_ref) with zero altitude.

Returns:

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

Return type:

ndarray

See also

sez2ecef

Inverse conversion.

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

Convert local SEZ coordinates to ECEF coordinates.

Parameters:
  • sez (array_like) – SEZ coordinates [south, east, zenith] in meters, shape (3,) or (3, N).

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • ecef_ref (array_like, optional) – Reference ECEF position. If None, the reference point is at (lat_ref, lon_ref) with zero altitude.

Returns:

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

Return type:

ndarray

See also

ecef2sez

Forward conversion.

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

Convert local SEZ coordinates to geodetic coordinates.

Parameters:
  • sez (array_like) – SEZ coordinates [south, east, zenith] in meters.

  • lat_ref (float) – Reference point latitude in radians.

  • lon_ref (float) – Reference point longitude in radians.

  • alt_ref (float) – Reference point altitude in meters.

  • a (float, optional) – Semi-major axis.

  • f (float, optional) – Flattening.

Returns:

  • lat (ndarray) – Geodetic latitude in radians.

  • lon (ndarray) – Geodetic longitude in radians.

  • alt (ndarray) – Altitude in meters.

Return type:

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

Examples

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

See also

geodetic2sez

Forward conversion.

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

Compute the geocentric radius at a given geodetic latitude.

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • a (float, optional) – Semi-major axis.

  • f (float, optional) – Flattening.

Returns:

r – Geocentric radius in meters.

Return type:

ndarray

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

Compute the prime vertical radius of curvature.

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • a (float, optional) – Semi-major axis.

  • f (float, optional) – Flattening.

Returns:

N – Prime vertical radius of curvature in meters.

Return type:

ndarray

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

Compute the meridional radius of curvature.

Parameters:
  • lat (array_like) – Geodetic latitude in radians.

  • a (float, optional) – Semi-major axis.

  • f (float, optional) – Flattening.

Returns:

M – Meridional radius of curvature in meters.

Return type:

ndarray

Rotations

Rotation representations and conversions.

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

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

Create rotation matrix for rotation about the x-axis.

Parameters:

angle (float) – Rotation angle in radians.

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

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

Create rotation matrix for rotation about the y-axis.

Parameters:

angle (float) – Rotation angle in radians.

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

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

Create rotation matrix for rotation about the z-axis.

Parameters:

angle (float) – Rotation angle in radians.

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

>>> R = rotz(np.pi/2)  # 90 degree rotation about z
>>> bool(np.allclose(R @ [1, 0, 0], [0, 1, 0]))  # x-axis maps to y-axis
True
pytcl.coordinate_systems.rotations.euler2rotmat(angles, sequence='ZYX')[source]

Convert Euler angles to rotation matrix.

Parameters:
  • angles (array_like) – Three Euler angles in radians [angle1, angle2, angle3].

  • sequence (str, optional) – Rotation sequence (e.g., ‘ZYX’, ‘XYZ’, ‘ZXZ’). Default is ‘ZYX’ (aerospace convention: yaw-pitch-roll).

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

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

Notes

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

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

Convert rotation matrix to Euler angles.

Parameters:
  • R (array_like) – 3x3 rotation matrix.

  • sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

angles – Three Euler angles in radians.

Return type:

ndarray

Examples

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

Notes

May have singularities (gimbal lock) at certain angles.

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

Convert axis-angle representation to rotation matrix.

Parameters:
  • axis (array_like) – Unit vector defining the rotation axis [ax, ay, az].

  • angle (float) – Rotation angle in radians.

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

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

Notes

Uses Rodrigues’ rotation formula.

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

Convert rotation matrix to axis-angle representation.

Parameters:

R (array_like) – 3x3 rotation matrix.

Returns:

  • axis (ndarray) – Unit vector rotation axis.

  • angle (float) – Rotation angle in radians [0, π].

Return type:

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

Examples

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

Convert quaternion to rotation matrix.

Parameters:

q (array_like) – Quaternion [qw, qx, qy, qz] (scalar-first convention).

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

>>> q = [1, 0, 0, 0]  # Identity quaternion
>>> quat2rotmat(q)
array([[1., 0., 0.],
       [0., 1., 0.],
       [0., 0., 1.]])
pytcl.coordinate_systems.rotations.rotmat2quat(R)[source]

Convert rotation matrix to quaternion.

Parameters:

R (array_like) – 3x3 rotation matrix.

Returns:

q – Quaternion [qw, qx, qy, qz] (scalar-first, positive qw).

Return type:

ndarray

Examples

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

Notes

Uses Shepperd’s method for numerical stability.

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

Convert Euler angles to quaternion.

Parameters:
  • angles (array_like) – Three Euler angles in radians.

  • sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

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

Return type:

ndarray

Examples

Convert yaw-pitch-roll angles to quaternion:

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

See also

quat2euler

Inverse conversion.

euler2rotmat

Convert to rotation matrix instead.

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

Convert quaternion to Euler angles.

Parameters:
  • q (array_like) – Quaternion [qw, qx, qy, qz].

  • sequence (str, optional) – Rotation sequence. Default is ‘ZYX’.

Returns:

angles – Three Euler angles in radians.

Return type:

ndarray

Examples

>>> q = [1, 0, 0, 0]  # Identity quaternion
>>> angles = quat2euler(q, 'ZYX')
>>> np.allclose(angles, [0, 0, 0])
True
pytcl.coordinate_systems.rotations.quat_multiply(q1, q2)[source]

Multiply two quaternions.

Parameters:
  • q1 (array_like) – First quaternion [qw, qx, qy, qz].

  • q2 (array_like) – Second quaternion.

Returns:

q – Product quaternion q1 * q2.

Return type:

ndarray

Notes

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

Examples

Combine two rotations using quaternion multiplication:

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

See also

quat_inverse

Compute quaternion inverse.

quat_rotate

Rotate a vector by a quaternion.

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

Compute quaternion conjugate.

Parameters:

q (array_like) – Quaternion [qw, qx, qy, qz].

Returns:

q_conj – Conjugate quaternion [qw, -qx, -qy, -qz].

Return type:

ndarray

Examples

>>> quat_conjugate([0.707, 0.707, 0, 0])
array([ 0.707, -0.707, -0.   , -0.   ])
pytcl.coordinate_systems.rotations.quat_inverse(q)[source]

Compute quaternion inverse.

Parameters:

q (array_like) – Quaternion [qw, qx, qy, qz].

Returns:

q_inv – Inverse quaternion.

Return type:

ndarray

Examples

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

Notes

For unit quaternions, inverse equals conjugate.

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

Rotate a vector using a quaternion.

Parameters:
  • q (array_like) – Quaternion [qw, qx, qy, qz].

  • v (array_like) – Vector to rotate [x, y, z].

Returns:

v_rot – Rotated vector.

Return type:

ndarray

Examples

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

Notes

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

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

Spherical linear interpolation between two quaternions.

Parameters:
  • q1 (array_like) – Start quaternion.

  • q2 (array_like) – End quaternion.

  • t (float) – Interpolation parameter in [0, 1].

Returns:

q – Interpolated quaternion.

Return type:

ndarray

Examples

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

Convert Rodrigues vector to rotation matrix.

Parameters:

rvec (array_like) – Rodrigues vector (axis * angle).

Returns:

R – 3x3 rotation matrix.

Return type:

ndarray

Examples

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

Notes

The Rodrigues vector encodes both the rotation axis and angle:

rvec = axis * angle, where |rvec| = angle
pytcl.coordinate_systems.rotations.rotmat2rodrigues(R)[source]

Convert rotation matrix to Rodrigues vector.

Parameters:

R (array_like) – 3x3 rotation matrix.

Returns:

rvec – Rodrigues vector (axis * angle).

Return type:

ndarray

Examples

>>> R = rotz(np.pi/2)  # 90 deg about Z
>>> rvec = rotmat2rodrigues(R)
>>> np.linalg.norm(rvec)  # magnitude is the angle
1.5707...
pytcl.coordinate_systems.rotations.dcm_rate(R, omega)[source]

Compute the time derivative of a rotation matrix.

Parameters:
  • R (array_like) – Current rotation matrix.

  • omega (array_like) – Angular velocity vector [wx, wy, wz] in body frame.

Returns:

R_dot – Time derivative of R.

Return type:

ndarray

Examples

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

Notes

R_dot = R @ skew(omega)

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

Check if a matrix is a valid rotation matrix.

Parameters:
  • R (array_like) – Matrix to check.

  • tol (float, optional) – Tolerance for numerical checks.

Returns:

valid – True if R is a valid rotation matrix.

Return type:

bool

Examples

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

Jacobians

Jacobian matrices for coordinate transformations.

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

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

Compute Jacobian of Cartesian to spherical transformation.

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

Parameters:
  • cart_point (array_like) – Cartesian coordinates [x, y, z].

  • system_type ({'standard', 'az-el'}, optional) – Spherical coordinate convention. Default is ‘az-el’.

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

Notes

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

Examples

>>> J = spherical_jacobian([1, 1, 1])
>>> J.shape
(3, 3)
pytcl.coordinate_systems.jacobians.spherical_jacobian_inv(r, az, el, system_type='az-el')[source]

Compute Jacobian of spherical to Cartesian transformation.

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

Parameters:
  • r (float) – Range.

  • az (float) – Azimuth in radians.

  • el (float) – Elevation in radians.

  • system_type ({'standard', 'az-el'}, optional) – Spherical coordinate convention.

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

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

Compute Jacobian of 2D Cartesian to polar transformation.

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

Parameters:

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

Returns:

J – 2x2 Jacobian matrix.

Return type:

ndarray

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

Compute Jacobian of polar to 2D Cartesian transformation.

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

Parameters:
  • r (float) – Radial distance.

  • theta (float) – Angle in radians.

Returns:

J – 2x2 Jacobian matrix.

Return type:

ndarray

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

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

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

Parameters:

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

Returns:

J – 3x3 Jacobian matrix where d[r, u, v] = J @ d[x, y, z].

Return type:

ndarray

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

Compute Jacobian of ECEF to ENU transformation.

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

Parameters:
  • lat (float) – Reference latitude in radians.

  • lon (float) – Reference longitude in radians.

Returns:

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

Return type:

ndarray

Notes

Uses cached computation with quantized inputs for performance.

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

Compute Jacobian of ECEF to NED transformation.

Parameters:
  • lat (float) – Reference latitude in radians.

  • lon (float) – Reference longitude in radians.

Returns:

J – 3x3 rotation matrix.

Return type:

ndarray

Notes

Uses cached computation with quantized inputs for performance.

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

Compute Jacobian of geodetic to ECEF transformation.

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Geodetic longitude in radians.

  • alt (float) – Altitude above ellipsoid in meters.

  • a (float, optional) – Semi-major axis (default: WGS84).

  • f (float, optional) – Flattening (default: WGS84).

Returns:

J – 3x3 Jacobian matrix.

Return type:

ndarray

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

Transform a covariance matrix through a Jacobian.

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

Parameters:
  • J (array_like) – Jacobian matrix of the transformation.

  • P (array_like) – Original covariance matrix.

Returns:

P_new – Transformed covariance matrix.

Return type:

ndarray

Examples

>>> # Transform spherical covariance to Cartesian
>>> P_sph = np.diag([1, 0.01, 0.01])  # [r, az, el] variances
>>> r, az, el = 1000, np.radians(45), np.radians(30)
>>> J = spherical_jacobian_inv(r, az, el)
>>> P_cart = cross_covariance_transform(J, P_sph)
pytcl.coordinate_systems.jacobians.numerical_jacobian(func, x, dx=1e-07)[source]

Compute Jacobian numerically using central differences.

Parameters:
  • func (callable) – Function f(x) -> y.

  • x (array_like) – Point at which to compute Jacobian.

  • dx (float, optional) – Step size for finite differences.

Returns:

J – Jacobian matrix.

Return type:

ndarray

Projections

Map Projections for Tracking and Navigation.

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

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

All angles are in radians unless otherwise noted.

References

  • Snyder, J. P. “Map Projections: A Working Manual.” U.S. Geological Survey Professional Paper 1395, 1987.

  • NIMA Technical Report 8350.2, “Department of Defense World Geodetic System 1984,” Third Edition, 2000.

  • Karney, C. F. F. “Transverse Mercator with an accuracy of a few nanometers.” Journal of Geodesy 85.8 (2011): 475-485.

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

Bases: NamedTuple

Result of a map projection.

Parameters:
  • x (float) – Easting coordinate in meters.

  • y (float) – Northing coordinate in meters.

  • scale (float) – Scale factor at the point (ratio of map distance to ground distance).

  • convergence (float) – Grid convergence angle in radians (angle from true north to grid north).

x: float

Alias for field number 0

y: float

Alias for field number 1

scale: float

Alias for field number 2

convergence: float

Alias for field number 3

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

Bases: NamedTuple

Result of UTM projection.

Parameters:
  • easting (float) – UTM easting coordinate in meters.

  • northing (float) – UTM northing coordinate in meters.

  • zone (int) – UTM zone number (1-60).

  • hemisphere (str) – ‘N’ for northern hemisphere, ‘S’ for southern.

  • scale (float) – Scale factor at the point.

  • convergence (float) – Grid convergence angle in radians.

easting: float

Alias for field number 0

northing: float

Alias for field number 1

zone: int

Alias for field number 2

hemisphere: str

Alias for field number 3

scale: float

Alias for field number 4

convergence: float

Alias for field number 5

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

Ellipsoidal Mercator projection (forward).

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • lon0 (float, optional) – Central meridian in radians. Default is 0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e (float, optional) – First eccentricity. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

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

Examples

>>> import numpy as np
>>> result = mercator(np.radians(45), np.radians(-75))
>>> print(f"x={result.x:.1f}, y={result.y:.1f}")
x=-8348961.8, y=5591295.9
pytcl.coordinate_systems.projections.projections.mercator_inverse(x, y, lon0=0.0, a=6378137.0, e=np.float64(0.08181919084262149), tol=1e-12, max_iter=10)[source]

Ellipsoidal Mercator projection (inverse).

Parameters:
  • x (float) – Easting coordinate in meters.

  • y (float) – Northing coordinate in meters.

  • lon0 (float, optional) – Central meridian in radians. Default is 0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e (float, optional) – First eccentricity. Default is WGS84.

  • tol (float, optional) – Convergence tolerance. Default is 1e-12.

  • max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

Examples

>>> import numpy as np
>>> lat, lon = mercator_inverse(1000000, 5000000)
pytcl.coordinate_systems.projections.projections.transverse_mercator(lat, lon, lat0=0.0, lon0=0.0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Transverse Mercator projection (forward).

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • lat0 (float, optional) – Origin latitude in radians. Default is 0.

  • lon0 (float, optional) – Central meridian in radians. Default is 0.

  • k0 (float, optional) – Scale factor at central meridian. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

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

From the meridian

Difference

2 degrees

< 0.1 mm

3 degrees

0.1 mm

4 degrees

0.7 mm

5 degrees

3.3 mm

6 degrees

11.5 mm

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

Examples

>>> import numpy as np
>>> result = transverse_mercator(np.radians(45), np.radians(-75),
...                              lon0=np.radians(-75))
pytcl.coordinate_systems.projections.projections.transverse_mercator_inverse(x, y, lat0=0.0, lon0=0.0, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Transverse Mercator projection (inverse).

Parameters:
  • x (float) – Easting coordinate in meters.

  • y (float) – Northing coordinate in meters.

  • lat0 (float, optional) – Origin latitude in radians. Default is 0.

  • lon0 (float, optional) – Central meridian in radians. Default is 0.

  • k0 (float, optional) – Scale factor at central meridian. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

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

Determine UTM zone number from longitude.

Parameters:
  • lon (float) – Longitude in radians.

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

Returns:

UTM zone number (1-60).

Return type:

int

Notes

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

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

Get central meridian for UTM zone.

Parameters:

zone (int) – UTM zone number (1-60).

Returns:

Central meridian in radians.

Return type:

float

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

Convert geodetic coordinates to UTM.

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • zone (int, optional) – Force specific UTM zone. If None, computed from longitude.

Returns:

UTM coordinates with zone information.

Return type:

UTMResult

Examples

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

Convert UTM coordinates to geodetic.

Parameters:
  • easting (float) – UTM easting in meters.

  • northing (float) – UTM northing in meters.

  • zone (int) – UTM zone number (1-60).

  • hemisphere (str, optional) – ‘N’ for northern, ‘S’ for southern hemisphere.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

Examples

>>> lat, lon = utm2geodetic(500000, 5000000, 18, 'N')
>>> print(f"Lat: {np.degrees(lat):.4f}, Lon: {np.degrees(lon):.4f}")
Lat: 45.1535, Lon: -75.0000
pytcl.coordinate_systems.projections.projections.geodetic2utm_batch(lats, lons, zone=None)[source]

Batch convert geodetic coordinates to UTM.

Parameters:
  • lats (ndarray) – Geodetic latitudes in radians.

  • lons (ndarray) – Longitudes in radians.

  • zone (int, optional) – Force specific UTM zone. If None, computed per point.

Returns:

  • eastings (ndarray) – UTM eastings in meters.

  • northings (ndarray) – UTM northings in meters.

  • zones (ndarray) – UTM zone numbers.

  • hemispheres (ndarray) – Hemisphere characters (‘N’ or ‘S’).

Return type:

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

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

Oblique stereographic projection (forward).

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • lat0 (float) – Center latitude in radians.

  • lon0 (float) – Center longitude in radians.

  • k0 (float, optional) – Scale factor at center. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

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

Examples

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

Notes

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

lat0

At 400 km

At 1,738 km

0 deg

1.34 km

5.9 km

30 deg

0.03 km

0.73 km

45 deg

1.38 km

6.7 km

60 deg

2.73 km

12.5 km

75 deg

3.71 km

16.5 km

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

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

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

Oblique stereographic projection (inverse).

Parameters:
  • x (float) – Easting in meters.

  • y (float) – Northing in meters.

  • lat0 (float) – Center latitude in radians.

  • lon0 (float) – Center longitude in radians.

  • k0 (float, optional) – Scale factor at center. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

  • tol (float, optional) – Convergence tolerance. Default is 1e-12.

  • max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

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

Polar stereographic projection (forward).

Standard polar stereographic used for polar regions.

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • north (bool, optional) – True for North Pole center, False for South Pole. Default is True.

  • k0 (float, optional) – Scale factor at pole. Default is 0.994 (UPS standard).

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

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

Examples

>>> import numpy as np
>>> # Arctic location
>>> result = polar_stereographic(np.radians(80), np.radians(45))
pytcl.coordinate_systems.projections.projections.lambert_conformal_conic(lat, lon, lat0, lon0, lat1, lat2, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Lambert Conformal Conic projection (forward).

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • lat0 (float) – Origin latitude in radians.

  • lon0 (float) – Central meridian in radians.

  • lat1 (float) – First standard parallel in radians.

  • lat2 (float) – Second standard parallel in radians.

  • k0 (float, optional) – Scale factor. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Examples

>>> import numpy as np
>>> # Continental US projection
>>> result = lambert_conformal_conic(
...     np.radians(40), np.radians(-100),
...     lat0=np.radians(39), lon0=np.radians(-96),
...     lat1=np.radians(33), lat2=np.radians(45)
... )
pytcl.coordinate_systems.projections.projections.lambert_conformal_conic_inverse(x, y, lat0, lon0, lat1, lat2, k0=1.0, a=6378137.0, e2=np.float64(0.006694379990141316), tol=1e-12, max_iter=10)[source]

Lambert Conformal Conic projection (inverse).

Parameters:
  • x (float) – Easting in meters.

  • y (float) – Northing in meters.

  • lat0 (float) – Origin latitude in radians.

  • lon0 (float) – Central meridian in radians.

  • lat1 (float) – First standard parallel in radians.

  • lat2 (float) – Second standard parallel in radians.

  • k0 (float, optional) – Scale factor. Default is 1.0.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

  • tol (float, optional) – Convergence tolerance. Default is 1e-12.

  • max_iter (int, optional) – Maximum iterations. Default is 10.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]

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

Azimuthal equidistant projection (forward).

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

Parameters:
  • lat (float) – Geodetic latitude in radians.

  • lon (float) – Longitude in radians.

  • lat0 (float) – Center latitude in radians.

  • lon0 (float) – Center longitude in radians.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

Projected coordinates with scale and convergence.

Return type:

ProjectionResult

Notes

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

Examples

>>> import numpy as np
>>> result = azimuthal_equidistant(np.radians(40), np.radians(-75),
...                                np.radians(38.9), np.radians(-77))
pytcl.coordinate_systems.projections.projections.azimuthal_equidistant_inverse(x, y, lat0, lon0, a=6378137.0, e2=np.float64(0.006694379990141316))[source]

Azimuthal equidistant projection (inverse).

Parameters:
  • x (float) – Easting in meters.

  • y (float) – Northing in meters.

  • lat0 (float) – Center latitude in radians.

  • lon0 (float) – Center longitude in radians.

  • a (float, optional) – Semi-major axis in meters. Default is WGS84.

  • e2 (float, optional) – First eccentricity squared. Default is WGS84.

Returns:

(latitude, longitude) in radians.

Return type:

Tuple[float, float]