%pip install altair plotly polars vegafusion vl-convert-python

GNSS Satellite Ephemeris Positions Tutorial¶

This tutorial demonstrates how to retrieve GNSS satellite positions using the EarthScope SDK.

import datetime as dt
import math

import altair as alt
import numpy as np
import plotly.colors as pcolors
import plotly.graph_objects as go
import plotly.io as pio
import polars as pl

from earthscope_sdk import AsyncEarthScopeClient
from earthscope_sdk.client.data_access.models import (
    FloatFilter,
    GeodeticCoordinate,
    SatelliteSystem,
)

alt.data_transformers.enable("vegafusion")  # Enable rust backend for altair
pio.renderers.default = "notebook"  # required for plotly to render in jupyterbook

es = AsyncEarthScopeClient()

Retrieving Satellite Positions¶

You can easily retrieve satellite positions in an Earth-Centered, Earth-Fixed reference frame using the SDK as demonstrated in the following cell.

Choose your desired time range, constellation(s), satellite(s) (optional), and desired sample period.

start_datetime = dt.datetime(2025, 9, 22)
end_datetime = dt.datetime(2025, 9, 23)
system = ["G", "R", "E", "C", "S"]

table = await es.data.gnss_ephemeris_positions(
    start_datetime=start_datetime,
    end_datetime=end_datetime,
    field=["x", "y", "z"],
    system=system,
    sample_interval=dt.timedelta(minutes=5),
).fetch()
df = pl.from_arrow(table)
df.sort("timestamp")

Plotting Satellite Positions¶

The following cell defines a function that uses plotly to plot the satellite orbits in an Earth-Centered, Earth-Fixed (ECEF) reference frame.

def plot_satellite_positions(df: pl.DataFrame, system: SatelliteSystem):
    """
    Plot satellite positions in ECEF reference frame with an earth-sized ellipsoid.
    """
    # Group coordinates by svid
    coords = df.filter(pl.col("system") == system).to_dicts()
    coords_by_svid = {}
    for coord in coords:
        svid = coord["satellite"]
        if svid not in coords_by_svid:
            coords_by_svid[svid] = {"x": [], "y": [], "z": [], "timestamp": []}
        coords_by_svid[svid]["x"].append(coord["x"])
        coords_by_svid[svid]["y"].append(coord["y"])
        coords_by_svid[svid]["z"].append(coord["z"])
        coords_by_svid[svid]["timestamp"].append(coord["timestamp"])

    # Earth's radii in meters (WGS84)
    equatorial_radius = 6378137.0
    polar_radius = 6356752.3

    # Create the Earth's ellipsoid wireframe
    u = np.linspace(0, 2 * np.pi, 100)
    v = np.linspace(0, np.pi, 100)
    x_earth = equatorial_radius * np.outer(np.cos(u), np.sin(v))
    y_earth = equatorial_radius * np.outer(np.sin(u), np.sin(v))
    z_earth = polar_radius * np.outer(np.ones(np.size(u)), np.cos(v))

    earth_wireframe = go.Surface(
        x=x_earth,
        y=y_earth,
        z=z_earth,
        colorscale="Blues",
        showscale=False,
        opacity=0.2,
        name="Earth",
    )

    # Create a 3D scatter plot for each satellite trajectory
    satellite_traces = []
    colors = pcolors.qualitative.Plotly
    for i, (svid, coords_list) in enumerate(coords_by_svid.items()):
        satellite_traces.append(
            go.Scatter3d(
                x=coords_list["x"],
                y=coords_list["y"],
                z=coords_list["z"],
                mode="lines",
                text=coords_list["timestamp"],
                hoverinfo="all",
                line=dict(color=colors[i % len(colors)], width=2),
                name=f"SVID {svid}",
            )
        )

    fig = go.Figure(data=[earth_wireframe] + satellite_traces)

    # Update the layout for a cleaner look
    fig.update_layout(
        title="Satellite XYZ Coordinates (System: {})".format(system),
        scene=dict(
            xaxis_title="X Coordinate",
            yaxis_title="Y Coordinate",
            zaxis_title="Z Coordinate",
            aspectmode="data",
        ),
        margin=dict(r=20, b=10, l=10, t=40),
    )

    return fig

fig = plot_satellite_positions(df, "C")
fig.show(height=1500)

Elevation and Azimuth¶

You can also provide a reference point to compute the elevation and azimuth to each satellite at a given sample interval.

Optionally, filters may be provided for elevation and azimuth.

In the following examples, we use elevation and azimuth to visualize the SNR in multiple ways.

Joining Elevation and Azimuth with Observations¶

With satellite elevation and azimuth, we can analyze observations in a spatial context.

First, let’s query the desired GNSS observations:

system = "G"
satellite = [5, 7, 11, 13]
station_name = "P041"
station_location = GeodeticCoordinate(
    latitude=39.7392358,
    longitude=-104.990251,
    height=1000,
)

obs_df = await es.data.gnss_observations(
    start_datetime=start_datetime,
    end_datetime=end_datetime,
    system=system,
    satellite=satellite,
    station_name=station_name,
    session_name="A",
    obs_code="1C",
    field="snr",
).fetch()
obs_df = pl.from_arrow(obs_df)
obs_df.sort("timestamp")

Next, we retrieve satellite elevation and azimuth for the same system and satellites.

table = await es.data.gnss_ephemeris_positions(
    start_datetime=start_datetime,
    end_datetime=end_datetime,
    field=["elevation", "azimuth"],
    system=system,
    satellite=satellite,
    reference_point=station_location,
    sample_interval=dt.timedelta(seconds=15),
    elevation_filter=FloatFilter(min=0),
).fetch()
azel_df = pl.from_arrow(table)
azel_df.sort("timestamp")

Finally, we join the two together and visualize it in different ways.

join_df = obs_df.join(
    azel_df,
    on=["timestamp", "satellite", "system"],
    how="inner",
)
join_df.sort("timestamp")

Elevation vs. SNR¶

join_df.filter(pl.col("elevation") < 30).with_columns(
    pl.col("satellite").cast(pl.Utf8)
).plot.point(x="elevation", y="snr", color="satellite").properties(
    title=f"{station_name}: Mean SNR by Elevation and Azimuth",
    width=650,
)

Polar Heatmap¶

# The main heatmap of SNR by elevation and azimuth
heatmap = (
    alt.Chart(join_df)
    .mark_arc(tooltip=True)
    .encode(
        alt.Theta(
            "azimuth:Q",
            bin=alt.Bin(maxbins=180),
            scale=alt.Scale(domain=[0, 360]),
        ),
        alt.Radius(
            "elevation:Q",
            bin=alt.Bin(step=2),
            scale=alt.Scale(domain=[90, 0]),  # Zenith at center
        ),
        alt.Color("mean(snr):Q", title="Mean SNR", scale=alt.Scale(scheme="viridis")),
    )
)

# --- Axis Layers ---

# Radial axis (for elevation)
elevation_rings_df = pl.DataFrame({"elevation": [75, 60, 45, 30, 15, 0]})
radius_scale = alt.Scale(domain=[90, 0])

# Concentric circles as gridlines
elevation_rings = (
    alt.Chart(elevation_rings_df)
    .mark_arc(fill=None, stroke="lightgrey", strokeDash=[2, 2])
    .encode(
        theta=alt.value(2 * math.pi),
        radius=alt.Radius("elevation:Q", scale=radius_scale).stack(False),
        detail="elevation:Q",  # Draw a separate arc for each elevation value
    )
)

# Text labels for elevation
elevation_labels = (
    alt.Chart(elevation_rings_df)
    .mark_text(color="grey", align="right", size=14)
    .encode(
        text=alt.Text("elevation:Q", format=".0f"),
        theta=alt.value(math.pi / 12),  # Place labels at 15 degrees
        radius=alt.Radius("elevation:Q", scale=radius_scale),
        detail="elevation:Q",  # Ensure a label for each elevation value
    )
)

# Azimuth axis labels
azimuth_labels_df = pl.DataFrame(
    {
        "azimuth": [0, 30, 60, 90, 120, 150, 180, 210, 240, 270, 300, 330],
        "label": [
            "N",
            "30°",
            "60°",
            "E",
            "120°",
            "150°",
            "S",
            "210°",
            "240°",
            "W",
            "300°",
            "330°",
        ],
    }
)

azimuth_labels = (
    alt.Chart(azimuth_labels_df)
    .mark_text(
        radius=265,  # Place labels just outside the plot
        color="grey",
        size=12,
    )
    .encode(text="label:N", theta=alt.Theta("azimuth:Q"))
)

# Layer all the charts
chart = (
    alt.layer(heatmap, elevation_rings, elevation_labels, azimuth_labels)
    .properties(
        title=f"{station_name}: Mean SNR by Elevation and Azimuth",
        width=500,
        height=500,
    )
    .configure_title(fontSize=16)
)

chart