Satellite Ground Contacts¶

This tutorial computes communication contact intervals between a satellite and a set of ground stations over a single day.

A satellite is "in view" of a ground station when it rises above a minimum elevation angle (here, 5 degrees). For each second of the day, we compute the satellite's position relative to each ground station in the North-East-Down (NED) frame, determine the elevation angle, and identify contiguous time intervals where the satellite is visible.

The example uses the Landsat-7 satellite and three ground stations: Svalbard, Alice Springs, and Sioux Falls.

In [1]:

import satkit as sk
import numpy as np
import plotly.express as px
from datetime import datetime, timedelta

# The TLE for landsat-7
tle_lines = [
    "0 LANDSAT-7",
    "1 25682U 99020A   24099.90566066  .00000551  00000-0  12253-3 0  9992",
    "2 25682  97.8952 129.9471 0001421 108.5441  14.5268 14.60548156329087"
]
landsat7 = sk.TLE.from_lines(tle_lines)

# The minimum elevation for a contact
min_elevation_deg = 5


# The date to compute ground contacts: April 9, 2024
date = datetime(2024, 4, 9)
# Any array of times representing every second of the day
time_array = np.array([date + timedelta(seconds=x) for x in range(86400)])

# Get satellite positions in TEME frame (pseudo-inertial) via SGP4
pTEME, _vTEME = sk.sgp4(landsat7, time_array)

# Get ITRF coordinates (Earth-Fixed) by rotating the position in the TEME frame
# to ITRF frame using the frametransform module
pITRF = np.array([q*x for q,x in zip(sk.frametransform.qteme2itrf(time_array), pTEME)])

# Setup some ground stations
ground_stations = [
    {'name': 'Svalbard', 'lat': 78.2232, 'lon': 15.6267, 'alt': 0},
    {'name': 'Alice Springs', 'lat': -23.6980, 'lon': 133.8807, 'alt': 0},
    {'name': 'Sioux Falls', 'lat': 43.5446, 'lon': -96.7311, 'alt': 0},
]

import satkit as sk import numpy as np import plotly.express as px from datetime import datetime, timedelta # The TLE for landsat-7 tle_lines = [ "0 LANDSAT-7", "1 25682U 99020A 24099.90566066 .00000551 00000-0 12253-3 0 9992", "2 25682 97.8952 129.9471 0001421 108.5441 14.5268 14.60548156329087" ] landsat7 = sk.TLE.from_lines(tle_lines) # The minimum elevation for a contact min_elevation_deg = 5 # The date to compute ground contacts: April 9, 2024 date = datetime(2024, 4, 9) # Any array of times representing every second of the day time_array = np.array([date + timedelta(seconds=x) for x in range(86400)]) # Get satellite positions in TEME frame (pseudo-inertial) via SGP4 pTEME, _vTEME = sk.sgp4(landsat7, time_array) # Get ITRF coordinates (Earth-Fixed) by rotating the position in the TEME frame # to ITRF frame using the frametransform module pITRF = np.array([q*x for q,x in zip(sk.frametransform.qteme2itrf(time_array), pTEME)]) # Setup some ground stations ground_stations = [ {'name': 'Svalbard', 'lat': 78.2232, 'lon': 15.6267, 'alt': 0}, {'name': 'Alice Springs', 'lat': -23.6980, 'lon': 133.8807, 'alt': 0}, {'name': 'Sioux Falls', 'lat': 43.5446, 'lon': -96.7311, 'alt': 0}, ]

Compute Satellite Positions¶

Propagate the TLE with SGP4 to get TEME positions at every second of the day, then rotate to the ITRF (Earth-fixed) frame. Define the ground stations by their geodetic coordinates.

Contact Detection Algorithm¶

For each ground station, the algorithm:

1. Computes the satellite position in the station's North-East-Down (NED) frame by subtracting the station position and rotating with qned2itrf
2. Computes elevation as the arcsine of the upward (-Down) component of the normalized NED vector
3. Finds all times where elevation exceeds the minimum threshold
4. Groups consecutive qualifying indices into individual contact intervals
5. For each contact, records the start/end times, duration, range, elevation, and heading

Contact Summary¶

Display all contacts for the day sorted by start time. Higher maximum elevation generally corresponds to better signal quality.

Contact Detail¶

Plot the range, elevation, and heading over time for an individual contact pass.

In [2]:

def calc_contacts(ground_station, pITRF, time_array):
    """
    Compute contact times for a single ground station given satellite position in Earth-fixed frame
    """

    # Create an "itrfcoord" object for the ground station
    coord = sk.itrfcoord(latitude_deg=ground_station['lat'], longitude_deg=ground_station['lon'], altitude_m=ground_station['alt'])

    # Get the North-East-Down coordinates of the satellite relative to the ground station
    # at all times by taking the difference between the satellite position and the ground
    # coordinates, then rotating to the "North-East-Down" frame relative to the ground station
    pNED = np.array([coord.qned2itrf.conj * (x - coord.vector) for x in pITRF])

    # Normalize the NED coordinates
    pNED_hat = pNED / np.linalg.norm(pNED, axis=1)[:, None]

    # Find the elevation from the ground station at all times
    # This is the arcsine of the "up" portion of the NED-hat vector
    elevation_deg = np.degrees(np.arcsin(-pNED_hat[:,2]))

    # We can see ground station when elevation is greater than min_elevation_deg
    inview_idx = np.argwhere(elevation_deg > min_elevation_deg).flatten().astype(int)

    # split indices into groups of consecutive indices
    # This indicates contiguous contacts
    inview_idx = np.split(inview_idx, np.where(np.diff(inview_idx) != 1)[0]+1)

    def get_single_contacts(inview_idx):
        for cidx in inview_idx:
            # cidx are indices to the time array for this contact

            # the North-East-Down position of the satellite relative to
            # ground station over the single contact
            cpNED = pNED[cidx,:]

            # Compute the range in meters
            range = np.linalg.norm(cpNED, axis=1)

            # elevation in degrees over the contact
            contact_elevation_deg = elevation_deg[cidx]

            # Heading clockwise from North is arctangent of east/north'
            heading_deg = np.degrees(np.arctan2(cpNED[:,1], cpNED[:,0]))

            # Yield a dictionary describing the results
            yield {
                'groundstation': ground_station['name'],
                'timearr': time_array[cidx],
                'range_km': range*1.0e-3,
                'elevation_deg': contact_elevation_deg,
                'heading_deg': heading_deg,
                'start': time_array[cidx[0]],
                'end': time_array[cidx[-1]],
                'max_elevation_deg': np.max(contact_elevation_deg),
                'duration': time_array[cidx[-1]] - time_array[cidx[0]]
            }
    return list(get_single_contacts(inview_idx))

# Calculate all the contacts
contacts = [calc_contacts(g, pITRF, time_array) for g in ground_stations]

# Flatten contacts into 1D list
contacts = [item for sublist in contacts for item in sublist]

def calc_contacts(ground_station, pITRF, time_array): """ Compute contact times for a single ground station given satellite position in Earth-fixed frame """ # Create an "itrfcoord" object for the ground station coord = sk.itrfcoord(latitude_deg=ground_station['lat'], longitude_deg=ground_station['lon'], altitude_m=ground_station['alt']) # Get the North-East-Down coordinates of the satellite relative to the ground station # at all times by taking the difference between the satellite position and the ground # coordinates, then rotating to the "North-East-Down" frame relative to the ground station pNED = np.array([coord.qned2itrf.conj * (x - coord.vector) for x in pITRF]) # Normalize the NED coordinates pNED_hat = pNED / np.linalg.norm(pNED, axis=1)[:, None] # Find the elevation from the ground station at all times # This is the arcsine of the "up" portion of the NED-hat vector elevation_deg = np.degrees(np.arcsin(-pNED_hat[:,2])) # We can see ground station when elevation is greater than min_elevation_deg inview_idx = np.argwhere(elevation_deg > min_elevation_deg).flatten().astype(int) # split indices into groups of consecutive indices # This indicates contiguous contacts inview_idx = np.split(inview_idx, np.where(np.diff(inview_idx) != 1)[0]+1) def get_single_contacts(inview_idx): for cidx in inview_idx: # cidx are indices to the time array for this contact # the North-East-Down position of the satellite relative to # ground station over the single contact cpNED = pNED[cidx,:] # Compute the range in meters range = np.linalg.norm(cpNED, axis=1) # elevation in degrees over the contact contact_elevation_deg = elevation_deg[cidx] # Heading clockwise from North is arctangent of east/north' heading_deg = np.degrees(np.arctan2(cpNED[:,1], cpNED[:,0])) # Yield a dictionary describing the results yield { 'groundstation': ground_station['name'], 'timearr': time_array[cidx], 'range_km': range*1.0e-3, 'elevation_deg': contact_elevation_deg, 'heading_deg': heading_deg, 'start': time_array[cidx[0]], 'end': time_array[cidx[-1]], 'max_elevation_deg': np.max(contact_elevation_deg), 'duration': time_array[cidx[-1]] - time_array[cidx[0]] } return list(get_single_contacts(inview_idx)) # Calculate all the contacts contacts = [calc_contacts(g, pITRF, time_array) for g in ground_stations] # Flatten contacts into 1D list contacts = [item for sublist in contacts for item in sublist]

In [3]:

# Convert to pandas dataframe for nice table display

import pandas as pd

data = pd.DataFrame(contacts)
data.sort_values(by='start', inplace=True)
data.reset_index(drop=True, inplace=True)
# Get nicer column names for display
data.rename(columns={"max_elevation_deg": "Max Elevation (deg)",
                     "duration": "Duration (s)",
                     "start": "Start (UTC)",
                     "end": "End (UTC)",
                     "groundstation": "Ground Station"}, inplace=True)
data.style \
    .hide(subset=["timearr", "range_km", "elevation_deg", "heading_deg"], axis=1) \
    .format({"Max Elevation (deg)": "{:.1f}",
             "Start (UTC)": lambda x: x.strftime('%H:%M:%S'),
             "End (UTC)": lambda x: x.strftime('%H:%M:%S'),
               "Duration (s)": lambda x: x.seconds  })

# Convert to pandas dataframe for nice table display import pandas as pd data = pd.DataFrame(contacts) data.sort_values(by='start', inplace=True) data.reset_index(drop=True, inplace=True) # Get nicer column names for display data.rename(columns={"max_elevation_deg": "Max Elevation (deg)", "duration": "Duration (s)", "start": "Start (UTC)", "end": "End (UTC)", "groundstation": "Ground Station"}, inplace=True) data.style \ .hide(subset=["timearr", "range_km", "elevation_deg", "heading_deg"], axis=1) \ .format({"Max Elevation (deg)": "{:.1f}", "Start (UTC)": lambda x: x.strftime('%H:%M:%S'), "End (UTC)": lambda x: x.strftime('%H:%M:%S'), "Duration (s)": lambda x: x.seconds })

Out[3]:

	Ground Station	Start (UTC)	End (UTC)	Max Elevation (deg)	Duration (s)
0	Sioux Falls	00:35:13	00:45:45	32.6	632
1	Svalbard	00:50:29	00:55:06	7.1	277
2	Sioux Falls	02:12:52	02:23:13	27.6	621
3	Svalbard	02:29:38	02:36:04	9.5	386
4	Svalbard	04:08:08	04:16:54	15.6	526
5	Svalbard	05:46:24	05:56:51	27.4	627
6	Svalbard	07:24:29	07:35:47	50.7	678
7	Svalbard	09:02:21	09:13:52	88.4	691
8	Alice Springs	10:08:08	10:18:52	33.7	644
9	Svalbard	10:39:59	10:51:21	63.4	682
10	Alice Springs	11:46:26	11:54:59	15.7	513
11	Svalbard	12:17:23	12:28:38	54.7	675
12	Sioux Falls	12:34:01	12:37:54	6.5	233
13	Svalbard	13:54:40	14:06:01	62.1	681
14	Sioux Falls	14:08:38	14:20:01	69.3	683
15	Svalbard	15:32:08	15:43:38	88.8	690
16	Sioux Falls	15:47:13	15:55:26	14.7	493
17	Svalbard	17:10:09	17:21:29	53.1	680
18	Svalbard	18:49:01	18:59:33	28.7	632
19	Svalbard	20:28:53	20:37:48	16.2	535
20	Alice Springs	21:04:50	21:12:26	12.5	456
21	Svalbard	22:09:40	22:16:17	9.8	397
22	Alice Springs	22:40:14	22:51:17	42.2	663
23	Sioux Falls	23:39:43	23:46:08	9.9	385
24	Svalbard	23:50:41	23:55:21	7.1	280

In [4]:

# Plot one of the contacts
contact = data.iloc[5]
import plotly.graph_objects as go
from plotly.subplots import make_subplots

fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.1,
                        subplot_titles=('Range [km]', 'Elevation [deg]', 'Heading [deg]'))
fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['range_km'], name='Range [km]'), row=1, col=1)
fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['elevation_deg'], name='Elevation [deg]'), row=2, col=1)
fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['heading_deg'], name='Heading [deg]'), row=3, col=1)
fig.update_layout(yaxis = {'title': 'Range (km)'},
                yaxis2={'title': 'Elevation [deg]'},
                yaxis3={'title': 'Heading [deg]'},
                title=f'Landsat 7 to {contact["Ground Station"]} on {contact["Start (UTC)"]}',
                width=800,
                height=600
                )

# Plot one of the contacts contact = data.iloc[5] import plotly.graph_objects as go from plotly.subplots import make_subplots fig = make_subplots(rows=3, cols=1, shared_xaxes=True, vertical_spacing=0.1, subplot_titles=('Range [km]', 'Elevation [deg]', 'Heading [deg]')) fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['range_km'], name='Range [km]'), row=1, col=1) fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['elevation_deg'], name='Elevation [deg]'), row=2, col=1) fig.add_trace(go.Scatter(x=contact['timearr'], y=contact['heading_deg'], name='Heading [deg]'), row=3, col=1) fig.update_layout(yaxis = {'title': 'Range (km)'}, yaxis2={'title': 'Elevation [deg]'}, yaxis3={'title': 'Heading [deg]'}, title=f'Landsat 7 to {contact["Ground Station"]} on {contact["Start (UTC)"]}', width=800, height=600 )