Satellite Ground Contacts

Compute communications contact intervals of “LandSat 7” saetllite over a single day

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 mininum 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},
]

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 statoin
    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
    # coordinated, 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 arcsign of the "up" portion of the NED-hat vetor
    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]

# 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  })

	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

# 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
                )