GNSS Satellite Ephemeris Positions Tutorial

%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")
shape: (31_422, 6)
timestampsystemsatellitexyz
datetime[ms, UTC]stri8f64f64f64
2025-09-22 00:00:00 UTC"S"48-2.8819e7-3.0791e75330.304
2025-09-22 00:00:00 UTC"S"285.1385e64.1866e7-23594.416
2025-09-22 00:00:00 UTC"R"59.2172e6-7.9632e6-2.2404e7
2025-09-22 00:00:00 UTC"R"6-763982.297457-2.4121e7-8.2510e6
2025-09-22 00:00:00 UTC"R"7-7.3418e6-2.3132e77.7931e6
………………
2025-09-23 00:00:00 UTC"E"23-867060.4977732.0556e7-2.1280e7
2025-09-23 00:00:00 UTC"E"21-1.7422e75.1905e6-2.3357e7
2025-09-23 00:00:00 UTC"E"251.7758e7-4.6943e62.3221e7
2025-09-23 00:00:00 UTC"E"29-1.7563e71.6301e71.7372e7
2025-09-23 00:00:00 UTC"E"8-4.7838e6-2.3962e71.6705e7

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:

Note

Learn more: check out our tutorial on retrieving 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")
shape: (7_842, 6)
timestampsatelliteobs_codesnrsystemigs
datetime[ms, UTC]u8strf32strstr
2025-09-22 00:00:00 UTC7"1C"48.0"G""P04100USA"
2025-09-22 00:00:15 UTC7"1C"48.0"G""P04100USA"
2025-09-22 00:00:30 UTC7"1C"47.75"G""P04100USA"
2025-09-22 00:00:45 UTC7"1C"47.75"G""P04100USA"
2025-09-22 00:01:00 UTC7"1C"47.75"G""P04100USA"
………………
2025-09-22 23:58:45 UTC7"1C"48.0"G""P04100USA"
2025-09-22 23:59:00 UTC7"1C"48.25"G""P04100USA"
2025-09-22 23:59:15 UTC7"1C"48.5"G""P04100USA"
2025-09-22 23:59:30 UTC7"1C"48.0"G""P04100USA"
2025-09-22 23:59:45 UTC7"1C"48.25"G""P04100USA"

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")
shape: (8_019, 5)
timestampsystemsatelliteelevationazimuth
datetime[ms, UTC]stri8f32f32
2025-09-22 00:00:00 UTC"G"748.274212298.52771
2025-09-22 00:00:15 UTC"G"748.352226298.636108
2025-09-22 00:00:30 UTC"G"748.43026298.744598
2025-09-22 00:00:45 UTC"G"748.508312298.853241
2025-09-22 00:01:00 UTC"G"748.58638298.962006
……………
2025-09-22 23:59:00 UTC"G"749.233974299.867706
2025-09-22 23:59:15 UTC"G"749.312218299.977661
2025-09-22 23:59:30 UTC"G"749.39048300.087738
2025-09-22 23:59:45 UTC"G"749.468758300.197968
2025-09-23 00:00:00 UTC"G"749.547054300.308319

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")
shape: (7_810, 8)
timestampsatelliteobs_codesnrsystemigselevationazimuth
datetime[ms, UTC]u8strf32strstrf32f32
2025-09-22 00:00:00 UTC7"1C"48.0"G""P04100USA"48.274212298.52771
2025-09-22 00:00:15 UTC7"1C"48.0"G""P04100USA"48.352226298.636108
2025-09-22 00:00:30 UTC7"1C"47.75"G""P04100USA"48.43026298.744598
2025-09-22 00:00:45 UTC7"1C"47.75"G""P04100USA"48.508312298.853241
2025-09-22 00:01:00 UTC7"1C"47.75"G""P04100USA"48.58638298.962006
……………………
2025-09-22 23:58:45 UTC7"1C"48.0"G""P04100USA"49.15575299.757874
2025-09-22 23:59:00 UTC7"1C"48.25"G""P04100USA"49.233974299.867706
2025-09-22 23:59:15 UTC7"1C"48.5"G""P04100USA"49.312218299.977661
2025-09-22 23:59:30 UTC7"1C"48.0"G""P04100USA"49.39048300.087738
2025-09-22 23:59:45 UTC7"1C"48.25"G""P04100USA"49.468758300.197968

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
Contents