Coordinate Systems: Transformations and Projections

This notebook covers coordinate system conversions and map projections commonly used in tracking and navigation. We explore:

Geodetic Coordinates - Latitude, longitude, altitude (WGS84)
ECEF Coordinates - Earth-Centered Earth-Fixed Cartesian system
Local Frames - ENU (East-North-Up) and NED (North-East-Down)
Rotation Representations - Euler angles, quaternions, rotation matrices
Map Projections - UTM, Mercator, stereographic

Prerequisites

pip install nrl-tracker plotly numpy

[1]:

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from pytcl.coordinate_systems import (
    # Geodetic conversions
    geodetic2ecef, ecef2geodetic,
    ecef2enu, enu2ecef, ecef2ned, ned2ecef,
    enu2ned, ned2enu,
    # Rotation operations
    rotx, roty, rotz,
    euler2rotmat, rotmat2euler,
    euler2quat, quat2euler,
    quat_multiply, quat_rotate, slerp,
    axisangle2rotmat, rotmat2axisangle,
    # Projections
    geodetic2utm, utm2geodetic,
    mercator, mercator_inverse,
    stereographic, polar_stereographic,
    lambert_conformal_conic,
)
from pytcl.core.constants import WGS84

np.random.seed(42)

# Plotly dark theme template
dark_template = go.layout.Template()
dark_template.layout = go.Layout(
    paper_bgcolor='#0d1117',
    plot_bgcolor='#0d1117',
    font=dict(color='#e6edf3'),
    xaxis=dict(gridcolor='#30363d', zerolinecolor='#30363d'),
    yaxis=dict(gridcolor='#30363d', zerolinecolor='#30363d'),
)

1. Geodetic and ECEF Coordinates

Geodetic Coordinates (LLA)

Geodetic coordinates describe position using:

Latitude (φ): Angle from equatorial plane (-90° to +90°)
Longitude (λ): Angle from prime meridian (-180° to +180°)
Altitude (h): Height above reference ellipsoid

ECEF Coordinates

Earth-Centered Earth-Fixed (ECEF) is a Cartesian coordinate system:

Origin: Earth’s center of mass
X-axis: Points to (0°N, 0°E) - equator/prime meridian intersection
Y-axis: Points to (0°N, 90°E) - equator/90° east
Z-axis: Points to North Pole (90°N)

Conversion Formulas

The WGS84 ellipsoid is defined by:

Semi-major axis: a = 6,378,137 m
Flattening: f = 1/298.257223563

[2]:

# Define some example locations
locations = {
    'Washington DC': (38.9072, -77.0369, 0),      # lat, lon (deg), alt (m)
    'London': (51.5074, -0.1278, 0),
    'Tokyo': (35.6762, 139.6503, 0),
    'Sydney': (-33.8688, 151.2093, 0),
    'North Pole': (90.0, 0.0, 0),
    'Equator/Prime': (0.0, 0.0, 0),
}

print("Location Coordinates")
print("=" * 80)
print(f"{'Location':15s} | {'Lat (°)':>10s} | {'Lon (°)':>10s} | "
      f"{'X (km)':>10s} | {'Y (km)':>10s} | {'Z (km)':>10s}")
print("-" * 80)

ecef_coords = {}
for name, (lat_deg, lon_deg, alt) in locations.items():
    lat = np.radians(lat_deg)
    lon = np.radians(lon_deg)

    ecef = geodetic2ecef(lat, lon, alt)
    ecef_coords[name] = ecef

    print(f"{name:15s} | {lat_deg:10.4f} | {lon_deg:10.4f} | "
          f"{ecef[0]/1e3:10.1f} | {ecef[1]/1e3:10.1f} | {ecef[2]/1e3:10.1f}")

Location Coordinates
================================================================================
Location        |    Lat (°) |    Lon (°) |     X (km) |     Y (km) |     Z (km)
--------------------------------------------------------------------------------
Washington DC   |    38.9072 |   -77.0369 |     1114.8 |    -4843.1 |     3984.3
London          |    51.5074 |    -0.1278 |     3978.0 |       -8.9 |     4968.9
Tokyo           |    35.6762 |   139.6503 |    -3953.1 |     3358.4 |     3699.1
Sydney          |   -33.8688 |   151.2093 |    -4646.1 |     2553.2 |    -3534.4
North Pole      |    90.0000 |     0.0000 |        0.0 |        0.0 |     6356.8
Equator/Prime   |     0.0000 |     0.0000 |     6378.1 |        0.0 |        0.0

[3]:

# Visualize locations in 3D ECEF
# Draw Earth ellipsoid (simplified as sphere)
u = np.linspace(0, 2 * np.pi, 50)
v = np.linspace(0, np.pi, 30)
R = WGS84.a / 1e6  # Scale to millions of meters

x_sphere = R * np.outer(np.cos(u), np.sin(v))
y_sphere = R * np.outer(np.sin(u), np.sin(v))
z_sphere = R * np.outer(np.ones(np.size(u)), np.cos(v)) * (1 - WGS84.f)

fig = go.Figure()

# Earth surface
fig.add_trace(
    go.Surface(x=x_sphere, y=y_sphere, z=z_sphere,
               colorscale=[[0, 'lightblue'], [1, 'lightblue']],
               opacity=0.3, showscale=False, name='Earth')
)

# Plot locations
colors = ['#00d4ff', '#ff4757', '#00ff88', '#ffb800', '#a855f7', '#ec4899']
for (name, ecef), color in zip(ecef_coords.items(), colors):
    fig.add_trace(
        go.Scatter3d(x=[ecef[0]/1e6], y=[ecef[1]/1e6], z=[ecef[2]/1e6],
                     mode='markers+text', marker=dict(size=8, color=color),
                     text=[name], textposition='top center',
                     name=name)
    )

# Draw axes
axis_length = 8
fig.add_trace(
    go.Scatter3d(x=[0, axis_length], y=[0, 0], z=[0, 0],
                 mode='lines+text', line=dict(color='red', width=4),
                 text=['', 'X'], textposition='top center',
                 showlegend=False)
)
fig.add_trace(
    go.Scatter3d(x=[0, 0], y=[0, axis_length], z=[0, 0],
                 mode='lines+text', line=dict(color='green', width=4),
                 text=['', 'Y'], textposition='top center',
                 showlegend=False)
)
fig.add_trace(
    go.Scatter3d(x=[0, 0], y=[0, 0], z=[0, axis_length],
                 mode='lines+text', line=dict(color='blue', width=4),
                 text=['', 'Z (North)'], textposition='top center',
                 showlegend=False)
)

fig.update_layout(
    template=dark_template,
    title='ECEF Coordinate System',
    height=600,
    scene=dict(
        xaxis_title='X (million m)',
        yaxis_title='Y (million m)',
        zaxis_title='Z (million m)',
        aspectmode='data',
        xaxis=dict(gridcolor='#30363d'),
        yaxis=dict(gridcolor='#30363d'),
        zaxis=dict(gridcolor='#30363d'),
    )
)
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

[4]:

# Verify round-trip conversion accuracy
print("Round-trip Conversion Accuracy (ECEF → Geodetic → ECEF)")
print("=" * 60)

for name, ecef_orig in ecef_coords.items():
    # ECEF to geodetic
    lat, lon, alt = ecef2geodetic(ecef_orig)

    # Back to ECEF
    ecef_back = geodetic2ecef(lat, lon, alt)

    # Compute error
    error = np.linalg.norm(ecef_orig - ecef_back)

    print(f"{name:15s}: Round-trip error = {error:.6e} m")

Round-trip Conversion Accuracy (ECEF → Geodetic → ECEF)
============================================================
Washington DC  : Round-trip error = 4.656613e-10 m
London         : Round-trip error = 1.818989e-12 m
Tokyo          : Round-trip error = 4.656613e-10 m
Sydney         : Round-trip error = 9.313226e-10 m
North Pole     : Round-trip error = 0.000000e+00 m
Equator/Prime  : Round-trip error = 0.000000e+00 m

2. Local Tangent Plane Frames: ENU and NED

For local operations near a reference point, we use tangent plane coordinate systems:

ENU (East-North-Up)

East: Points east (increasing longitude)
North: Points north (increasing latitude)
Up: Points away from Earth center (opposite gravity)

NED (North-East-Down)

North: Points north
East: Points east
Down: Points toward Earth center (with gravity)

NED is common in aerospace (aircraft body frame), while ENU is common in robotics.

[5]:

# Reference point: Washington DC airport
ref_lat = np.radians(38.9)
ref_lon = np.radians(-77.0)
ref_alt = 0.0

print(f"Reference Point: {np.degrees(ref_lat):.2f}°N, {np.degrees(ref_lon):.2f}°W")
print("\nSimulated aircraft positions relative to reference:")

# Simulate aircraft at various positions
aircraft_offsets = [
    ('Aircraft 1', 1000, 2000, 500),     # 1km E, 2km N, 500m Up
    ('Aircraft 2', -500, 3000, 1000),    # 500m W, 3km N, 1km Up
    ('Aircraft 3', 2000, -1000, 2000),   # 2km E, 1km S, 2km Up
]

print(f"{'Aircraft':12s} | {'East (m)':>10s} | {'North (m)':>10s} | {'Up (m)':>10s}")
print("-" * 55)

ecef_aircraft = []
enu_aircraft = []

for name, east, north, up in aircraft_offsets:
    print(f"{name:12s} | {east:10.0f} | {north:10.0f} | {up:10.0f}")

    # Convert ENU offset to ECEF
    enu = np.array([east, north, up])
    ecef = enu2ecef(enu, ref_lat, ref_lon)

    enu_aircraft.append(enu)
    ecef_aircraft.append(ecef)

Reference Point: 38.90°N, -77.00°W

Simulated aircraft positions relative to reference:
Aircraft     |   East (m) |  North (m) |     Up (m)
-------------------------------------------------------
Aircraft 1   |       1000 |       2000 |        500
Aircraft 2   |       -500 |       3000 |       1000
Aircraft 3   |       2000 |      -1000 |       2000

[6]:

# Convert same positions to NED
print("\nComparison: ENU vs NED frames")
print("=" * 80)
print(f"{'Aircraft':12s} | {'E (m)':>8s} {'N (m)':>8s} {'U (m)':>8s} | "
      f"{'N (m)':>8s} {'E (m)':>8s} {'D (m)':>8s}")
print(f"{'':12s} | {'--- ENU ---':^26s} | {'--- NED ---':^26s}")
print("-" * 80)

for i, (name, _, _, _) in enumerate(aircraft_offsets):
    enu = enu_aircraft[i]
    ned = enu2ned(enu)

    print(f"{name:12s} | {enu[0]:8.0f} {enu[1]:8.0f} {enu[2]:8.0f} | "
          f"{ned[0]:8.0f} {ned[1]:8.0f} {ned[2]:8.0f}")

print("\nNote: NED swaps E↔N and negates Up to get Down")


Comparison: ENU vs NED frames
================================================================================
Aircraft     |    E (m)    N (m)    U (m) |    N (m)    E (m)    D (m)
             |        --- ENU ---         |        --- NED ---
--------------------------------------------------------------------------------
Aircraft 1   |     1000     2000      500 |     2000     1000     -500
Aircraft 2   |     -500     3000     1000 |     3000     -500    -1000
Aircraft 3   |     2000    -1000     2000 |    -1000     2000    -2000

Note: NED swaps E↔N and negates Up to get Down

[7]:

# Visualize ENU and NED frames
fig = make_subplots(rows=1, cols=2, subplot_titles=('ENU Frame (looking down)', 'NED Frame (X=North, Y=East)'),
                    horizontal_spacing=0.1)

colors = ['#00d4ff', '#ff4757', '#00ff88']

# ENU plot
for i, (name, _, _, _) in enumerate(aircraft_offsets):
    enu = enu_aircraft[i]
    fig.add_trace(
        go.Scatter(x=[enu[0]], y=[enu[1]], mode='markers+text',
                   marker=dict(size=12, color=colors[i]),
                   text=[f"{name}<br>Up={enu[2]:.0f}m"], textposition='top right',
                   name=name, showlegend=True),
        row=1, col=1
    )

# Reference point for ENU
fig.add_trace(
    go.Scatter(x=[0], y=[0], mode='markers',
               marker=dict(size=15, color='white', symbol='star'),
               name='Reference', showlegend=True),
    row=1, col=1
)

# NED plot
for i, (name, _, _, _) in enumerate(aircraft_offsets):
    enu = enu_aircraft[i]
    ned = enu2ned(enu)
    fig.add_trace(
        go.Scatter(x=[ned[1]], y=[ned[0]], mode='markers+text',
                   marker=dict(size=12, color=colors[i]),
                   text=[f"Down={ned[2]:.0f}m"], textposition='top right',
                   showlegend=False),
        row=1, col=2
    )

# Reference point for NED
fig.add_trace(
    go.Scatter(x=[0], y=[0], mode='markers',
               marker=dict(size=15, color='white', symbol='star'),
               showlegend=False),
    row=1, col=2
)

fig.update_layout(
    template=dark_template,
    height=450,
)
fig.update_xaxes(title_text='East (m)', row=1, col=1, scaleanchor='y', scaleratio=1)
fig.update_xaxes(title_text='East (m)', row=1, col=2)
fig.update_yaxes(title_text='North (m)', row=1, col=1)
fig.update_yaxes(title_text='North (m)', row=1, col=2, scaleanchor='x2', scaleratio=1)

fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

3. Rotation Representations

Rotations in 3D can be represented several ways:

Representation	Parameters	Singularities	Use Case
Euler Angles	3	Gimbal lock	Human intuition
Rotation Matrix	9 (6 DOF)	None	Direct application
Quaternion	4 (3 DOF)	None	Interpolation, composition
Axis-Angle	3	Small angles	Visualization

Euler Angles (ZYX Convention)

The aerospace convention uses:

Yaw (ψ): Rotation about Z-axis (heading)
Pitch (θ): Rotation about Y-axis (nose up/down)
Roll (φ): Rotation about X-axis (bank)

[8]:

# Demonstrate rotation representations
yaw = np.radians(45)    # 45° heading
pitch = np.radians(15)  # 15° nose up
roll = np.radians(10)   # 10° bank right

angles = np.array([yaw, pitch, roll])

print("Input Euler angles (ZYX convention):")
print(f"  Yaw:   {np.degrees(yaw):6.1f}°")
print(f"  Pitch: {np.degrees(pitch):6.1f}°")
print(f"  Roll:  {np.degrees(roll):6.1f}°")

# Convert to rotation matrix
R = euler2rotmat(angles, 'ZYX')
print(f"\nRotation Matrix:\n{R}")

# Convert to quaternion
q = euler2quat(angles, 'ZYX')
print(f"\nQuaternion [w, x, y, z]: {q}")

# Convert back to verify
angles_back = rotmat2euler(R, 'ZYX')
print(f"\nRecovered Euler angles:")
print(f"  Yaw:   {np.degrees(angles_back[0]):6.1f}°")
print(f"  Pitch: {np.degrees(angles_back[1]):6.1f}°")
print(f"  Roll:  {np.degrees(angles_back[2]):6.1f}°")

Input Euler angles (ZYX convention):
  Yaw:     45.0°
  Pitch:   15.0°
  Roll:    10.0°

Rotation Matrix:
[[ 0.6830127  -0.66458442  0.30302013]
 [ 0.6830127   0.72814406  0.05744452]
 [-0.25881905  0.16773126  0.95125124]]

Quaternion [w, x, y, z]: [0.9168435  0.0300724  0.15319931 0.3674556 ]

Recovered Euler angles:
  Yaw:     45.0°
  Pitch:   15.0°
  Roll:    10.0°

[9]:

# Visualize rotation effect on aircraft axes
# Define aircraft body axes (before rotation)
body_axes = np.array([
    [1, 0, 0],  # Forward (nose)
    [0, 1, 0],  # Right wing
    [0, 0, 1],  # Down
]).T

# Apply rotation
rotated_axes = R @ body_axes

fig = make_subplots(rows=1, cols=2, specs=[[{'type': 'scene'}, {'type': 'scene'}]],
                    subplot_titles=['Body Axes (No Rotation)',
                                   f'Rotated (Yaw={np.degrees(yaw):.0f}°, Pitch={np.degrees(pitch):.0f}°, Roll={np.degrees(roll):.0f}°)'])

colors = ['red', 'green', 'blue']
labels = ['Forward (X)', 'Right (Y)', 'Down (Z)']

# Plot before rotation (left)
for i, (color, label) in enumerate(zip(colors, labels)):
    fig.add_trace(
        go.Scatter3d(x=[0, body_axes[0, i]], y=[0, body_axes[1, i]], z=[0, body_axes[2, i]],
                     mode='lines', line=dict(color=color, width=6),
                     name=label),
        row=1, col=1
    )

# Plot after rotation (right)
for i, (color, label) in enumerate(zip(colors, labels)):
    fig.add_trace(
        go.Scatter3d(x=[0, rotated_axes[0, i]], y=[0, rotated_axes[1, i]], z=[0, rotated_axes[2, i]],
                     mode='lines', line=dict(color=color, width=6),
                     showlegend=False),
        row=1, col=2
    )

scene_layout = dict(
    xaxis=dict(range=[-1, 1], gridcolor='#30363d'),
    yaxis=dict(range=[-1, 1], gridcolor='#30363d'),
    zaxis=dict(range=[-1, 1], gridcolor='#30363d'),
    aspectmode='cube'
)

fig.update_layout(
    template=dark_template,
    height=450,
    scene=scene_layout,
    scene2=scene_layout,
)
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

4. Quaternion Operations

Quaternions are excellent for:

Composition: Multiplying quaternions combines rotations
Interpolation: SLERP provides smooth rotation interpolation
Numerical stability: No gimbal lock issues

[10]:

# Demonstrate quaternion operations

# Two rotations to compose
q1 = euler2quat(np.radians([30, 0, 0]), 'ZYX')  # 30° yaw
q2 = euler2quat(np.radians([0, 20, 0]), 'ZYX')  # 20° pitch

print("Composing rotations:")
print(f"  q1 (30° yaw):   {q1}")
print(f"  q2 (20° pitch): {q2}")

# Multiply quaternions (composition)
q_composed = quat_multiply(q1, q2)
print(f"  q1 * q2:        {q_composed}")

# Convert result back to Euler
euler_composed = quat2euler(q_composed, 'ZYX')
print(f"\nComposed rotation (Euler):")
print(f"  Yaw: {np.degrees(euler_composed[0]):.1f}°, "
      f"Pitch: {np.degrees(euler_composed[1]):.1f}°, "
      f"Roll: {np.degrees(euler_composed[2]):.1f}°")

Composing rotations:
  q1 (30° yaw):   [0.96592583 0.         0.         0.25881905]
  q2 (20° pitch): [0.98480775 0.         0.17364818 0.        ]
  q1 * q2:        [ 0.95125124 -0.04494346  0.16773126  0.254887  ]

Composed rotation (Euler):
  Yaw: 30.0°, Pitch: 20.0°, Roll: 0.0°

[11]:

# SLERP interpolation demonstration
q_start = euler2quat(np.radians([0, 0, 0]), 'ZYX')    # No rotation
q_end = euler2quat(np.radians([180, 45, 30]), 'ZYX')  # Complex rotation

# Interpolate
t_values = np.linspace(0, 1, 11)
interpolated = []

print("SLERP Interpolation:")
print(f"{'t':>5s} | {'Yaw (°)':>10s} | {'Pitch (°)':>10s} | {'Roll (°)':>10s}")
print("-" * 45)

for t in t_values:
    q_interp = slerp(q_start, q_end, t)
    euler = quat2euler(q_interp, 'ZYX')
    interpolated.append(np.degrees(euler))
    print(f"{t:5.2f} | {np.degrees(euler[0]):10.2f} | "
          f"{np.degrees(euler[1]):10.2f} | {np.degrees(euler[2]):10.2f}")

interpolated = np.array(interpolated)

SLERP Interpolation:
    t |    Yaw (°) |  Pitch (°) |   Roll (°)
---------------------------------------------
 0.00 |       0.00 |      -0.00 |       0.00
 0.10 |      14.91 |       4.82 |      -5.67
 0.20 |      29.46 |      10.92 |      -9.97
 0.30 |      43.94 |      17.90 |     -12.65
 0.40 |      58.79 |      25.30 |     -13.44
 0.50 |      74.62 |      32.63 |     -12.03
 0.60 |      92.16 |      39.31 |      -7.96
 0.70 |     112.04 |      44.57 |      -0.89
 0.80 |     134.26 |      47.60 |       8.92
 0.90 |     157.59 |      47.76 |      19.98
 1.00 |     180.00 |      45.00 |      30.00

[12]:

# Visualize SLERP interpolation
fig = go.Figure()

fig.add_trace(
    go.Scatter(x=t_values, y=interpolated[:, 0], mode='lines+markers',
               name='Yaw', line=dict(color='#00d4ff', width=2),
               marker=dict(size=8))
)
fig.add_trace(
    go.Scatter(x=t_values, y=interpolated[:, 1], mode='lines+markers',
               name='Pitch', line=dict(color='#ff4757', width=2),
               marker=dict(size=8, symbol='square'))
)
fig.add_trace(
    go.Scatter(x=t_values, y=interpolated[:, 2], mode='lines+markers',
               name='Roll', line=dict(color='#00ff88', width=2),
               marker=dict(size=8, symbol='triangle-up'))
)

fig.update_layout(
    template=dark_template,
    title='SLERP Quaternion Interpolation',
    xaxis_title='Interpolation Parameter t',
    yaxis_title='Angle (degrees)',
    height=400,
)
fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

5. Map Projections

Map projections convert 3D geodetic coordinates to 2D plane coordinates. Different projections preserve different properties:

Projection	Preserves	Distorts	Use Case
Mercator	Angles (conformal)	Area	Navigation, web maps
UTM	Local distances	At zone edges	Military, surveying
Lambert	Angles (conformal)	Poles	Mid-latitude regions
Stereographic	Angles (conformal)	Edges	Polar regions

[13]:

# UTM projection example
print("UTM Projection Examples")
print("=" * 70)

for name, (lat_deg, lon_deg, _) in locations.items():
    lat = np.radians(lat_deg)
    lon = np.radians(lon_deg)

    # Skip poles (UTM not defined there)
    if abs(lat_deg) > 84:
        print(f"{name:15s}: UTM not defined (use UPS for polar regions)")
        continue

    utm = geodetic2utm(lat, lon)

    print(f"{name:15s}: Zone {utm.zone:2d}{utm.hemisphere}, "
          f"E={utm.easting:10.1f}m, N={utm.northing:11.1f}m, "
          f"Scale={utm.scale:.6f}")

UTM Projection Examples
======================================================================
Washington DC  : Zone 18N, E=  323383.2m, N=  4318900.7m, Scale=0.999984
London         : Zone 30N, E=  699316.2m, N=  5720598.4m, Scale=1.000088
Tokyo          : Zone 54N, E=  377855.8m, N=  3958999.5m, Scale=0.999784
Sydney         : Zone 56S, E=  334368.6m, N=  6241061.1m, Scale=0.999938
North Pole     : UTM not defined (use UPS for polar regions)
Equator/Prime  : Zone 31N, E=  166021.4m, N=        0.0m, Scale=1.000981

[14]:

# Compare projections for a region
# Generate a grid of points around Washington DC
center_lat = np.radians(38.9)
center_lon = np.radians(-77.0)

# Create a 10° x 10° grid
lat_range = np.linspace(center_lat - np.radians(5), center_lat + np.radians(5), 20)
lon_range = np.linspace(center_lon - np.radians(5), center_lon + np.radians(5), 20)

LAT, LON = np.meshgrid(lat_range, lon_range)

# Project using different methods
mercator_x = np.zeros_like(LAT)
mercator_y = np.zeros_like(LAT)
stereo_x = np.zeros_like(LAT)
stereo_y = np.zeros_like(LAT)

for i in range(LAT.shape[0]):
    for j in range(LAT.shape[1]):
        # Mercator
        merc = mercator(LAT[i, j], LON[i, j], center_lon)
        mercator_x[i, j] = merc.x / 1e6
        mercator_y[i, j] = merc.y / 1e6

        # Stereographic
        ster = stereographic(LAT[i, j], LON[i, j], center_lat, center_lon)
        stereo_x[i, j] = ster.x / 1e6
        stereo_y[i, j] = ster.y / 1e6

[15]:

# Visualize projection comparison
fig = make_subplots(rows=1, cols=3,
                    subplot_titles=['Geographic (Lat/Lon)', 'Mercator Projection', 'Stereographic Projection'],
                    horizontal_spacing=0.08)

# Original lat/lon grid
for i in range(LAT.shape[0]):
    fig.add_trace(
        go.Scatter(x=np.degrees(LON[i, :]), y=np.degrees(LAT[i, :]),
                   mode='lines', line=dict(color='#00d4ff', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=1
    )
for j in range(LAT.shape[1]):
    fig.add_trace(
        go.Scatter(x=np.degrees(LON[:, j]), y=np.degrees(LAT[:, j]),
                   mode='lines', line=dict(color='#00d4ff', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=1
    )
fig.add_trace(
    go.Scatter(x=[np.degrees(center_lon)], y=[np.degrees(center_lat)],
               mode='markers', marker=dict(color='#ff4757', size=10),
               name='Center', showlegend=True),
    row=1, col=1
)

# Mercator
for i in range(LAT.shape[0]):
    fig.add_trace(
        go.Scatter(x=mercator_x[i, :], y=mercator_y[i, :],
                   mode='lines', line=dict(color='#00ff88', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=2
    )
for j in range(LAT.shape[1]):
    fig.add_trace(
        go.Scatter(x=mercator_x[:, j], y=mercator_y[:, j],
                   mode='lines', line=dict(color='#00ff88', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=2
    )
fig.add_trace(
    go.Scatter(x=[0], y=[mercator_y[10, 10]],
               mode='markers', marker=dict(color='#ff4757', size=10),
               showlegend=False),
    row=1, col=2
)

# Stereographic
for i in range(LAT.shape[0]):
    fig.add_trace(
        go.Scatter(x=stereo_x[i, :], y=stereo_y[i, :],
                   mode='lines', line=dict(color='#a855f7', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=3
    )
for j in range(LAT.shape[1]):
    fig.add_trace(
        go.Scatter(x=stereo_x[:, j], y=stereo_y[:, j],
                   mode='lines', line=dict(color='#a855f7', width=1), opacity=0.5,
                   showlegend=False),
        row=1, col=3
    )
fig.add_trace(
    go.Scatter(x=[0], y=[0],
               mode='markers', marker=dict(color='#ff4757', size=10),
               showlegend=False),
    row=1, col=3
)

fig.update_layout(
    template=dark_template,
    height=400,
)
fig.update_xaxes(title_text='Longitude (°)', row=1, col=1, scaleanchor='y', scaleratio=1)
fig.update_xaxes(title_text='Easting (million m)', row=1, col=2, scaleanchor='y2', scaleratio=1)
fig.update_xaxes(title_text='Easting (million m)', row=1, col=3, scaleanchor='y3', scaleratio=1)
fig.update_yaxes(title_text='Latitude (°)', row=1, col=1)
fig.update_yaxes(title_text='Northing (million m)', row=1, col=2)
fig.update_yaxes(title_text='Northing (million m)', row=1, col=3)

fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

[16]:

# UTM zone visualization
print("UTM Zone Coverage")
print("=" * 50)

# Show UTM zones for various longitudes
longitudes = np.arange(-180, 181, 30)
print(f"{'Longitude':>12s} | {'UTM Zone':>10s} | {'Central Meridian':>18s}")
print("-" * 50)

for lon_deg in longitudes:
    lon = np.radians(lon_deg)
    lat = 0  # Equator

    utm_result = geodetic2utm(lat, lon)

    from pytcl.coordinate_systems import utm_central_meridian
    central = np.degrees(utm_central_meridian(utm_result.zone))

    print(f"{lon_deg:12.0f}° | {utm_result.zone:10d} | {central:18.0f}°")

UTM Zone Coverage
==================================================
   Longitude |   UTM Zone |   Central Meridian
--------------------------------------------------
        -180° |          1 |               -177°
        -150° |          6 |               -147°
        -120° |         11 |               -117°
         -90° |         16 |                -87°
         -60° |         21 |                -57°
         -30° |         26 |                -27°
           0° |         31 |                  3°
          30° |         36 |                 33°
          60° |         41 |                 63°
          90° |         46 |                 93°
         120° |         51 |                123°
         150° |         56 |                153°
         180° |          1 |               -177°

6. Practical Example: Aircraft Tracking

Let’s put it all together with a complete aircraft tracking example.

[17]:

# Simulate an aircraft trajectory
# Takeoff from Washington DC, fly northeast

# Radar station location
radar_lat = np.radians(38.9)
radar_lon = np.radians(-77.0)
radar_alt = 100.0  # 100m tower

# Generate flight path (geodetic)
n_points = 50
t = np.linspace(0, 1, n_points)

# Start position
start_lat = np.radians(38.95)
start_lon = np.radians(-77.05)

# End position (100km northeast)
end_lat = np.radians(39.7)
end_lon = np.radians(-76.2)

# Interpolate position
flight_lat = start_lat + t * (end_lat - start_lat)
flight_lon = start_lon + t * (end_lon - start_lon)
flight_alt = 1000 + 9000 * np.sin(np.pi * t)  # Climb to 10km, then descend

print(f"Flight from ({np.degrees(start_lat):.2f}°, {np.degrees(start_lon):.2f}°) "
      f"to ({np.degrees(end_lat):.2f}°, {np.degrees(end_lon):.2f}°)")
print(f"Max altitude: {max(flight_alt):.0f} m")

Flight from (38.95°, -77.05°) to (39.70°, -76.20°)
Max altitude: 9995 m

[18]:

# Convert to various coordinate systems
flight_ecef = np.array([geodetic2ecef(lat, lon, alt)
                        for lat, lon, alt in zip(flight_lat, flight_lon, flight_alt)])

flight_enu = np.array([ecef2enu(ecef, radar_lat, radar_lon)
                       for ecef in flight_ecef])

flight_ned = np.array([ecef2ned(ecef, radar_lat, radar_lon)
                       for ecef in flight_ecef])

# Compute range and angles from radar
ranges = np.linalg.norm(flight_enu, axis=1)
azimuths = np.degrees(np.arctan2(flight_enu[:, 0], flight_enu[:, 1]))  # From north
elevations = np.degrees(np.arcsin(flight_enu[:, 2] / ranges))

print("\nFlight Path Statistics:")
print(f"  Range: {ranges[0]/1e3:.1f} km to {ranges.max()/1e3:.1f} km")
print(f"  Azimuth: {azimuths[0]:.1f}° to {azimuths[-1]:.1f}°")
print(f"  Elevation: {elevations[0]:.1f}° to {elevations.max():.1f}°")


Flight Path Statistics:
  Range: 7.1 km to 112.5 km
  Azimuth: -38.0° to 37.6°
  Elevation: 8.0° to 14.8°

[19]:

# Comprehensive visualization
fig = make_subplots(
    rows=2, cols=2,
    specs=[[{'type': 'scene'}, {'type': 'xy'}],
           [{'type': 'xy'}, {'type': 'xy'}]],
    subplot_titles=['ENU Coordinates from Radar', 'Overhead View',
                    'Range-Altitude Profile', 'Radar Angles'],
    vertical_spacing=0.12,
    horizontal_spacing=0.1
)

# 3D ENU view
fig.add_trace(
    go.Scatter3d(x=flight_enu[:, 0]/1e3, y=flight_enu[:, 1]/1e3, z=flight_enu[:, 2]/1e3,
                 mode='lines', line=dict(color='#00d4ff', width=4),
                 name='Flight path'),
    row=1, col=1
)
fig.add_trace(
    go.Scatter3d(x=[0], y=[0], z=[0.1],
                 mode='markers', marker=dict(color='#ff4757', size=8, symbol='diamond'),
                 name='Radar'),
    row=1, col=1
)
fig.add_trace(
    go.Scatter3d(x=[flight_enu[0, 0]/1e3], y=[flight_enu[0, 1]/1e3], z=[flight_enu[0, 2]/1e3],
                 mode='markers', marker=dict(color='#00ff88', size=8),
                 name='Takeoff'),
    row=1, col=1
)
fig.add_trace(
    go.Scatter3d(x=[flight_enu[-1, 0]/1e3], y=[flight_enu[-1, 1]/1e3], z=[flight_enu[-1, 2]/1e3],
                 mode='markers', marker=dict(color='#ffb800', size=8, symbol='square'),
                 name='Landing'),
    row=1, col=1
)

# 2D overhead view with altitude coloring
fig.add_trace(
    go.Scatter(x=flight_enu[:, 0]/1e3, y=flight_enu[:, 1]/1e3,
               mode='markers+lines',
               marker=dict(color=flight_alt/1e3, colorscale='Viridis', size=6,
                          colorbar=dict(title='Alt (km)', x=0.45, len=0.4, y=0.85)),
               line=dict(color='rgba(0,212,255,0.3)', width=1),
               name='Flight', showlegend=False),
    row=1, col=2
)
fig.add_trace(
    go.Scatter(x=[0], y=[0], mode='markers',
               marker=dict(color='#ff4757', size=12, symbol='triangle-up'),
               showlegend=False),
    row=1, col=2
)

# Range-altitude profile
horizontal_range = np.sqrt(flight_enu[:, 0]**2 + flight_enu[:, 1]**2) / 1e3
fig.add_trace(
    go.Scatter(x=horizontal_range, y=flight_alt/1e3, mode='lines',
               fill='tozeroy', fillcolor='rgba(0,212,255,0.3)',
               line=dict(color='#00d4ff', width=2),
               name='Altitude', showlegend=False),
    row=2, col=1
)
fig.add_trace(
    go.Scatter(x=[horizontal_range[0], horizontal_range[-1]], y=[0, 0],
               mode='lines', line=dict(color='#8B4513', width=4),
               name='Ground', showlegend=False),
    row=2, col=1
)

# Radar angles
fig.add_trace(
    go.Scatter(x=t * 100, y=azimuths, mode='lines',
               name='Azimuth', line=dict(color='#00ff88', width=2)),
    row=2, col=2
)
fig.add_trace(
    go.Scatter(x=t * 100, y=elevations, mode='lines',
               name='Elevation', line=dict(color='#ff4757', width=2)),
    row=2, col=2
)

fig.update_layout(
    template=dark_template,
    height=700,
    scene=dict(
        xaxis_title='East (km)',
        yaxis_title='North (km)',
        zaxis_title='Up (km)',
        xaxis=dict(gridcolor='#30363d'),
        yaxis=dict(gridcolor='#30363d'),
        zaxis=dict(gridcolor='#30363d'),
    ),
)

fig.update_xaxes(title_text='East (km)', row=1, col=2)
fig.update_yaxes(title_text='North (km)', row=1, col=2, scaleanchor='x2', scaleratio=1)
fig.update_xaxes(title_text='Horizontal Range (km)', row=2, col=1)
fig.update_yaxes(title_text='Altitude (km)', row=2, col=1)
fig.update_xaxes(title_text='Flight Progress (%)', row=2, col=2)
fig.update_yaxes(title_text='Angle (degrees)', row=2, col=2)

fig.show()

Data type cannot be displayed: application/vnd.plotly.v1+json

Summary

Key takeaways:

Geodetic coordinates are intuitive but require conversion for calculations
ECEF is ideal for global calculations but unintuitive for local work
ENU/NED local frames are best for tracking near a reference point
Quaternions avoid gimbal lock and provide smooth interpolation
Map projections trade off different properties - choose based on use case

Exercises

Implement a coordinate converter that handles all transitions between geodetic, ECEF, ENU, and NED
Compare quaternion SLERP to linear interpolation of Euler angles
Visualize how projection distortion varies with distance from the projection center
Track a satellite in ECEF and display its ground track on a map projection

References

Groves, P. D. (2013). Principles of GNSS, Inertial, and Multisensor Integrated Navigation Systems
Snyder, J. P. (1987). Map Projections: A Working Manual. USGS Professional Paper 1395.
Kuipers, J. B. (1999). Quaternions and Rotation Sequences.