Example 13: Visualizing Constellations in Cesium

This tutorial shows how to turn LuPNT trajectory samples into interactive CesiumJS scenes. Visualization is not just presentation: it is a useful debugging tool for frame conventions, body-fixed coordinates, surface locations, and moving link geometry.

We build three scene layers:

  1. Earth GNSS satellites from TLEs in an Earth-fixed view;

  2. lunar LCRNS relay satellites propagated in MOON_CI and rendered in the Moon-fixed MOON_PA frame, with lunar surface stations;

  3. Earth GNSS + LCRNS relay in Earth centered inertial frame

pnt.plot.CesiumScene writes a standalone HTML file containing CZML, Cesium’s time-dynamic JSON format. You provide Cartesian samples in the correct frame, or pass inertial samples to add_trajectory(...) and let LuPNT rotate them to the desired frame before rendering.

No Cesium ion account is required for this tutorial; LuPNT embeds a textured ellipsoid and loads only the CesiumJS library from a public CDN.

1. Setup

[1]:
import os, sys
from pathlib import Path
from datetime import datetime, timezone

import numpy as np

# Use THIS repo's pylupnt/data, not a stale copy on a global PYTHONPATH (~/.zshrc) or kernelspec.
for _b in [Path.cwd(), *Path.cwd().parents]:
    if (_b / "python/pylupnt/__init__.py").exists():
        sys.path.insert(0, str((_b / "python").resolve()))
        if (_b / "data/LuPNT_data/ephemeris").is_dir():
            os.environ["LUPNT_DATA_PATH"] = str((_b / "data/LuPNT_data").resolve())
        REPO_ROOT = _b
        break
else:
    REPO_ROOT = Path.cwd()
import pylupnt as pnt

# Cesium scenes are written under the repo-root output dir (output/python_examples/).
SCENES_DIR = REPO_ROOT / "output" / "python_examples" / "cesium_scenes"
SCENES_DIR.mkdir(parents=True, exist_ok=True)


def find_lupnt_data_dir():
    env = os.environ.get("LUPNT_DATA_PATH")
    if env and Path(env).is_dir():
        return Path(env).resolve()
    for base in [Path.cwd(), *Path.cwd().parents]:
        for cand in (base / "data" / "LuPNT_data", base / "LuPNT_data"):
            if cand.is_dir():
                return cand.resolve()
    raise FileNotFoundError("Could not locate LuPNT_data; set LUPNT_DATA_PATH.")


DATA_DIR = find_lupnt_data_dir()
print(f"LuPNT_data: {DATA_DIR}")

# Host heavy interactive figures on lupnt-doc-assets instead of embedding them.
import sys as _sys
from pathlib import Path as _Path

for _p in [_Path.cwd(), *_Path.cwd().parents]:
    if (_p / "python" / "examples" / "_doc_assets.py").exists():
        _sys.path.insert(0, str(_p / "python" / "examples"))
        break
from _doc_assets import embed_plotly, embed_cesium_scene
[00.00][PyLuPNT] Initializing
LuPNT_data: /Users/keidaiiiyama/Documents/sw_navlab/LuPNT/data/LuPNT_data

2. Earth GNSS constellation (GPS + Galileo)

We load the GPS and Galileo two-line elements archived at 2025-01-01, propagate each satellite with two-body dynamics over 16 hours (more than one full orbit for both systems), and hand the inertial (ECI) states to add_trajectory, which rotates them into the Earth-fixed frame for us.

[2]:
T0_TDB = pnt.gregorian_to_time(2025, 1, 1, 0, 0, 0)  # LuPNT internal time (~TDB)
T0_TAI = pnt.convert_time(T0_TDB, pnt.Time.TDB, pnt.Time.TAI)
EPOCH_EARTH = datetime(2025, 1, 1, tzinfo=timezone.utc)

# pnt.GM_EARTH is SI (m^3/s^2), consistent with the metre states used here, so it can be
# used directly for the Earth two-body propagation.
GM_EARTH_SI = pnt.GM_EARTH  # [m^3/s^2]

DT_E, SPAN_E = 240.0, 16 * 3600  # 4-min steps over 16 h (> 1 orbit for GPS & Galileo)
t_tai = T0_TAI + np.arange(0, SPAN_E + DT_E, DT_E)
t_tdb = pnt.convert_time(t_tai, pnt.Time.TAI, pnt.Time.TDB)


def load_tles(filepath):
    lines = [l.strip() for l in open(filepath) if l.strip()]
    tles, i = [], 0
    while i + 2 < len(lines):
        if lines[i + 1].startswith("1") and lines[i + 2].startswith("2"):
            tles.append(pnt.TLE.from_lines(lines[i], lines[i + 1], lines[i + 2]))
            i += 3
        else:
            i += 1
    return tles


TLE_DIR = DATA_DIR / "tle/2025_01_01"
gps_tles = load_tles(TLE_DIR / "gps_2025_01_01.txt")
gal_tles = load_tles(TLE_DIR / "galileo_2025_01_01.txt")

_dyn_e = pnt.CartesianTwoBodyDynamics(GM_EARTH_SI)
_dyn_e.set_integrator(pnt.IntegratorType.RKF45)


def gnss_eci(tle):
    coe = np.array(pnt.tle_to_classical(tle, T0_TAI))  # classical elements in ECI
    rv0 = pnt.classical_to_cart(coe, GM_EARTH_SI)  # -> ECI state at epoch
    return np.array(_dyn_e.propagate(rv0, t_tai))  # [N, 6] ECI


print(f"GPS: {len(gps_tles)} sats   Galileo: {len(gal_tles)} sats")
GPS: 31 sats   Galileo: 31 sats
[3]:
# Build the Earth scene. add_trajectory() converts ECI -> ECEF (the scene's fixed frame).
earth = pnt.plot.CesiumScene(
    body="EARTH",
    name="Earth GNSS constellation",
    epoch=EPOCH_EARTH,
    multiplier=120,
    texture=True,
)
for i, tle in enumerate(gps_tles):
    earth.add_trajectory(
        f"GPS-{i+1}", t_tdb, gnss_eci(tle), frame_in=pnt.ECI, color=(0, 180, 255)
    )  # cyan
for i, tle in enumerate(gal_tles):
    earth.add_trajectory(
        f"GAL-{i+1}", t_tdb, gnss_eci(tle), frame_in=pnt.ECI, color=(255, 180, 0)
    )  # amber

# Cyan = GPS, amber = Galileo. Press play (bottom-left) or scrub the timeline; drag to
# rotate, scroll to zoom, click a satellite for its name.
embed_cesium_scene(earth, "earth_gnss_constellation")
[doc-asset] cesium/earth_gnss_constellation.html  (2044 kB)  ->  https://stanford-navlab.github.io/lupnt-doc-assets/cesium/earth_gnss_constellation.html

3. Lunar navigation constellation (LCRNS relays)

Five satellites in Elliptical Lunar Frozen Orbits (ELFOs) from the NASA LCRNS Reference Constellation 3.1 (Moon-centered ICRF states at 2027-03-01). We propagate them for 31 hours (one full ELFO period ≈ 30 h) with multibody dynamics — Moon 8×8 gravity + Earth + Sun third bodies — and let add_trajectory rotate the inertial MOON_CI states into the Moon-fixed MOON_PA frame. The displayed orbit path is the finite propagated arc, not a mathematically closed conic, so the top of the ELFO may look clipped where the first and last samples do not wrap into one another.

[4]:
M0_TDB = pnt.gregorian_to_time(2027, 3, 1, 0, 0, 0)
EPOCH_MOON = datetime(2027, 3, 1, tzinfo=timezone.utc)

# LCRNS Reference Constellation 3.1: Moon-centered inertial (MOON_CI) states, km & km/s.
lcrns_states = {
    "SV-1": ([-198.931445, 386.023132, 3458.355412], [-1.539867, 0.022243, -0.091059]),
    "SV-2": (
        [1089.980598, -2260.918271, -18984.562823],
        [0.280301, -0.004049, 0.016575],
    ),
    "SV-3": ([9187.665852, -260.470721, -8543.222759], [0.273629, 0.417588, -0.313828]),
    "SV-4": (
        [6484.019002, 14329.883921, -7966.345594],
        [-0.258433, 0.033951, 0.235264],
    ),
    "SV-5": (
        [-5074.242314, 14473.022751, -8638.530033],
        [-0.229931, -0.026926, -0.264381],
    ),
}

DT_M, SPAN_M = 600.0, 31 * 3600  # 10-min steps over 31 h (~one full ELFO period)
t_moon_tdb = M0_TDB + np.arange(0, SPAN_M + DT_M, DT_M)

dyn_m = pnt.NBodyDynamics()
dyn_m.set_frame(pnt.MOON_CI)
dyn_m.set_units(pnt.SI_UNITS)
dyn_m.add_body(pnt.create_body(pnt.MOON, 8, 8))
dyn_m.add_body(pnt.create_body(pnt.EARTH))
dyn_m.add_body(pnt.create_body(pnt.SUN))
dyn_m.set_integrator(pnt.IntegratorType.RKF45)


def relay_mci(r_km, v_kmps):
    rv0 = np.hstack([np.array(r_km) * 1e3, np.array(v_kmps) * 1e3])
    return np.array(dyn_m.propagate(rv0, t_moon_tdb))  # [N, 6] MOON_CI


relay_states = {name: relay_mci(r, v) for name, (r, v) in lcrns_states.items()}
for name, rv in relay_states.items():
    alt = np.linalg.norm(rv[:, :3], axis=1) / 1e3 - pnt.R_MOON / 1e3
    print(
        f"{name}: periapsis {alt.min():7.0f} km   apoapsis {alt.max():7.0f} km  altitude"
    )
SV-1: periapsis    1728 km   apoapsis   17406 km  altitude
SV-2: periapsis    1704 km   apoapsis   17433 km  altitude
SV-3: periapsis    1890 km   apoapsis   17325 km  altitude
SV-4: periapsis    1802 km   apoapsis   17289 km  altitude
SV-5: periapsis    1723 km   apoapsis   17379 km  altitude

4. Lunar surface stations

Ground assets are added with add_station(lat=..., lon=...) (degrees), placed on a sphere of the Moon’s radius in the same MOON_PA frame as the relay orbits — so stations and satellites share one consistent scene. These sites cluster around the lunar south pole (the Artemis region of interest) plus a near/far-side pair. In the render cell below, always_visible=True keeps every HQ marker/label visible through the Moon, which is useful for constellation-level link geometry even though the far-side sites would be physically occulted from the current camera view.

[5]:
surface_sites = {  # name: (lat_deg, lon_deg)
    "Shackleton HQ": (-89.9, 0.0),
    "Connecting Ridge HQ": (-89.45, 222.8),
    "de Gerlache HQ": (-88.5, 290.0),
    "Nearside HQ": (0.0, 0.0),
    "Farside HQ": (10.0, 180.0),
}

# Spread labels apart; several far-side/south-pole HQs project into the same screen area.
site_label_offsets = {
    "Shackleton HQ": (-82, -24),
    "Connecting Ridge HQ": (-118, 16),
    "de Gerlache HQ": (18, -42),
    "Nearside HQ": (14, 2),
    "Farside HQ": (20, -26),
}

5. Render the lunar scene (relays + surface stations)

The render uses dashed ELFO trajectories and dashed dynamic relay-to-HQ links. moon.show() writes output/python_examples/cesium_scenes/lunar_navigation_constellation.html and displays it in an iframe when the notebook front-end allows CDN JavaScript. For the most reliable view with the embedded Moon texture, serve the repo root and open the saved HTML directly:

cd /Users/keidaiiiyama/Documents/sw_navlab/LuPNT
python -m http.server 8765

Then open http://127.0.0.1:8765/output/python_examples/cesium_scenes/lunar_navigation_constellation.html

[6]:
moon = pnt.plot.CesiumScene(
    body="MOON",
    name="Lunar navigation constellation",
    epoch=EPOCH_MOON,
    multiplier=400,
    texture=True,
)

relay_colors = [
    (255, 90, 90),
    (90, 255, 120),
    (255, 210, 0),
    (180, 120, 255),
    (0, 220, 255),
]
for (name, rv), col in zip(relay_states.items(), relay_colors):
    moon.add_trajectory(
        name,
        t_moon_tdb,
        rv,
        frame_in=pnt.MOON_CI,
        color=col,
        label=True,
        dashed=False,
        width=0.7,
    )

for name, (lat, lon) in surface_sites.items():
    moon.add_station(
        name,
        lat=lat,
        lon=lon,
        always_visible=True,
        label_offset=site_label_offsets[name],
    )

for relay_name in relay_states:
    for site_name in surface_sites:
        moon.add_link(
            relay_name, site_name, color=(255, 255, 255, 95), width=0.45, dashed=False
        )

# Solid ELFO relay trajectories + white HQs + solid relay-to-HQ links, all in MOON_PA.
# Press play to watch the relays sweep over the south-pole sites; click any entity for its name.
# If the inline iframe is blocked, run `python -m http.server 8765` from the repo root and open
# http://127.0.0.1:8765/docs/_static/cesium_scenes/lunar_navigation_constellation.html in a browser.
embed_cesium_scene(moon, "lunar_navigation_constellation")
[doc-asset] cesium/lunar_navigation_constellation.html  (948 kB)  ->  https://stanford-navlab.github.io/lupnt-doc-assets/cesium/lunar_navigation_constellation.html

6. Combined Earth GPS + lunar constellation scene

This scene uses one Earth-centered inertial frame so Earth, GPS/Galileo, the moving Moon, and the LCRNS relay orbits can share a single Cesium clock. The Cesium camera is kept in an inertial view so the Earth texture rotates under the inertial trajectories, and the Moon ellipsoid uses LuPNT MOON_PA -> ECI orientation quaternions so its texture follows the Moon’s body-fixed attitude. The Earth/Moon geometry is true scale: GPS and Galileo stay close to Earth, while the Moon and its relays appear much farther away. The GPS/Galileo TLE constellation is reused as a visual Earth-orbit layer on the combined timeline; the lunar constellation uses the 2027 LCRNS epoch above.

Similarly, you can render this by opening http://127.0.0.1:8765/output/python_examples/cesium_scenes/earth_moon_gps_lcrns.html

[7]:
COMBO_DT, COMBO_SPAN = 600.0, SPAN_M
t_combo_tdb = M0_TDB + np.arange(0, COMBO_SPAN + COMBO_DT, COMBO_DT)
t_combo_tai = pnt.convert_time(t_combo_tdb, pnt.Time.TDB, pnt.Time.TAI)


def gnss_eci_combo(tle):
    coe = np.array(pnt.tle_to_classical(tle, T0_TAI))
    rv0 = pnt.classical_to_cart(coe, GM_EARTH_SI)
    return np.array(_dyn_e.propagate(rv0, t_combo_tai))


combined = pnt.plot.CesiumScene(
    body="EARTH",
    name="Earth GPS + lunar navigation constellation",
    epoch=EPOCH_MOON,
    multiplier=800,
    texture=True,
    inertial_view=True,
)

for i, tle in enumerate(gps_tles):
    rv_eci = gnss_eci_combo(tle)
    combined.add_satellite(
        f"GPS-{i+1}",
        rv_eci[:, :3],
        offsets_s=t_combo_tdb - t_combo_tdb[0],
        color=(0, 180, 255),
        dashed=False,
        width=0.5,
        pixel_size=5,
        reference_frame="INERTIAL",
    )
for i, tle in enumerate(gal_tles):
    rv_eci = gnss_eci_combo(tle)
    combined.add_satellite(
        f"GAL-{i+1}",
        rv_eci[:, :3],
        offsets_s=t_combo_tdb - t_combo_tdb[0],
        color=(255, 180, 0),
        dashed=False,
        width=0.5,
        pixel_size=5,
        reference_frame="INERTIAL",
    )

moon_eci = np.asarray(pnt.get_body_pos_vel(t_combo_tdb, pnt.EARTH, pnt.MOON, pnt.ECI))[
    :, :3
]
moon_pa_to_eci = pnt.plot.frame_orientation_quaternions(
    t_combo_tdb, pnt.MOON_PA, pnt.ECI
)
combined.add_body_trajectory(
    "Moon",
    moon_eci,
    offsets_s=t_combo_tdb - t_combo_tdb[0],
    radius=pnt.R_MOON,
    color=(170, 170, 180, 230),
    path_color=(210, 210, 230, 120),
    path_width=0.6,
    reference_frame="INERTIAL",
    texture="MOON",
    orientation_quat=moon_pa_to_eci,
)

for (name, rv), col in zip(relay_states.items(), relay_colors):
    relay_eci = np.asarray(pnt.convert_frame(t_combo_tdb, rv, pnt.MOON_CI, pnt.ECI))[
        :, :3
    ]
    combined.add_satellite(
        f"LCRNS {name}",
        relay_eci,
        offsets_s=t_combo_tdb - t_combo_tdb[0],
        color=col,
        label=True,
        dashed=False,
        width=0.7,
        pixel_size=7,
        reference_frame="INERTIAL",
    )

# True-scale Earth-Moon view: scroll toward Earth or Moon to inspect either constellation.
embed_cesium_scene(combined, "earth_moon_gps_lcrns", height=620)
[doc-asset] cesium/earth_moon_gps_lcrns.html  (2481 kB)  ->  https://stanford-navlab.github.io/lupnt-doc-assets/cesium/earth_moon_gps_lcrns.html

7. Summary

This notebook created interactive Earth, Moon, and combined Earth-Moon scenes from LuPNT trajectory data.

Key points:

  • CesiumScene expects body-fixed Cartesian samples for rendering.

  • add_trajectory(...) can rotate inertial LuPNT states into the render frame.

  • add_station(lat, lon) places surface assets consistently with the selected body frame.

  • The generated output/python_examples/cesium_scenes/*.html files are standalone and can be opened or hosted locally.

These scenes are useful companions to the ODTS examples: plot truth and estimated trajectories together, add dynamic relay-to-user links, or visualize when geometry becomes poor.