Notebook

Working with satellites¶

License¶

Working with the satellites module of pycraf.
Copyright (C) 2015+  Benjamin Winkel (bwinkel@mpifr.de)

This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program; if not, write to the Free Software
Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

Note: GitHub's integrated notebook viewer won't show the embedded videos. You can use the nbviewer service to open the notebook. Or follow the direct link.

In [1]:

%matplotlib inline

In [2]:

import requests
import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import animation, rc

from pycraf import satellite
from astropy import time
from astropy import units as u
from astropy.coordinates import EarthLocation

rc('animation', html='html5')

Two-line elements sets (TLE)¶

To compute positions of a satellite, two-line element sets (TLE) have to be used. These look like the following

ISS (ZARYA)
1 25544U 98067A   13165.59097222  .00004759  00000-0  88814-4 0    47
2 25544  51.6478 121.2152 0011003  68.5125 263.9959 15.50783143834295

which is the latest TLE for the international space station. The first line has only the name (or identifier) of a satellite. The following two lines contain all the necessary information about the satellite orbit to calculate its position for a certain time. Note that the TLEs are usually published once a day, because the contained parameters quickly change; drag forces cause rapid changes in the orbits of almost all satellites.

Of course, the position of a satellite depends on the choosen coordinate system. With the help of the sgp4 package pycraf calculates the ECI position of the satellite (in Cartesian coordinates). ECI is the Earth-centered inertial frame. But often we are interested in the apparent sky position (horizontal frame) of a satellite with respect to an observer. The latter is not provided by the sgp4 library, but pycraf adds this functionality.

Fetching latest satellites from Celestrak¶

As an example how to use the pycraf.satellites module, we will download the latest TLE data of all satellites launched in the last 30 days.

In [3]:

celestrak_latest = requests.get('http://celestrak.com/NORAD/elements/tle-new.txt')

In [4]:

celestrak_latest.text[:1000]

Out[4]:

'FENGYUN 4B              \r\n1 48808U 21047A   21179.92803759 -.00000367  00000-0  00000-0 0  9992\r\n2 48808   0.0541 226.6478 0000694 252.0694 256.3446  1.00272877   387\r\nCZ-3B R/B               \r\n1 48809U 21047B   21179.71501365  .00048070  97412-6  15215-2 0  9995\r\n2 48809  28.6217 122.6465 7215152 195.7071 116.8720  2.40687722   638\r\nDRAGON CRS-22           \r\n1 48831U 21048A   21179.59233343  .00001944  00000-0  43787-4 0  9992\r\n2 48831  51.6454 286.3205 0002279 118.9439   5.7591 15.48743536290365\r\nSXM-8                   \r\n1 48838U 21049A   21180.57409069  .00000001  00000-0  00000+0 0  9999\r\n2 48838   0.0430 263.7940 0002147 166.6985 293.5692  1.00273981   441\r\nFALCON 9 R/B            \r\n1 48839U 21049B   21180.25559672  .00006015  00000-0  42660-3 0  9994\r\n2 48839  27.0231 135.4209 5912344 205.6044 106.4870  4.22186156   979\r\n2021-050A               \r\n1 48840U 21050A   21180.43817509 -.00000073  00000-0  00000+0 0  9991\r\n2 48840  97.5118 250.1699 0011777 203.6227 211.1975 15.23623308'

The 'tle-new.txt' file apparently simply contains a list of TLEs. Let's fix the line endings (\r\n to \n) and split into a list of TLEs:

In [5]:

all_lines = celestrak_latest.text.split('\r\n')
tle_list = []
for idx in range(len(all_lines) // 3):
    tle_list.append('\n'.join(all_lines[3 * idx: 3 * (idx + 1)]))

In [6]:

print(tle_list[0])
print(tle_list[1])

FENGYUN 4B              
1 48808U 21047A   21179.92803759 -.00000367  00000-0  00000-0 0  9992
2 48808   0.0541 226.6478 0000694 252.0694 256.3446  1.00272877   387
CZ-3B R/B               
1 48809U 21047B   21179.71501365  .00048070  97412-6  15215-2 0  9995
2 48809  28.6217 122.6465 7215152 195.7071 116.8720  2.40687722   638

Plotting satellite position¶

Before we start plotting, we'll use have sgp4 parse all the satellite strings and construct sgp4.io.Satellite instances. These merely contain the orbital parameters, but it will be faster to work the the already parsed data instead of letting the TLEs being parsed in each time step (yes, we are going to do a nice animation).

In [7]:

sat_list = [satellite.get_sat(tle_string) for tle_string in tle_list]
# contains a list of tuples: (sat_name, sat_instance)
sat_list[0]

Out[7]:

('FENGYUN 4B              ', <sgp4.model.Satrec at 0x55a5f3a1a640>)

The second ingredient, we need is an array of "observing" times. pycraf wants an astropy.time.Time object for this, because it has (accurate!) built-in conversion between UTC and sidereal time.

In [8]:

start_time = time.Time.now()
print('Start time:', start_time)
td = time.TimeDelta(np.arange(0, 3600 * 24, 60 * 1), format='sec')
times = start_time + td  # one day in steps of 1 min

Start time: 2021-06-29 18:32:03.306425

ECI position¶

For this we make use of sgp4 built-in propagate method. This is easy to use, but it doesn't work with numpy arrays. Let's vectorize it:

In [9]:

def _propagate(sat, dt):
    '''
    True equator mean equinox (TEME) position from `sgp4` at given time.

    Parameters
    ----------
    sat : `sgp4.io.Satellite` instance
        Satellite object filled from TLE
    dt : `~datetime.datetime`
        Time

    Returns
    -------
    xs, ys, zs : float
        TEME (=True equator mean equinox) position of satellite [km]
    '''

    # pos [km], vel [km/s]

    try:
        from sgp4.api import jday, SGP4_ERRORS
    except ImportError:
        raise ImportError(
            'The "sgp4" package (version 2+) is necessary to use this '
            'function.'
            )

    jd, fr = jday(
        dt.year, dt.month, dt.day,
        dt.hour, dt.minute, dt.second + dt.microsecond / 1e6
        )
    err_code, position, velocity = sat.sgp4(jd, fr)

    if err_code:
        raise ValueError(
            'Satellite propagation error', err_code,
            '({})'.format(SGP4_ERRORS[err_code])
            )

    return position


vec_propagate = np.vectorize(
    _propagate, excluded=['sat'], otypes=[np.float64] * 3
    )

In [10]:

pos_eci = np.array([vec_propagate(sat[1], times.datetime) for sat in sat_list])

In [11]:

plim = 6000

fig = plt.figure(figsize=(12, 9))
ax = Axes3D(fig)
ax.view_init(azim=60, elev=30)
# Aspect ratio is  not implemented;
# see https://github.com/matplotlib/matplotlib/issues/1077/
# ax.set_aspect('equal')
# need to manually stretch to make it approx. right
ax.set_xlim((-2 * plim, 2 * plim))
ax.set_ylim((-2 * plim, 2 * plim))
ax.set_zlim((-plim, plim))
ax.set_xlabel('x [km]')
ax.set_ylabel('y [km]')
ax.set_zlabel('z [km]')

points = ax.scatter(pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0], marker='o')
title = ax.set_title('{:%y/%m/%d %H:%M}'.format(times[0].datetime), loc='center', fontsize=20)

def init():
    points._offsets3d = (pos_eci[:, 0, 0], pos_eci[:, 1, 0], pos_eci[:, 2, 0])
    title.set_text('{:%y/%m/%d %H:%M}'.format(times[0].datetime))
    return points, title

def animate(i):
    points._offsets3d = (pos_eci[:, 0, i], pos_eci[:, 1, i], pos_eci[:, 2, i])
    title.set_text('{:%y/%m/%d %H:%M}'.format(times[i].datetime))
    return points, title

# call the animator. blit=True means only re-draw the parts that have changed.
anim = animation.FuncAnimation(
    fig, animate, init_func=init, frames=pos_eci.shape[2], interval=20, blit=True
    )

# this takes a while!
plt.close(anim._fig)
anim

Out[11]:

Horizontal position¶

For the next plot, we need to define an observer position, say the famous Parkes 64-m dish and the Effelsberg 100-m radio telescope.

In [12]:

obs_loc_parkes = EarthLocation(148.25738, -32.9933, 414.8)
obs_loc_effbg = EarthLocation(6.88375, 50.525, 366.)

# create a SatelliteObserver instance
sat_obs_parkes = satellite.SatelliteObserver(obs_loc_parkes)
sat_obs_effbg = satellite.SatelliteObserver(obs_loc_effbg)
pos_horiz_parkes = np.array([
    sat_obs_parkes.azel_from_sat(sat[1], times)
    for sat in sat_list
    ])
pos_horiz_effbg = np.array([
    sat_obs_effbg.azel_from_sat(sat[1], times)
    for sat in sat_list
    ])

In [13]:

vmin, vmax = np.log10(100), np.log10(100000)

fig = plt.figure(figsize=(12, 7))
ax1 = fig.add_axes((0.1, 0.5, 0.8, 0.35))
ax2 = fig.add_axes((0.1, 0.1, 0.8, 0.35))
cax = fig.add_axes((0.91, 0.2, 0.02, 0.5))
ax2.set_xlabel('Azimuth [deg]')
ax1.set_ylabel('Elevation [deg]')
for ax in [ax1, ax2]:
    ax.set_xlim((-180, 180))
    ax.set_ylim((0, 90))
    ax.set_xticks(range(-150, 180, 30))
    ax.set_yticks(range(0, 91, 30))
    ax.set_aspect('equal')
    ax.grid()

points1 = ax1.scatter(
    pos_horiz_parkes[:, 0, 0], pos_horiz_parkes[:, 1, 0],
    c=np.log10(pos_horiz_parkes[:, 2, 0]),
    cmap='viridis', vmin=vmin, vmax=vmax,
    )
points2 = ax2.scatter(
    pos_horiz_effbg[:, 0, 0], pos_horiz_effbg[:, 1, 0],
    c=np.log10(pos_horiz_effbg[:, 2, 0]),
    cmap='viridis', vmin=vmin, vmax=vmax,
    )
cbar = fig.colorbar(points1, cax=cax)
cbar.set_label('Distance (km)')
cbar.set_ticks([2, 3, 4, 5])
cbar.set_ticklabels([100, 1000, 10000, 100000])

ax1.text(-170, 75, 'Parkes 64-m', fontsize=16)
ax2.text(-170, 75, 'Effelsberg 100-m', fontsize=16)
title = ax1.text(
    174, 75, '{:%y/%m/%d %H:%M}'.format(times[0].datetime),
    fontsize=15, ha='right'
    )

def init():
    points1.set_offsets(pos_horiz_parkes[:, 0:2, 0])
    points1.set_array(np.log10(pos_horiz_parkes[:, 2, 0]))
    points2.set_offsets(pos_horiz_effbg[:, 0:2, 0])
    points2.set_array(np.log10(pos_horiz_effbg[:, 2, 0]))
    title.set_text('{:%y/%m/%d %H:%M}'.format(times[0].datetime))
    return points, title

def animate(i):
    points1.set_offsets(pos_horiz_parkes[:, 0:2, i])
    points1.set_array(np.log10(pos_horiz_parkes[:, 2, i]))
    points2.set_offsets(pos_horiz_effbg[:, 0:2, i])
    points2.set_array(np.log10(pos_horiz_effbg[:, 2, i]))
    title.set_text('{:%y/%m/%d %H:%M}'.format(times[i].datetime))
    return points, title

anim = animation.FuncAnimation(
    fig, animate, init_func=init, frames=pos_horiz_parkes.shape[2], interval=20, blit=True
    )

# this takes a while!
plt.close(anim._fig)
anim

Out[13]:

In [ ]: