Skip to content

Predicting Ground Contacts

In this example we'll predict ground contacts between the NISAR radar satellite and the NASA Near Earth Network ground stations using Brahe. We'll download the latest TLE data for the satellite from CelesTrak, load the NASA Near Earth Network ground station data, and compute the ground contacts between the satellite and ground stations over a 7-day period. We'll then analyze the statistics of the contact duration and number of contacts per ground station.

Setup

First, we'll import the necessary libraries, initialize Earth orientation parameters, download the TLE for NISAR (NORAD ID 65053) from CelesTrak, and load the NASA Near Earth Network ground station network.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
import time
import csv
import os
import pathlib
import sys
import brahe as bh
import numpy as np
import plotly.graph_objects as go

bh.initialize_eop()

We download the NISAR TLE directly by NORAD ID and load all NASA NEN ground stations:

1
2
nisar = bh.datasets.celestrak.get_tle_by_id_as_propagator(65053, 60.0)
nisar = nisar.with_name("NISAR")

1
nen_stations = bh.datasets.groundstations.load("nasa nen")

We then propagate NISAR for 3 orbital periods to prepare for ground track visualization:

1
2
orbital_period = bh.orbital_period(nisar.semi_major_axis)
nisar.propagate_to(nisar.epoch + 3 * orbital_period)

Ground Track Visualization

Next we'll visualize the ground track and communication cones for NISAR over a 3-orbit period. The communication cones show the coverage area of each ground station based on a 5° minimum elevation angle:

 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
fig_groundtrack = bh.plot_groundtrack(
    trajectories=[
        {
            "trajectory": nisar.trajectory,
            "color": "red",
            "line_width": 2.0,
            "track_length": 3,
            "track_units": "orbits",
        }
    ],
    ground_stations=[
        {
            "stations": nen_stations,
            "color": "blue",
            "alpha": 0.15,
            "point_size": 5.0,
        }
    ],
    gs_cone_altitude=nisar.semi_major_axis - bh.R_EARTH,
    gs_min_elevation=5.0,
    basemap="natural_earth",
    show_borders=True,
    show_coastlines=True,
    show_legend=False,
    backend="plotly",
)

The resulting plot shows NISAR's ground track in red and the NASA Near Earth Network stations with their communication cones in blue:

Compute Ground Contacts

We'll compute the ground contacts between NISAR and the NASA Near Earth Network ground stations over a 7-day period using Brahe's access computation tools. We use an elevation constraint of 5° minimum elevation:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
nisar.reset()

epoch_start = nisar.epoch
epoch_end = epoch_start + 7 * 86400.0  # 7 days in seconds

# Propagate for full 7-day period
nisar.propagate_to(epoch_end)

# Compute access windows with 5 degree minimum elevation
constraint = bh.ElevationConstraint(min_elevation_deg=5.0)
windows = bh.location_accesses(
    nen_stations, [nisar], epoch_start, epoch_end, constraint
)

Below is the table of the first 20 contact windows. Click on any column header to sort:

Station Start Time (UTC) End Time (UTC) Duration (min) Max Elevation (deg)
Svalbard 2026-01-03 07:09:50 2026-01-03 07:19:11 9.3 16.2
Kiruna 2026-01-03 07:10:43 2026-01-03 07:13:31 2.8 5.7
Gilmore Creek 2026-01-03 07:16:23 2026-01-03 07:26:42 10.3 23.3
Alaska Satellite Facility 2026-01-03 07:16:26 2026-01-03 07:26:45 10.3 23.3
North Pole 2026-01-03 07:16:27 2026-01-03 07:26:43 10.3 23
Awarua 2026-01-03 07:51:02 2026-01-03 07:58:47 7.7 11.4
McMurdo 2026-01-03 07:58:36 2026-01-03 08:09:20 10.7 23.2
Troll 2026-01-03 08:05:15 2026-01-03 08:17:31 12.3 71.1
Wallops 2026-01-03 08:43:14 2026-01-03 08:47:26 4.2 6.7
Svalbard 2026-01-03 08:51:55 2026-01-03 08:58:58 7.1 10
Gilmore Creek 2026-01-03 08:54:39 2026-01-03 09:03:00 8.4 13.8
Alaska Satellite Facility 2026-01-03 08:54:43 2026-01-03 09:03:02 8.3 13.7
North Pole 2026-01-03 08:54:43 2026-01-03 09:02:59 8.3 13.5
Dongara 2026-01-03 09:25:00 2026-01-03 09:36:38 11.6 42.1
McMurdo 2026-01-03 09:40:12 2026-01-03 09:48:53 8.7 13.5
Troll 2026-01-03 09:43:53 2026-01-03 09:55:31 11.6 38.8
Santiago 2026-01-03 09:58:29 2026-01-03 10:08:27 10 19.6
Merrit Island 2026-01-03 10:15:30 2026-01-03 10:26:20 10.8 28.6
Ponce de Leon 2026-01-03 10:15:40 2026-01-03 10:26:29 10.8 28.4
Wallops 2026-01-03 10:17:03 2026-01-03 10:29:01 12 65.5

Analyze Contact Statistics

Finally, we'll analyze the contact statistics, including the average daily contacts per station and distribution of contact durations.

We group the contact windows by station and compute the average daily contacts:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
# Group contacts by station
station_contacts = {}
for window in windows:
    station_name = window.location_name
    if station_name not in station_contacts:
        station_contacts[station_name] = []
    station_contacts[station_name].append(window)

# Calculate average daily contacts per station
days = 7.0
station_daily_avg = {}
for station, contacts in station_contacts.items():
    avg_per_day = len(contacts) / days
    station_daily_avg[station] = avg_per_day

# Sort by average daily contacts
sorted_stations = sorted(station_daily_avg.items(), key=lambda x: x[1], reverse=True)

Then we create two visualizations: a bar chart of average daily contacts per station, and a histogram of contact duration distribution:

 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
# Figure 1: Daily contacts per station (bar chart)
stations_list = [s[0] for s in sorted_stations]
daily_avgs = [s[1] for s in sorted_stations]

fig_daily = go.Figure(
    data=[
        go.Bar(
            x=stations_list,
            y=daily_avgs,
            marker_color="steelblue",
            hovertemplate="<b>%{x}</b><br>%{y:.1f} contacts/day<extra></extra>",
        )
    ]
)
fig_daily.update_layout(
    title="NISAR Average Daily Contacts by Station (7-day period)",
    xaxis_title="Ground Station",
    yaxis_title="Average Daily Contacts",
    xaxis_tickangle=-45,
    height=700,
    margin=dict(l=60, r=40, t=80, b=120),
)

# Figure 2: Contact duration distribution (histogram)
durations = [w.duration / 60.0 for w in windows]  # Convert to minutes

mean_duration = np.mean(durations)
median_duration = np.median(durations)
max_duration = np.max(durations)

fig_duration = go.Figure(
    data=[
        go.Histogram(
            x=durations,
            nbinsx=30,
            marker_color="coral",
            marker_line_color="black",
            marker_line_width=1,
            hovertemplate="Duration: %{x:.1f} min<br>Count: %{y}<extra></extra>",
        )
    ]
)
fig_duration.update_layout(
    title="NISAR Contact Duration Distribution",
    xaxis_title="Contact Duration (minutes)",
    yaxis_title="Frequency",
    height=700,
    margin=dict(l=60, r=40, t=80, b=60),
    annotations=[
        dict(
            text=f"Mean: {mean_duration:.1f} min<br>Median: {median_duration:.1f} min<br>Max: {max_duration:.1f} min",
            xref="paper",
            yref="paper",
            x=0.05,
            y=0.97,
            xanchor="left",
            yanchor="top",
            showarrow=False,
            bordercolor="grey",
            borderwidth=1,
            borderpad=8,
        )
    ],
)

The daily contacts bar chart shows which stations have the most frequent visibility to NISAR:

The duration histogram shows the distribution of contact lengths, with statistics overlay:

Full Code Example

ground_contacts.py
  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
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
import time
import csv
import os
import pathlib
import sys
import brahe as bh
import numpy as np
import plotly.graph_objects as go

bh.initialize_eop()

# Configuration for output files
SCRIPT_NAME = pathlib.Path(__file__).stem
OUTDIR = pathlib.Path(os.getenv("BRAHE_FIGURE_OUTPUT_DIR", "./docs/figures/"))
os.makedirs(OUTDIR, exist_ok=True)

# Download TLE data for NISAR from CelesTrak
# NISAR (NASA-ISRO SAR) is a joint Earth observation satellite
# NORAD ID: 65053
print("Downloading NISAR TLE from CelesTrak...")
start_time = time.time()
nisar = bh.datasets.celestrak.get_tle_by_id_as_propagator(65053, 60.0)
nisar = nisar.with_name("NISAR")
elapsed = time.time() - start_time
print(f"Downloaded NISAR TLE in {elapsed:.2f} seconds.")
print(f"Epoch: {nisar.epoch}")
print(f"Semi-major axis: {nisar.semi_major_axis / 1000:.1f} km")

# Load NASA Near Earth Network ground stations
print("\nLoading NASA Near Earth Network ground stations...")
start_time = time.time()
nen_stations = bh.datasets.groundstations.load("nasa nen")
elapsed = time.time() - start_time
print(f"Loaded {len(nen_stations)} NASA NEN ground stations in {elapsed:.2f} seconds.")

# Propagate NISAR for 3 orbits to create ground track visualization
print("\nPropagating NISAR for 3 orbits...")
start_time = time.time()
orbital_period = bh.orbital_period(nisar.semi_major_axis)
nisar.propagate_to(nisar.epoch + 3 * orbital_period)
elapsed = time.time() - start_time
print(f"Orbital period: {orbital_period / 60:.1f} minutes")
print(f"Propagated NISAR for 3 orbits in {elapsed:.2f} seconds.")

# Create ground track visualization with communication cones
print("\nCreating ground track visualization with communication cones...")
start_time = time.time()
fig_groundtrack = bh.plot_groundtrack(
    trajectories=[
        {
            "trajectory": nisar.trajectory,
            "color": "red",
            "line_width": 2.0,
            "track_length": 3,
            "track_units": "orbits",
        }
    ],
    ground_stations=[
        {
            "stations": nen_stations,
            "color": "blue",
            "alpha": 0.15,
            "point_size": 5.0,
        }
    ],
    gs_cone_altitude=nisar.semi_major_axis - bh.R_EARTH,
    gs_min_elevation=5.0,
    basemap="natural_earth",
    show_borders=True,
    show_coastlines=True,
    show_legend=False,
    backend="plotly",
)
elapsed = time.time() - start_time
print(f"Created ground track visualization in {elapsed:.2f} seconds.")

start_time = time.time()
# Reset propagator and compute 7-day access windows
print("\nComputing 7-day ground contacts...")
nisar.reset()

epoch_start = nisar.epoch
epoch_end = epoch_start + 7 * 86400.0  # 7 days in seconds

# Propagate for full 7-day period
nisar.propagate_to(epoch_end)

# Compute access windows with 5 degree minimum elevation
constraint = bh.ElevationConstraint(min_elevation_deg=5.0)
windows = bh.location_accesses(
    nen_stations, [nisar], epoch_start, epoch_end, constraint
)
elapsed = time.time() - start_time
print(f"Computed {len(windows)} contact windows in {elapsed:.2f} seconds.")

# Print sample of contact windows
print("\n" + "=" * 80)
print("Sample Contact Windows (first 10)")
print("=" * 80)
print(
    f"{'Station':<20} {'Start Time':<25} {'End Time':<25} {'Duration':>10} {'Max Elev':>10}"
)
print("-" * 80)
for i, window in enumerate(windows[:10]):
    duration_min = window.duration / 60.0
    max_elev = window.properties.elevation_max
    start_str = str(window.start).split(".")[0]  # Remove fractional seconds
    end_str = str(window.end).split(".")[0]
    print(
        f"{window.location_name:<20} {start_str:<25} {end_str:<25} {duration_min:>8.1f} m {max_elev:>8.1f}°"
    )
print("=" * 80)

# Export first 20 contact windows to CSV for documentation
csv_path = OUTDIR / f"{SCRIPT_NAME}_windows.csv"
with open(csv_path, "w", newline="") as csvfile:
    writer = csv.writer(csvfile)
    writer.writerow(
        [
            "Station",
            "Start Time (UTC)",
            "End Time (UTC)",
            "Duration (min)",
            "Max Elevation (deg)",
        ]
    )
    for window in windows[:20]:  # Only export first 20 for documentation
        duration_min = window.duration / 60.0
        max_elev = window.properties.elevation_max
        start_str = str(window.start).split(".")[0]  # Remove fractional seconds
        end_str = str(window.end).split(".")[0]
        writer.writerow(
            [
                window.location_name,
                start_str,
                end_str,
                f"{duration_min:.1f}",
                f"{max_elev:.1f}",
            ]
        )
print(f"✓ Exported first 20 contact windows to {csv_path}")

# Analyze contact statistics
print("\nAnalyzing contact statistics...")

# Group contacts by station
station_contacts = {}
for window in windows:
    station_name = window.location_name
    if station_name not in station_contacts:
        station_contacts[station_name] = []
    station_contacts[station_name].append(window)

# Calculate average daily contacts per station
days = 7.0
station_daily_avg = {}
for station, contacts in station_contacts.items():
    avg_per_day = len(contacts) / days
    station_daily_avg[station] = avg_per_day

# Sort by average daily contacts
sorted_stations = sorted(station_daily_avg.items(), key=lambda x: x[1], reverse=True)

print("\nAverage Daily Contacts per Station:")
print("-" * 40)
for station, avg in sorted_stations:
    total = len(station_contacts[station])
    print(f"{station:<20}: {avg:>5.1f} contacts/day ({total} total)")


# Figure 1: Daily contacts per station (bar chart)
stations_list = [s[0] for s in sorted_stations]
daily_avgs = [s[1] for s in sorted_stations]

fig_daily = go.Figure(
    data=[
        go.Bar(
            x=stations_list,
            y=daily_avgs,
            marker_color="steelblue",
            hovertemplate="<b>%{x}</b><br>%{y:.1f} contacts/day<extra></extra>",
        )
    ]
)
fig_daily.update_layout(
    title="NISAR Average Daily Contacts by Station (7-day period)",
    xaxis_title="Ground Station",
    yaxis_title="Average Daily Contacts",
    xaxis_tickangle=-45,
    height=700,
    margin=dict(l=60, r=40, t=80, b=120),
)

# Figure 2: Contact duration distribution (histogram)
durations = [w.duration / 60.0 for w in windows]  # Convert to minutes

mean_duration = np.mean(durations)
median_duration = np.median(durations)
max_duration = np.max(durations)

fig_duration = go.Figure(
    data=[
        go.Histogram(
            x=durations,
            nbinsx=30,
            marker_color="coral",
            marker_line_color="black",
            marker_line_width=1,
            hovertemplate="Duration: %{x:.1f} min<br>Count: %{y}<extra></extra>",
        )
    ]
)
fig_duration.update_layout(
    title="NISAR Contact Duration Distribution",
    xaxis_title="Contact Duration (minutes)",
    yaxis_title="Frequency",
    height=700,
    margin=dict(l=60, r=40, t=80, b=60),
    annotations=[
        dict(
            text=f"Mean: {mean_duration:.1f} min<br>Median: {median_duration:.1f} min<br>Max: {max_duration:.1f} min",
            xref="paper",
            yref="paper",
            x=0.05,
            y=0.97,
            xanchor="left",
            yanchor="top",
            showarrow=False,
            bordercolor="grey",
            borderwidth=1,
            borderpad=8,
        )
    ],
)

print("\nContact Duration Statistics:")
print(f"  Mean: {mean_duration:.1f} minutes")
print(f"  Median: {median_duration:.1f} minutes")
print(f"  Min: {np.min(durations):.1f} minutes")
print(f"  Max: {max_duration:.1f} minutes")

See Also