Seminar 6¶
In this seminar, we will take a look on how to compute the path of active satellites and space stations.
In [1]:
import wget
import numpy as np
import ephem
import datetime
import progressbar
from matplotlib import pylab as plt
ISS Two Line Elements (TLE)¶
We will download the two-line elements for all space stations in the current directory, using pywget as follows.
In [2]:
wget.download("https://celestrak.com/NORAD/elements/stations.txt")
Out[2]:
Fortunately, the ISS information (NAME: Zarya -> First module of the ISS and two-line elements) are contained in the first three lines. We can extract the first three lines as follows:
In [3]:
f = open("stations.txt")
L1= f.readline()
L2= f.readline()
L3= f.readline()
f.close()
print L1
print L2
print L3
We can parse the TLE into an ephem object by using the readtle() function, to quickly compute the horizontal coordinates for a given time and observer. Besides the ISS we will also initialize Jupyter as an ephem object.
In [4]:
ISS = ephem.readtle(L1,L2,L3)
JUPITER = ephem.Jupiter()
Define Observer:
In [5]:
nhao = ephem.Observer()
nhao.lat = "35.02528"
nhao.lon = "134.33556"
nhao.elevation = 449
nhao.date = datetime.datetime.now()-datetime.timedelta(hours=9)
In [6]:
NOW = datetime.datetime.now()-datetime.timedelta(hours=9)
str(NOW)
Out[6]:
Define a time range $\pm 1 $ Day from now on.
In [7]:
dt = []
for i in range(-60*24,60*24):
dt.append(NOW+datetime.timedelta(minutes=1*i))
The following Function will obtain the horizontal coordinates (ALT-AZ) for a given OBJECT and a given time interval (DELTATIME).
In [8]:
def get_ALT_AZ(OBJ, DELTATIME):
ALT, AZ = [], []
for i in DELTATIME:
nhao.date = i
OBJ.compute(nhao)
ALT.append(np.rad2deg(OBJ.alt))
AZ.append(np.rad2deg(OBJ.az))
return ALT, AZ
JUP_alt, JUP_az = get_ALT_AZ(JUPITER, dt)
ISS_alt, ISS_az = get_ALT_AZ(ISS , dt)
Altitude-Time Diagramm¶
In [9]:
NOW = datetime.datetime.now()-datetime.timedelta(hours=9)
plt.figure(figsize = (12,6))
plt.plot(dt,ISS_alt, ".", label = "ISS")
plt.plot(dt,JUP_alt, ".", label = "Jupiter")
plt.axvline(NOW, color = "Red")
plt.text(NOW, 80, "NOW", color = "Red", rotation = 90, fontsize =18)
plt.axhline( 0, color = "k")
plt.axhline(30)
plt.axhspan(-90,0, color = "cornflowerblue", alpha = 0.7)
plt.axhspan(0,30, color = "cornflowerblue", alpha = 0.2)
plt.xlim([dt[0],dt[-1]])
plt.ylim([-90,90])
plt.grid(True)
plt.xlabel("Time [UTC]")
plt.ylabel("Elevation [deg]")
plt.legend()
plt.show()
Altitude-Azimuth Diagramm¶
In [10]:
plt.figure(figsize = (12,6))
plt.plot(ISS_az,ISS_alt, ".", label = "ISS")
plt.plot(JUP_az,JUP_alt, ".", label = "Jupiter")
plt.axhline( 0, color = "k")
plt.axhline(30)
plt.axhspan(-90,0, color = "cornflowerblue", alpha = 0.7)
plt.axhspan(0,30, color = "cornflowerblue", alpha = 0.2)
plt.xlim([0,360])
plt.ylim([-90,90])
plt.grid(True)
plt.xlabel("Azimuth [deg]")
plt.ylabel("Elevation [deg]")
plt.legend()
plt.show()
More Satellites¶
Download TLE list:
In [11]:
wget.download("https://celestrak.com/NORAD/elements/active.txt")
Out[11]:
In [12]:
satellite_path = "active.txt"
In [13]:
def get_all_TLE(PATH):
""" Loads a TLE file and creates a list of satellites.
Input: file path
Output: List of ephem objects"""
f = open(PATH)
satlist = []
l1 = f.readline()
while l1:
l2 = f.readline()
l3 = f.readline()
sat = ephem.readtle(l1,l2,l3)
satlist.append(sat)
l1 = f.readline()
f.close()
return satlist
SAT_LIST = get_all_TLE(satellite_path)
In [14]:
"Number of Satellites:", len(SAT_LIST)
Out[14]:
Since the number of satellites is quite high, let`s limit the time interval to $\pm$ 6 hours.
In [15]:
dt2 = []
for i in range(-60*6,60*6):
dt2.append(NOW+datetime.timedelta(minutes=1*i))
In [16]:
def get_all_ALT_AZ(OBJECTS, TIMEDELTA):
"""Computes ALT AZ for a List of pyephem objects
Input: List of pyephem objects
Output: two 2D numpy arrays (ALT[alt,time], AZ[az,time]) """
ALL_ALT, ALL_AZ = [],[]
bar = progressbar.ProgressBar(max_value=len(OBJECTS))
for i in range(len(OBJECTS)):
ALT, AZ = get_ALT_AZ(OBJECTS[i], TIMEDELTA)
ALL_ALT.append(ALT)
ALL_AZ.append(AZ)
bar.update(i+1)
return np.asarray(ALL_ALT),np.asarray(ALL_AZ)
ALL_ALT, ALL_AZ = get_all_ALT_AZ(SAT_LIST, dt2)
In [17]:
#### Get AZ and ALT for all stations and Jupiter
STA_LIST = get_all_TLE("stations.txt")
STA_ALT, STA_AZ = get_all_ALT_AZ(STA_LIST, dt2)
JUP_ALT, JUP_AZ = get_ALT_AZ(JUPITER, dt2)
Altitude-Azimuth plot of all active satellites¶
In [18]:
plt.figure(figsize = (12,8))
plt.axhline(30,linewidth = 1, linestyle = "--")
plt.axhline(0)
plt.axhspan(-90,0, color = "cornflowerblue", alpha = 0.3)
plt.axhspan(0,30 , color = "cornflowerblue", alpha = 0.1)
plt.scatter(ALL_AZ[:,0], ALL_ALT[:,0], 30 , color = "Darkred" , alpha = 0.2, label = "Satellites")
plt.scatter(STA_AZ[:,0], STA_ALT[:,0], 40 , color = "Darkgreen" , alpha = 1 , label = "Stations")
plt.scatter(STA_AZ[0,0], STA_ALT[0,0], 100, color = "cornflowerblue" , label = "ISS")
plt.scatter(JUP_AZ[0] , JUP_ALT[0] , 200, color = "Darkorange" , label = "Jupiter")
plt.title("Satellites at "+str(dt2[0])+"[UTC] at NHAO")
plt.xlim([0,360])
plt.ylim([-90,90])
plt.grid(True)
plt.xlabel("Azimuth [deg]")
plt.ylabel("Elevation [deg]")
plt.legend(loc="upper right")
plt.show()
Time propagation Animation¶
In [19]:
from matplotlib.animation import ArtistAnimation
from IPython.display import HTML
We can use the ArtistAnimation function included in the matplotlib function to produce an animation of sattelite movement from the time dependent coordinate layers.
In [20]:
fig , ax = plt.subplots(figsize = (12,8))
plt.axhline(30,linewidth = 1, linestyle = "--")
plt.axhline(0)
plt.axhspan(-90,0, color = "cornflowerblue", alpha = 0.2)
plt.axhspan(0,30 , color = "cornflowerblue", alpha = 0.1)
plt.xlim([0,360])
plt.ylim([-90,90])
plt.grid(True)
plt.xlabel("Azimuth [deg]")
plt.ylabel("Elevation [deg]")
FRAMES =[]
for i in range(len(dt2)):
SATS = ax.plot(ALL_AZ[:,i], ALL_ALT[:,i], ".", color = "Darkred", alpha = 0.2)
STAT = ax.plot(STA_AZ[:,i], STA_ALT[:,i], ".", color = "Darkgreen" )
SISS = ax.plot(STA_AZ[0,i], STA_ALT[0,i], "o", color = "Darkgreen")
SJUP = ax.plot(JUP_AZ[i], JUP_ALT[i], "o", color = "Darkorange")
JUPT = ax.text(JUP_AZ[i], JUP_ALT[i],"Jupiter", fontsize = 14)
ISST = ax.text(STA_AZ[0,i], STA_ALT[0,i],"ISS", fontsize = 14)
TEXT = ax.text(2,83,"Satellites on "+str(dt2[i])+"[UTC] at NHAO", fontsize = 14)
FRAMES.append(SATS+STAT+SISS+SJUP+[JUPT]+[ISST]+[TEXT])
ani = ArtistAnimation(fig, FRAMES, interval=50)
print "Done"
#ani.save('Satellites.gif', writer="imagemagick") ### saves the animation to a .gif file
plt.close()
HTML(ani.to_html5_video())
Out[20]: