Skip to content

Topocentric Coordinate Transformations

Topocentric coordinate systems are local horizon-based reference frames centered on an observer, such as a ground station or radar site. These coordinate systems are essential for satellite tracking, visibility analysis, and determining where to point antennas or telescopes.

Unlike global coordinate systems (ECEF, ECI), topocentric systems define positions relative to a specific location on Earth, making it easy to determine whether a satellite is visible and where to look in the sky.

For complete API details, see the Topocentric Coordinates API Reference.

Topocentric Coordinate Systems

Brahe supports two local horizon coordinate systems:

ENZ (East-North-Zenith)

  • East (E): Positive toward geographic east
  • North (N): Positive toward geographic north
  • Zenith (Z): Positive upward (toward the sky)

This is the most common topocentric system for satellite tracking and is aligned with geographic directions.

SEZ (South-East-Zenith)

  • South (S): Positive toward geographic south
  • East (E): Positive toward geographic east
  • Zenith (Z): Positive upward (toward the sky)

The SEZ system is sometimes used in radar and missile tracking applications. The main difference from ENZ is that the first two axes are rotated 180° around the zenith axis.

Info

Both ENZ and SEZ use a right-handed coordinate system with the zenith axis pointing up. The choice between them is typically driven by convention in your specific field or application.

Station Location Interpretation

When specifying the observer (ground station) location, you must choose whether the coordinates represent:

  • Geodetic (EllipsoidalConversionType.GEODETIC): Station coordinates use WGS84 ellipsoid (recommended for accuracy)
  • Geocentric (EllipsoidalConversionType.GEOCENTRIC): Station coordinates use spherical Earth model

For ground stations, geodetic interpretation is almost always preferred for accuracy.

ENZ Transformations

Converting ECEF to ENZ

To get the position of an object relative to a location, you need to convert the object's ECEF position to the local ENZ frame centered on the location:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define ground station location in geodetic coordinates
# Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
lon_deg = -122.17329
lat_deg = 37.42692
alt_m = 32.0

print("Ground Station (Stanford):")
print(f"Longitude: {lon_deg:.5f}° = {np.radians(lon_deg):.6f} rad")
print(f"Latitude:  {lat_deg:.5f}° = {np.radians(lat_deg):.6f} rad")
print(f"Altitude:  {alt_m:.1f} m\n")
# Longitude: -122.17329° = -2.132605 rad
# Latitude:  37.42692° = 0.653131 rad
# Altitude:  32.0 m

# Convert ground station to ECEF
geodetic_station = np.array([lon_deg, lat_deg, alt_m])
station_ecef = bh.position_geodetic_to_ecef(geodetic_station, bh.AngleFormat.DEGREES)

print("Ground Station ECEF:")
print(f"x = {station_ecef[0]:.3f} m")
print(f"y = {station_ecef[1]:.3f} m")
print(f"z = {station_ecef[2]:.3f} m\n")
# x = -2700691.122 m
# y = -4292566.016 m
# z = 3855395.780 m

# Define satellite in sun-synchronous orbit at 500 km altitude
# SSO orbit passes over Stanford at approximately 10:30 AM local time
# Orbital elements: [a, e, i, RAAN, omega, M]
oe = np.array([bh.R_EARTH + 500e3, 0.001, 97.8, 240.0, 0.0, 90.0])

# Define epoch when satellite passes near Stanford (Jan 1, 2024, 17:05 UTC)
epoch = bh.Epoch.from_datetime(2024, 1, 1, 17, 5, 0.0, 0.0, bh.TimeSystem.UTC)

# Convert orbital elements to ECI state
sat_state_eci = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Convert ECI state to ECEF at the given epoch
sat_state_ecef = bh.state_eci_to_ecef(epoch, sat_state_eci)
sat_ecef = sat_state_ecef[0:3]  # Extract position only

year, month, day, hour, minute, second, ns = epoch.to_datetime()
print(f"Epoch: {year}-{month:02d}-{day:02d} {hour:02d}:{minute:02d}:{second:06.3f} UTC")
print("Satellite ECEF:")
print(f"x = {sat_ecef[0]:.3f} m")
print(f"y = {sat_ecef[1]:.3f} m")
print(f"z = {sat_ecef[2]:.3f} m\n")

# Convert satellite position to ENZ coordinates relative to ground station
enz = bh.relative_position_ecef_to_enz(
    station_ecef, sat_ecef, bh.EllipsoidalConversionType.GEODETIC
)

print("Satellite position in ENZ frame (relative to Stanford):")
print(f"East:   {enz[0] / 1000:.3f} km")
print(f"North:  {enz[1] / 1000:.3f} km")
print(f"Zenith: {enz[2] / 1000:.3f} km")
print(f"Range:  {np.linalg.norm(enz) / 1000:.3f} km")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define ground station location in geodetic coordinates
    // Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
    let lon_deg = -122.17329_f64;
    let lat_deg = 37.42692_f64;
    let alt_m = 32.0;

    println!("Ground Station (Stanford):");
    println!("Longitude: {:.5}° = {:.6} rad", lon_deg, lon_deg.to_radians());
    println!("Latitude:  {:.5}° = {:.6} rad", lat_deg, lat_deg.to_radians());
    println!("Altitude:  {:.1} m\n", alt_m);
    // Expected output:
    // Longitude: -122.17329° = -2.132605 rad
    // Latitude:  37.42692° = 0.653131 rad
    // Altitude:  32.0 m

    // Convert ground station to ECEF
    let geodetic_station = na::Vector3::new(lon_deg, lat_deg, alt_m);
    let station_ecef = bh::position_geodetic_to_ecef(geodetic_station, bh::AngleFormat::Degrees).unwrap();

    println!("Ground Station ECEF:");
    println!("x = {:.3} m", station_ecef[0]);
    println!("y = {:.3} m", station_ecef[1]);
    println!("z = {:.3} m\n", station_ecef[2]);
    // Expected output:
    // x = -2700691.122 m
    // y = -4292566.016 m
    // z = 3855395.780 m

    // Define satellite in sun-synchronous orbit at 500 km altitude
    // SSO orbit passes over Stanford at approximately 10:30 AM local time
    // Orbital elements: [a, e, i, RAAN, omega, M]
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,
        0.001,
        97.8_f64,
        240.0_f64,
        0.0_f64,
        90.0_f64,
    );

    // Define epoch when satellite passes near Stanford (Jan 1, 2024, 17:05 UTC)
    let epoch = bh::Epoch::from_datetime(2024, 1, 1, 17, 5, 0.0, 0.0, bh::TimeSystem::UTC);

    // Convert orbital elements to ECI state
    let sat_state_eci = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);

    // Convert ECI state to ECEF at the given epoch
    let sat_state_ecef = bh::state_eci_to_ecef(epoch, sat_state_eci);
    let sat_ecef = na::Vector3::new(sat_state_ecef[0], sat_state_ecef[1], sat_state_ecef[2]);

    let (year, month, day, hour, minute, second, _ns) = epoch.to_datetime();
    println!("Epoch: {}-{:02}-{:02} {:02}:{:02}:{:06.3} UTC", year, month, day, hour, minute, second);
    println!("Satellite ECEF:");
    println!("x = {:.3} m", sat_ecef[0]);
    println!("y = {:.3} m", sat_ecef[1]);
    println!("z = {:.3} m\n", sat_ecef[2]);

    // Convert satellite position to ENZ coordinates relative to ground station
    let enz = bh::relative_position_ecef_to_enz(
        station_ecef,
        sat_ecef,
        bh::EllipsoidalConversionType::Geodetic,
    );

    println!("Satellite position in ENZ frame (relative to Stanford):");
    println!("East:   {:.3} km", enz[0] / 1000.0);
    println!("North:  {:.3} km", enz[1] / 1000.0);
    println!("Zenith: {:.3} km", enz[2] / 1000.0);
    println!("Range:  {:.3} km", enz.norm() / 1000.0);
}

Converting ENZ to ECEF

The reverse transformation converts a relative ENZ position back to an absolute ECEF position:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define ground station location in geodetic coordinates
# Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
lon_deg = -122.17329
lat_deg = 37.42692
alt_m = 32.0

print("Ground Station (Stanford):")
print(f"Longitude: {lon_deg:.5f}° = {np.radians(lon_deg):.6f} rad")
print(f"Latitude:  {lat_deg:.5f}° = {np.radians(lat_deg):.6f} rad")
print(f"Altitude:  {alt_m:.1f} m\n")
# Longitude: -122.17329° = -2.132605 rad
# Latitude:  37.42692° = 0.653131 rad
# Altitude:  32.0 m

# Convert ground station to ECEF
geodetic_station = np.array([lon_deg, lat_deg, alt_m])
station_ecef = bh.position_geodetic_to_ecef(geodetic_station, bh.AngleFormat.DEGREES)

print("Ground Station ECEF:")
print(f"x = {station_ecef[0]:.3f} m")
print(f"y = {station_ecef[1]:.3f} m")
print(f"z = {station_ecef[2]:.3f} m\n")
# x = -2700691.122 m
# y = -4292566.016 m
# z = 3855395.780 m

# Define relative position in ENZ coordinates
# Example: 50 km East, 100 km North, 200 km Up from station
enz = np.array([50e3, 100e3, 200e3])

print("Relative position in ENZ frame:")
print(f"East:   {enz[0] / 1000:.1f} km")
print(f"North:  {enz[1] / 1000:.1f} km")
print(f"Zenith: {enz[2] / 1000:.1f} km\n")
# East:   50.0 km
# North:  100.0 km
# Zenith: 200.0 km

# Convert ENZ relative position to absolute ECEF position
target_ecef = bh.relative_position_enz_to_ecef(
    station_ecef, enz, bh.EllipsoidalConversionType.GEODETIC
)

print("Target position in ECEF:")
print(f"x = {target_ecef[0]:.3f} m")
print(f"y = {target_ecef[1]:.3f} m")
print(f"z = {target_ecef[2]:.3f} m")
print(f"Distance from Earth center: {np.linalg.norm(target_ecef) / 1000:.3f} km")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define ground station location in geodetic coordinates
    // Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
    let lon_deg = -122.17329_f64;
    let lat_deg = 37.42692_f64;
    let alt_m = 32.0;

    println!("Ground Station (Stanford):");
    println!("Longitude: {:.5}° = {:.6} rad", lon_deg, lon_deg.to_radians());
    println!("Latitude:  {:.5}° = {:.6} rad", lat_deg, lat_deg.to_radians());
    println!("Altitude:  {:.1} m\n", alt_m);
    // Expected output:
    // Longitude: -122.17329° = -2.132605 rad
    // Latitude:  37.42692° = 0.653131 rad
    // Altitude:  32.0 m

    // Convert ground station to ECEF
    let geodetic_station = na::Vector3::new(lon_deg, lat_deg, alt_m);
    let station_ecef = bh::position_geodetic_to_ecef(geodetic_station, bh::AngleFormat::Degrees).unwrap();

    println!("Ground Station ECEF:");
    println!("x = {:.3} m", station_ecef[0]);
    println!("y = {:.3} m", station_ecef[1]);
    println!("z = {:.3} m\n", station_ecef[2]);
    // Expected output:
    // x = -2700691.122 m
    // y = -4292566.016 m
    // z = 3855395.780 m

    // Define relative position in ENZ coordinates
    // Example: 50 km East, 100 km North, 200 km Up from station
    let enz = na::Vector3::new(50e3, 100e3, 200e3);

    println!("Relative position in ENZ frame:");
    println!("East:   {:.1} km", enz[0] / 1000.0);
    println!("North:  {:.1} km", enz[1] / 1000.0);
    println!("Zenith: {:.1} km\n", enz[2] / 1000.0);
    // Expected output:
    // East:   50.0 km
    // North:  100.0 km
    // Zenith: 200.0 km

    // Convert ENZ relative position to absolute ECEF position
    let target_ecef = bh::relative_position_enz_to_ecef(
        station_ecef,
        enz,
        bh::EllipsoidalConversionType::Geodetic,
    );

    println!("Target position in ECEF:");
    println!("x = {:.3} m", target_ecef[0]);
    println!("y = {:.3} m", target_ecef[1]);
    println!("z = {:.3} m", target_ecef[2]);
    println!("Distance from Earth center: {:.3} km", target_ecef.norm() / 1000.0);
}

SEZ Transformations

Converting ECEF to SEZ

Similar to ENZ, you can convert ECEF positions to the SEZ frame:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define ground station location in geodetic coordinates
# Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
lon_deg = -122.17329
lat_deg = 37.42692
alt_m = 32.0

print("Ground Station (Stanford):")
print(f"Longitude: {lon_deg:.5f}° = {np.radians(lon_deg):.6f} rad")
print(f"Latitude:  {lat_deg:.5f}° = {np.radians(lat_deg):.6f} rad")
print(f"Altitude:  {alt_m:.1f} m\n")
# Longitude: -122.17329° = -2.132605 rad
# Latitude:  37.42692° = 0.653131 rad
# Altitude:  32.0 m

# Convert ground station to ECEF
geodetic_station = np.array([lon_deg, lat_deg, alt_m])
station_ecef = bh.position_geodetic_to_ecef(geodetic_station, bh.AngleFormat.DEGREES)

print("Ground Station ECEF:")
print(f"x = {station_ecef[0]:.3f} m")
print(f"y = {station_ecef[1]:.3f} m")
print(f"z = {station_ecef[2]:.3f} m\n")
# x = -2700691.122 m
# y = -4292566.016 m
# z = 3855395.780 m

# Define satellite in sun-synchronous orbit at 500 km altitude
# SSO orbit passes over Stanford at approximately 10:30 AM local time
# Orbital elements: [a, e, i, RAAN, omega, M]
oe = np.array([bh.R_EARTH + 500e3, 0.001, 97.8, 240.0, 0.0, 90.0])

# Define epoch when satellite passes near Stanford (Jan 1, 2024, 17:05 UTC)
epoch = bh.Epoch.from_datetime(2024, 1, 1, 17, 5, 0.0, 0.0, bh.TimeSystem.UTC)

# Convert orbital elements to ECI state
sat_state_eci = bh.state_koe_to_eci(oe, bh.AngleFormat.DEGREES)

# Convert ECI state to ECEF at the given epoch
sat_state_ecef = bh.state_eci_to_ecef(epoch, sat_state_eci)
sat_ecef = sat_state_ecef[0:3]  # Extract position only

year, month, day, hour, minute, second, ns = epoch.to_datetime()
print(f"Epoch: {year}-{month:02d}-{day:02d} {hour:02d}:{minute:02d}:{second:06.3f} UTC")
print("Satellite ECEF:")
print(f"x = {sat_ecef[0]:.3f} m")
print(f"y = {sat_ecef[1]:.3f} m")
print(f"z = {sat_ecef[2]:.3f} m\n")

# Convert satellite position to SEZ coordinates relative to ground station
sez = bh.relative_position_ecef_to_sez(
    station_ecef, sat_ecef, bh.EllipsoidalConversionType.GEODETIC
)

print("Satellite position in SEZ frame (relative to Stanford):")
print(f"South:  {sez[0] / 1000:.3f} km")
print(f"East:   {sez[1] / 1000:.3f} km")
print(f"Zenith: {sez[2] / 1000:.3f} km")
print(f"Range:  {np.linalg.norm(sez) / 1000:.3f} km")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define ground station location in geodetic coordinates
    // Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
    let lon_deg = -122.17329_f64;
    let lat_deg = 37.42692_f64;
    let alt_m = 32.0;

    println!("Ground Station (Stanford):");
    println!("Longitude: {:.5}° = {:.6} rad", lon_deg, lon_deg.to_radians());
    println!("Latitude:  {:.5}° = {:.6} rad", lat_deg, lat_deg.to_radians());
    println!("Altitude:  {:.1} m\n", alt_m);
    // Expected output:
    // Longitude: -122.17329° = -2.132605 rad
    // Latitude:  37.42692° = 0.653131 rad
    // Altitude:  32.0 m

    // Convert ground station to ECEF
    let geodetic_station = na::Vector3::new(lon_deg, lat_deg, alt_m);
    let station_ecef = bh::position_geodetic_to_ecef(geodetic_station, bh::AngleFormat::Degrees).unwrap();

    println!("Ground Station ECEF:");
    println!("x = {:.3} m", station_ecef[0]);
    println!("y = {:.3} m", station_ecef[1]);
    println!("z = {:.3} m\n", station_ecef[2]);
    // Expected output:
    // x = -2700691.122 m
    // y = -4292566.016 m
    // z = 3855395.780 m

    // Define satellite in sun-synchronous orbit at 500 km altitude
    // SSO orbit passes over Stanford at approximately 10:30 AM local time
    // Orbital elements: [a, e, i, RAAN, omega, M]
    let oe = na::SVector::<f64, 6>::new(
        bh::R_EARTH + 500e3,
        0.001,
        97.8_f64,
        240.0_f64,
        0.0_f64,
        90.0_f64,
    );

    // Define epoch when satellite passes near Stanford (Jan 1, 2024, 17:05 UTC)
    let epoch = bh::Epoch::from_datetime(2024, 1, 1, 17, 5, 0.0, 0.0, bh::TimeSystem::UTC);

    // Convert orbital elements to ECI state
    let sat_state_eci = bh::state_koe_to_eci(oe, bh::AngleFormat::Degrees);

    // Convert ECI state to ECEF at the given epoch
    let sat_state_ecef = bh::state_eci_to_ecef(epoch, sat_state_eci);
    let sat_ecef = na::Vector3::new(sat_state_ecef[0], sat_state_ecef[1], sat_state_ecef[2]);

    let (year, month, day, hour, minute, second, _ns) = epoch.to_datetime();
    println!("Epoch: {}-{:02}-{:02} {:02}:{:02}:{:06.3} UTC", year, month, day, hour, minute, second);
    println!("Satellite ECEF:");
    println!("x = {:.3} m", sat_ecef[0]);
    println!("y = {:.3} m", sat_ecef[1]);
    println!("z = {:.3} m\n", sat_ecef[2]);

    // Convert satellite position to SEZ coordinates relative to ground station
    let sez = bh::relative_position_ecef_to_sez(
        station_ecef,
        sat_ecef,
        bh::EllipsoidalConversionType::Geodetic,
    );

    println!("Satellite position in SEZ frame (relative to Stanford):");
    println!("South:  {:.3} km", sez[0] / 1000.0);
    println!("East:   {:.3} km", sez[1] / 1000.0);
    println!("Zenith: {:.3} km", sez[2] / 1000.0);
    println!("Range:  {:.3} km", sez.norm() / 1000.0);
}

Converting SEZ to ECEF

The reverse transformation converts a relative SEZ position back to an absolute ECEF position:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define ground station location in geodetic coordinates
# Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
lon_deg = -122.17329
lat_deg = 37.42692
alt_m = 32.0

print("Ground Station (Stanford):")
print(f"Longitude: {lon_deg:.5f}° = {np.radians(lon_deg):.6f} rad")
print(f"Latitude:  {lat_deg:.5f}° = {np.radians(lat_deg):.6f} rad")
print(f"Altitude:  {alt_m:.1f} m\n")
# Longitude: -122.17329° = -2.132605 rad
# Latitude:  37.42692° = 0.653131 rad
# Altitude:  32.0 m

# Convert ground station to ECEF
geodetic_station = np.array([lon_deg, lat_deg, alt_m])
station_ecef = bh.position_geodetic_to_ecef(geodetic_station, bh.AngleFormat.DEGREES)

print("Ground Station ECEF:")
print(f"x = {station_ecef[0]:.3f} m")
print(f"y = {station_ecef[1]:.3f} m")
print(f"z = {station_ecef[2]:.3f} m\n")
# x = -2700691.122 m
# y = -4292566.016 m
# z = 3855395.780 m

# Define relative position in SEZ coordinates
# Example: 30 km South, 50 km East, 100 km Up from station
sez = np.array([30e3, 50e3, 100e3])

print("Relative position in SEZ frame:")
print(f"South:  {sez[0] / 1000:.1f} km")
print(f"East:   {sez[1] / 1000:.1f} km")
print(f"Zenith: {sez[2] / 1000:.1f} km\n")
# South:  30.0 km
# East:   50.0 km
# Zenith: 100.0 km

# Convert SEZ relative position to absolute ECEF position
target_ecef = bh.relative_position_sez_to_ecef(
    station_ecef, sez, bh.EllipsoidalConversionType.GEODETIC
)

print("Target position in ECEF:")
print(f"x = {target_ecef[0]:.3f} m")
print(f"y = {target_ecef[1]:.3f} m")
print(f"z = {target_ecef[2]:.3f} m")
print(f"Distance from Earth center: {np.linalg.norm(target_ecef) / 1000:.3f} km")

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define ground station location in geodetic coordinates
    // Stanford University: (lon=-122.17329°, lat=37.42692°, alt=32.0m)
    let lon_deg = -122.17329_f64;
    let lat_deg = 37.42692_f64;
    let alt_m = 32.0;

    println!("Ground Station (Stanford):");
    println!("Longitude: {:.5}° = {:.6} rad", lon_deg, lon_deg.to_radians());
    println!("Latitude:  {:.5}° = {:.6} rad", lat_deg, lat_deg.to_radians());
    println!("Altitude:  {:.1} m\n", alt_m);
    // Expected output:
    // Longitude: -122.17329° = -2.132605 rad
    // Latitude:  37.42692° = 0.653131 rad
    // Altitude:  32.0 m

    // Convert ground station to ECEF
    let geodetic_station = na::Vector3::new(lon_deg, lat_deg, alt_m);
    let station_ecef = bh::position_geodetic_to_ecef(geodetic_station, bh::AngleFormat::Degrees).unwrap();

    println!("Ground Station ECEF:");
    println!("x = {:.3} m", station_ecef[0]);
    println!("y = {:.3} m", station_ecef[1]);
    println!("z = {:.3} m\n", station_ecef[2]);
    // Expected output:
    // x = -2700691.122 m
    // y = -4292566.016 m
    // z = 3855395.780 m

    // Define relative position in SEZ coordinates
    // Example: 30 km South, 50 km East, 100 km Up from station
    let sez = na::Vector3::new(30e3, 50e3, 100e3);

    println!("Relative position in SEZ frame:");
    println!("South:  {:.1} km", sez[0] / 1000.0);
    println!("East:   {:.1} km", sez[1] / 1000.0);
    println!("Zenith: {:.1} km\n", sez[2] / 1000.0);
    // Expected output:
    // South:  30.0 km
    // East:   50.0 km
    // Zenith: 100.0 km

    // Convert SEZ relative position to absolute ECEF position
    let target_ecef = bh::relative_position_sez_to_ecef(
        station_ecef,
        sez,
        bh::EllipsoidalConversionType::Geodetic,
    );

    println!("Target position in ECEF:");
    println!("x = {:.3} m", target_ecef[0]);
    println!("y = {:.3} m", target_ecef[1]);
    println!("z = {:.3} m", target_ecef[2]);
    println!("Distance from Earth center: {:.3} km", target_ecef.norm() / 1000.0);
}

Azimuth and Elevation from Topocentric Coordinates

For object tracking, it's often more intuitive to work with azimuth (compass direction) and elevation (angle above the horizon) rather than Cartesian ENZ or SEZ coordinates. Both ENZ and SEZ topocentric systems can be converted to azimuth-elevation-range format.

From ENZ Coordinates

Convert ENZ positions to azimuth (measured clockwise from North), elevation (angle above horizon), and range:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define several relative positions in ENZ coordinates
test_cases = [
    ("Directly overhead", np.array([0.0, 0.0, 100e3])),
    ("North horizon", np.array([0.0, 100e3, 0.0])),
    ("East horizon", np.array([100e3, 0.0, 0.0])),
    ("South horizon", np.array([0.0, -100e3, 0.0])),
    ("West horizon", np.array([-100e3, 0.0, 0.0])),
    ("Northeast at 45° elevation", np.array([50e3, 50e3, 70.7e3])),
]

print("Converting ENZ coordinates to Azimuth-Elevation-Range:\n")

for name, enz in test_cases:
    # Convert to azimuth-elevation-range
    azel = bh.position_enz_to_azel(enz, bh.AngleFormat.DEGREES)

    print(f"{name}:")
    print(
        f"  ENZ:   E={enz[0] / 1000:.1f}km, N={enz[1] / 1000:.1f}km, Z={enz[2] / 1000:.1f}km"
    )
    print(
        f"  Az/El: Az={azel[0]:.1f}°, El={azel[1]:.1f}°, Range={azel[2] / 1000:.1f}km\n"
    )

# Expected outputs:
# Directly overhead: Az=0.0°, El=90.0°, Range=100.0km
# North horizon: Az=0.0°, El=0.0°, Range=100.0km
# East horizon: Az=90.0°, El=0.0°, Range=100.0km
# South horizon: Az=180.0°, El=0.0°, Range=100.0km
# West horizon: Az=270.0°, El=0.0°, Range=100.0km

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define several relative positions in ENZ coordinates
    let test_cases = vec![
        ("Directly overhead", na::Vector3::new(0.0, 0.0, 100e3)),
        ("North horizon", na::Vector3::new(0.0, 100e3, 0.0)),
        ("East horizon", na::Vector3::new(100e3, 0.0, 0.0)),
        ("South horizon", na::Vector3::new(0.0, -100e3, 0.0)),
        ("West horizon", na::Vector3::new(-100e3, 0.0, 0.0)),
        ("Northeast at 45° elevation", na::Vector3::new(50e3, 50e3, 70.7e3)),
    ];

    println!("Converting ENZ coordinates to Azimuth-Elevation-Range:\n");

    for (name, enz) in test_cases {
        // Convert to azimuth-elevation-range
        let azel = bh::position_enz_to_azel(enz, bh::AngleFormat::Degrees);

        println!("{}:", name);
        println!("  ENZ:   E={:.1}km, N={:.1}km, Z={:.1}km",
                 enz[0] / 1000.0, enz[1] / 1000.0, enz[2] / 1000.0);
        println!("  Az/El: Az={:.1}°, El={:.1}°, Range={:.1}km\n",
                 azel[0], azel[1], azel[2] / 1000.0);
    }

    // Expected outputs:
    // Directly overhead: Az=0.0°, El=90.0°, Range=100.0km
    // North horizon: Az=0.0°, El=0.0°, Range=100.0km
    // East horizon: Az=90.0°, El=0.0°, Range=100.0km
    // South horizon: Az=180.0°, El=0.0°, Range=100.0km
    // West horizon: Az=270.0°, El=0.0°, Range=100.0km
}

Info

Azimuth is measured clockwise from North (0° = North, 90° = East, 180° = South, 270° = West). Elevation is the angle above the horizon (0° = horizon, 90° = directly overhead).

From SEZ Coordinates

The same conversion is available from SEZ coordinates:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
import brahe as bh
import numpy as np

bh.initialize_eop()

# Define several relative positions in SEZ coordinates
test_cases = [
    ("Directly overhead", np.array([0.0, 0.0, 100e3])),
    ("North horizon", np.array([-100e3, 0.0, 0.0])),
    ("East horizon", np.array([0.0, 100e3, 0.0])),
    ("South horizon", np.array([100e3, 0.0, 0.0])),
    ("West horizon", np.array([0.0, -100e3, 0.0])),
    ("Northeast at 45° elevation", np.array([-50e3, 50e3, 70.7e3])),
]

print("Converting SEZ coordinates to Azimuth-Elevation-Range:\n")

for name, sez in test_cases:
    # Convert to azimuth-elevation-range
    azel = bh.position_sez_to_azel(sez, bh.AngleFormat.DEGREES)

    print(f"{name}:")
    print(
        f"  SEZ:   S={sez[0] / 1000:.1f}km, E={sez[1] / 1000:.1f}km, Z={sez[2] / 1000:.1f}km"
    )
    print(
        f"  Az/El: Az={azel[0]:.1f}°, El={azel[1]:.1f}°, Range={azel[2] / 1000:.1f}km\n"
    )

# Expected outputs:
# Directly overhead: Az=0.0°, El=90.0°, Range=100.0km
# North horizon: Az=0.0°, El=0.0°, Range=100.0km
# East horizon: Az=90.0°, El=0.0°, Range=100.0km
# South horizon: Az=180.0°, El=0.0°, Range=100.0km
# West horizon: Az=270.0°, El=0.0°, Range=100.0km

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
use brahe as bh;
use nalgebra as na;

fn main() {
    bh::initialize_eop().unwrap();

    // Define several relative positions in SEZ coordinates
    let test_cases = vec![
        ("Directly overhead", na::Vector3::new(0.0, 0.0, 100e3)),
        ("North horizon", na::Vector3::new(-100e3, 0.0, 0.0)),
        ("East horizon", na::Vector3::new(0.0, 100e3, 0.0)),
        ("South horizon", na::Vector3::new(100e3, 0.0, 0.0)),
        ("West horizon", na::Vector3::new(0.0, -100e3, 0.0)),
        ("Northeast at 45° elevation", na::Vector3::new(-50e3, 50e3, 70.7e3)),
    ];

    println!("Converting SEZ coordinates to Azimuth-Elevation-Range:\n");

    for (name, sez) in test_cases {
        // Convert to azimuth-elevation-range
        let azel = bh::position_sez_to_azel(sez, bh::AngleFormat::Degrees);

        println!("{}:", name);
        println!("  SEZ:   S={:.1}km, E={:.1}km, Z={:.1}km",
                 sez[0] / 1000.0, sez[1] / 1000.0, sez[2] / 1000.0);
        println!("  Az/El: Az={:.1}°, El={:.1}°, Range={:.1}km\n",
                 azel[0], azel[1], azel[2] / 1000.0);
    }

    // Expected outputs:
    // Directly overhead: Az=0.0°, El=90.0°, Range=100.0km
    // North horizon: Az=0.0°, El=0.0°, Range=100.0km
    // East horizon: Az=90.0°, El=0.0°, Range=100.0km
    // South horizon: Az=180.0°, El=0.0°, Range=100.0km
    // West horizon: Az=270.0°, El=0.0°, Range=100.0km
}

Info

Both ENZ and SEZ produce identical azimuth-elevation-range results for the same physical position. The choice between them is purely a matter of intermediate representation.

See Also