-
Notifications
You must be signed in to change notification settings - Fork 0
/
genPlots.py
153 lines (121 loc) · 6.05 KB
/
genPlots.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
## Alex Philpott
## Generates fancy plots of debris distribution
import math
import argparse
import pickle
import os
import numpy as np
from bokeh.plotting import figure, output_file, show, ColumnDataSource
from bokeh.models import HoverTool, Legend, LogColorMapper, ColorBar, LogTicker
from bokeh.layouts import gridplot
from bokeh.palettes import Blues9
from astroUtils import parse_catalog, Re
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="Input arguments")
optional = parser.add_argument_group('optional arguments')
optional.add_argument("--parse_data", "-p", help="Re-parse sat cat file, rather than read stored data",
required=False, action="store_true", default=False)
#required = parser.add_argument_group('optional arguments')
args = vars(parser.parse_args())
parse = args['parse_data']
dataDir = 'Data'
savefile = 'satcat3.pickle'
if parse or savefile not in os.listdir(os.path.join(os.getcwd(),dataDir)):
print('Parsing tle file...')
objects = parse_catalog(3)
with open(os.path.join(os.getcwd(),dataDir,savefile), 'wb') as f:
pickle.dump(objects,f)
else:
print('Loading pickle file')
with open(os.path.join(os.getcwd(),dataDir,savefile),'rb') as f:
objects = pickle.load(f)
## Check that Plots folder is made
if not 'Plots' in os.listdir(os.getcwd()):
os.makedirs(os.getcwd()+'/Plots')
if not 'Data' in os.listdir(os.getcwd()):
os.makedirs(os.getcwd()+'/Data')
lower_alt = 200
upper_alt = 1000
alt_delim = 5
lower_i = 0
upper_i = 180
i_delim = 1
np_alt = np.linspace(lower_alt, upper_alt, int((upper_alt+alt_delim-lower_alt)/alt_delim))
altitudes = dict.fromkeys(list(np_alt),0)
deb_altitudes = dict.fromkeys(list(np_alt),0)
np_inc = np.linspace(lower_i, upper_i, int((upper_i+i_delim-lower_i)/i_delim))
inclinations = dict.fromkeys(list(np_inc),0) #Range from 0 to 180 deg inclin. in 5 deg steps
deb_inclinations = dict.fromkeys(list(np_inc),0)
z = np.zeros([len(np_inc),len(np_alt)])
deb_z = np.zeros([len(np_inc),len(np_alt)])
x_vals = []
y_vals = []
deb_x_vals = []
deb_y_vals = []
debCounter = 0
for obj in objects:
if obj.deb:
debCounter += 1
inc = math.degrees(obj.i)
alt = (obj.a - Re)/1000 #km
nearest_alt = int(alt//alt_delim)*alt_delim
nearest_inc = int(inc//i_delim)*i_delim
if nearest_alt >= lower_alt and nearest_alt <= upper_alt:
altitudes[nearest_alt]+= 1
if obj.deb:
deb_altitudes[nearest_alt]+=1
if nearest_inc >= lower_i and nearest_inc <= upper_i:
inclinations[nearest_inc]+=1
if obj.deb:
deb_inclinations[nearest_inc]+=1
## We need both within range to generate the image array
if nearest_alt < lower_alt or nearest_alt > upper_alt \
or nearest_inc < lower_i or nearest_inc > upper_i:
continue
x_coord = np.where(np_alt == nearest_alt)[0][0]
y_coord = np.where(np_inc == nearest_inc)[0][0]
z[y_coord][x_coord]+= 1
x_vals.append(x_coord)
y_vals.append(y_coord)
if obj.deb:
deb_z[y_coord][x_coord]+= 1
deb_x_vals.append(x_coord)
deb_y_vals.append(y_coord)
print('{} TLEs, with {} deb'.format(len(objects),debCounter))
Blues9.reverse()
cmap = LogColorMapper(palette=Blues9,low=1,high=z.max())
output_file('Plots/tle_distro.html')
## All TLE plots
p = figure(title='Distribution of TLEs', x_axis_label='Altitude (km)', y_axis_label='TLE Amount',
tooltips=[("x","$x"),("y","$y")])
p.vbar(x=list(altitudes.keys()), top=list(altitudes.values()), bottom=0, width=alt_delim)
p2 = figure(title='Distribution of TLEs', x_axis_label='Inclination (deg)', y_axis_label='TLE Amount',
tooltips=[("x","$x"),("y","$y")])
p2.vbar(x=list(inclinations.keys()), top=list(inclinations.values()), bottom=0, width=i_delim)
p3 = figure(title='Distribution of TLEs',x_axis_label='Altitude (km)',y_axis_label='Inclination (deg)',
x_range=(lower_alt,upper_alt),y_range=(lower_i,upper_i),
tooltips=[("x","$x"),("y","$y"),("value","@image")])
p3.image(image=[z],x=lower_alt,y=lower_i,dw=upper_alt-lower_alt,dh=upper_i-lower_i,color_mapper=cmap)
color_bar = ColorBar(color_mapper=cmap,location=(0,0),ticker=LogTicker())
p3.add_layout(color_bar,'right')
p4 = figure(title='Distribution of TLEs',x_axis_label='Altitude (km)',y_axis_label='Inclination (deg)',
tooltips=[("x","$x"),("y","$y")])
p4.scatter(x=x_vals,y=y_vals,marker='circle')
## Debris plots
p5 = figure(title='Distribution of Debris', x_axis_label='Altitude (km)', y_axis_label='Debris Amount',
tooltips=[("x","$x"),("y","$y")])
p5.vbar(x=list(deb_altitudes.keys()), top=list(deb_altitudes.values()), bottom=0, width=alt_delim)
p6 = figure(title='Distribution of Debris', x_axis_label='Inclination (deg)', y_axis_label='Debris Amount',
tooltips=[("x","$x"),("y","$y")])
p6.vbar(x=list(deb_inclinations.keys()), top=list(deb_inclinations.values()), bottom=0, width=i_delim)
p7 = figure(title='Distribution of Debris', x_axis_label='Altitude (km)', y_axis_label='Inclination (deg)',
x_range=(lower_alt, upper_alt), y_range=(lower_i, upper_i),
tooltips=[("x", "$x"), ("y", "$y"), ("value", "@image")])
p7.image(image=[z], x=lower_alt, y=lower_i, dw=upper_alt - lower_alt, dh=upper_i - lower_i, color_mapper=cmap)
color_bar = ColorBar(color_mapper=cmap, location=(0, 0), ticker=LogTicker())
p7.add_layout(color_bar, 'right')
p8 = figure(title='Distribution of Debris', x_axis_label='Altitude (km)', y_axis_label='Inclination (deg)',
tooltips=[("x", "$x"), ("y", "$y")])
p8.scatter(x=deb_x_vals, y=deb_y_vals, marker='circle')
grid = gridplot([[p, p2],[p3, p4], [p5, p6], [p7, p8]])
show(grid)